いまだゴールが見えないSICP。ようやく第3章が終わった。思い返せば第2章を終えたのが今年の3月22日で、それから5ヵ月もかかってる。おそすぎ。まあゆっくりやります。

最後の問題は、また乱数とモンテカルロだった。ところで第3章では疑似乱数を線形合同法で求める方法が紹介されているので（226ページの注6）、Gauche に用意されているメルセンヌ・ツイスターを使っちゃずるいと思う。いや、べつに誰がずるいというわけでなく、それでは自分的に勉強にならないと思う。ってな大口を叩いてたところで、Knuth本を持ってないからWikipediaに載っている係数をひろってくるだけなんだけど。

(define (monte-carlo experiment-stream passed failed)
  (define (next passed failed)
    (stream-cons
     (/ passed (+ passed failed))
     (monte-carlo
      (stream-cdr experiment-stream) passed failed)))
  (if (stream-car experiment-stream)
      (next (+ passed 1) failed)
      (next passed (+ failed 1))))

;; ex. 3.81
(define rand-init 1)
(define (rand-update x)
  (modulo (* x 1664525) 1013904223))
(define (random-numbers c)
  (define (rands init)
    (stream-cons init
                 (stream-map rand-update
                             (stream-delay (rands init)))))
  (cond ((equal? c 'generate)
         (rands rand-init))
        ((equal? c 'reset)
         (lambda (init)
           (rands init)))
        (else
         (error "Unknown argument" c))))

;; ex. 3.82
(define (random-numbers-in-range x1 x2 rand-init)
  (define (normalize x)
    (+ (* x (/ (- x2 x1) 1013904223)) x1))
  (stream-map normalize ((random-numbers 'reset) rand-init)))
(define (estimate-integral P x1 x2 y1 y2)
  (stream-scale 
   (monte-carlo
    (stream-map P 
                (random-numbers-in-range x1 x2 1)
                (random-numbers-in-range y1 y2 506952111))
    0 0)
   (* (- x2 x1) (- y2 y1))))

(define (P x y)
  (<= (+ (square (- x 5))
         (square (- y 7)))
      9))
(define ex-3-82
  (estimate-integral P 2 8 4 10))
(define (ex-3-82-pi howfar)
  (/ (stream-ref ex-3-82 howfar) 9))

gosh > (ex-3-82-pi 1000)
3.1128871128871127

最初、x軸とy軸について同じ random-numbers-in-range を使ってたら、どうにも想定している結果よりかなり小さい値しか求まらない。そりゃあ y=x 上の点しか拾ってないんだから、よく考えれば当り前の話だった。結果的に初期値として選んだ 1 と 506952111 は適当で、あまり背景はない（1013904223の約半分）。これくらいオーダーを違くすれば、そこそこバラけるような。