2005/07/15

スターウォーズ エピソード3を見にいったら、コルサント→モンテカルロみたいな連想が働いてしまって、長方形が重なるかどうかを面積の差で解決する問題(参照)ってモンテカルロで解けるじゃん。

(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

これなら長方形を実数で指定できる! ハエを殺すのに焼夷弾を持ち出すみたいな話だってのはわかってますって。

2 comments :

hisashim said...

; 無粋なへそ曲がりより

(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

k16 said...

0.00001の差を判定するためには、10^10 以上の試行が必要になるはずです(あやふやですが、理想的な乱数であれば、モンテカルロ法は n^(-1/2) のオーダーで誤差が収束していくんだったような)。

ちゃんとやるには、与えられた矩形の精度に合わせて試行回数を決定するなどの前処理が必要になるんでしょう。いずれにしろ、この実装では 10^10 回も試行したらフリーズしますが……