(use gauche.uvector)
(use srfi-27)
(use math.mt-random)
;; make random vector
(define (make-random-vector! n)
(mt-random-fill-f64vector! (make-random-source) (make-f64vector n)))
(define (ith-random randoms i)
(f64vector-ref randoms i))
(define (ith-random-pair randoms i)
(list (ith-random randoms (* 2 i))
(ith-random randoms (+ (* 2 i) 1))))
;; rect as ((x0 x1) (y0 y1))
(define (point-inner-closed? point closed)
(and (< (car closed) point) (< point (cadr closed))))
(define (pair-inner-rects? pair rects) ; not cps !
(let ((rect (car rects)))
(if (null? (cdr rects))
(and (point-inner-closed? (car pair) (car rect))
(point-inner-closed? (cadr pair) (cadr rect)))
(or (pair-inner-rects? pair (list rect))
(pair-inner-rects? pair (cdr rects))))))
;; monte-carlo
;; area as '(rect1 rect2 ...)
(define (monte-carlo area n)
(let ((randoms (make-random-vector! (* 2 (+ n 1)))))
(let f ((i 0))
(if (> i n)
0
(if (pair-inner-rects? (ith-random-pair randoms i) area)
(+ (f (+ i 1)) 1)
(f (+ i 1)))))))
;; main
(define times 10000)
(define rect1 (list (list 0.0 0.2) (list 0.0 0.2)))
(define rect2 (list (list 0.1 0.3) (list 0.1 0.3)))
(define s1 (monte-carlo (list rect1) times))
(define s2 (monte-carlo (list rect2) times))
(define s1-and-s2 (monte-carlo (list rect1 rect2) times))
(display (< s1-and-s2 (+ s1 s2)))
=> #t
これなら長方形を実数で指定できる! ハエを殺すのに焼夷弾を持ち出すみたいな話だってのはわかってますって。
; 無粋なへそ曲がりより
返信削除(define rect1 (list (list 0.0 0.500005) (list 0.0 0.500005)))
(define rect2 (list (list 0.499995 1.0) (list 0.499995 1.0)))
(define times 1000000)
(display (+ s1 s2))
(display s1-and-s2)
(display (< s1-and-s2 (+ s1 s2)))
;=> 500069500069#f
0.00001の差を判定するためには、10^10 以上の試行が必要になるはずです(あやふやですが、理想的な乱数であれば、モンテカルロ法は n^(-1/2) のオーダーで誤差が収束していくんだったような)。
返信削除ちゃんとやるには、与えられた矩形の精度に合わせて試行回数を決定するなどの前処理が必要になるんでしょう。いずれにしろ、この実装では 10^10 回も試行したらフリーズしますが……