Permalink
Find file
Fetching contributors…
Cannot retrieve contributors at this time
109 lines (90 sloc) 3.96 KB
;;;; $Id: 1979f6357459652341a9079baf09b3672d726bb3 $
;;;;
;;;; Copyright (c) 2010 Steve Knight <stkni@gmail.com>
;;;;
;;;; Permission is hereby granted, free of charge, to any person obtaining
;;;; a copy of this software and associated documentation files (the
;;;; "Software"), to deal in the Software without restriction, including
;;;; without limitation the rights to use, copy, modify, merge, publish,
;;;; distribute, sublicense, and/or sell copies of the Software, and to
;;;; permit persons to whom the Software is furnished to do so, subject to
;;;; the following conditions:
;;;;
;;;; The above copyright notice and this permission notice shall be
;;;; included in all copies or substantial portions of the Software.
;;;;
;;;; THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
;;;; EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
;;;; MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND
;;;; NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE
;;;; LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION
;;;; OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION
;;;; WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
;;;;
(defun make-point (x y)
(cons x y))
(defun point-x (p)
(car p))
(defun point-y (p)
(cdr p))
(defun sqr (x)
(* x x))
(defun hypoteneuse-length (point)
"The square of the hypoteneuse is equal to the sum of the square of
the two adjacent sides."
(sqrt (+ (sqr (point-x point))
(sqr (point-y point)))))
(defun reflect (x r)
"Given a point on an axis reflect it about r"
(+ r (- r x)))
(defun reflect-in-x (point r)
"Reflect x about the line rx"
(make-point (reflect (point-x point) r)
(point-y point)))
(defun reflect-in-y (point r)
"Reflect y about the line rx"
(make-point (point-x point)
(reflect (point-y point) r)))
(defun reflect-in-xy (point r)
(make-point (reflect (point-x point) r)
(reflect (point-y point) r)))
(defun inside-quadrant-p (point radius)
"We assume that there is an arc starting from 0 radians from the x-axis
going to pi/2 radians to meet the top of the y-axis.
Returns T if point is inside the quadrant and NIL otherwise"
(<= (hypoteneuse-length point) radius))
(defun inside-shape-p (point radius)
"The shape is 4 arcs succesively rotated through pi radians.
If we rotate the given point in each axis (by reflection) we know if it
lies inside the shape because it must be inside all 4 arcs."
(let ((r (/ radius 2)))
(and (inside-quadrant-p point radius)
(inside-quadrant-p (reflect-in-x point r) radius)
(inside-quadrant-p (reflect-in-y point r) radius)
(inside-quadrant-p (reflect-in-xy point r) radius))))
(defun monte-carlo (trials ftest farg)
"Execute ftest "
(let ((result (loop for trial from 0 to trials
count (apply ftest (funcall farg trial)))))
(/ result trials)))
(defun monte-carlo-test (radius trials f)
"Given radius run the test for n trials. returning the proportion
of passed tests against trials."
(monte-carlo trials
(lambda (x y)
(funcall f (make-point x y) radius))
(lambda (x)
(declare (ignore x))
(list (random radius) (random radius)))))
;(* 16.0 (coerce (monte-carlo-test 4.0 10000000 #'inside-shape-p) 'float))
(defun draw-ascii-shape (f radius width-in-chars)
"It's pretty helpful to get an ascii representation of the shape to be able
to visually confirm that the inside-shape-p works correctly."
(let ((step (/ radius width-in-chars)))
(loop for y from 0.0 to radius by step
do (format t "~A~%"
(loop for x from 0.0 to radius by step
collect (if (funcall f (make-point x y) radius)
"*"
"."))))))
;(draw-ascii-shape #'inside-shape-p 4 40)