#lang sicp
;余数
(remainder 9 7)
;求最大公约数
(define (gcd a b)
(if (= b 0)
a
(gcd b (remainder a b))))
; 线性同余法,a 和 m 是素数
(define (rand-update x)
(let ((a 48271) (b 19851020) (m 2147483647))
(modulo (+ (* a x) b) m)))
(define random-init 7)
(define rand
(let ((x random-init))
(lambda ()
(set! x (rand-update x))
x)))
;(rand)
;(rand)
;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
;6/(PI^2) = 两个整数最大公约数为1的概率
;也就是6/概率 = pi^2
;所以pi = (sqrt (/ 6 概率))
;trial试验,试验的次数
;sqrt求平方根
(define (estimate-pi trials)
(sqrt (/ 6 (monte-carlo trials cesaro-test))))
;返回值为true或者是false,判断两个随机数的最大公约数是否为1(互质)
(define (cesaro-test)
(= (gcd (rand) (rand)) 1))
;(cesaro-test)
;experiment 试验,最大公约数是否为1
;蒙特卡罗的参数,试验次数,每次试验的值(即:最大公约数是否为1)
;根据每次试验的值,进行统计,统计次数为试验次数
;返回值为一个概率统计
(define (monte-carlo trials experiment)
(define (iter trials-remaining trials-passed)
(cond ((= trials-remaining 0)
(/ trials-passed trials))
((experiment)
(iter (- trials-remaining 1) (+ trials-passed 1)))
(else
(iter (- trials-remaining 1) trials-passed))))
(iter trials 0))
;;;;;;;;;;;;;;;;;
(define (estimate-pi-2 trials)
(sqrt (/ 6 (random-gcd-test trials random-init))))
;需要将判断随机数的最大公约数,嵌入到程序内部,程序必须显示的去操作随机数
(define (random-gcd-test trials initial-x)
(define (iter trials-remaining trials-passed x)
(let ((x1 (rand-update x)))
(let ((x2 (rand-update x1)))
(cond ((= trials-remaining 0)
(/ trials-passed trials))
((= (gcd x1 x2) 1)
(iter (- trials-remaining 1) (+ trials-passed 1) x2))
(else
(iter (- trials-remaining 1) trials-passed x2))))))
(iter trials 0 initial-x))
;;;;;;;;;;;;;;;;;
(estimate-pi 1000000)
(estimate-pi-2 1000000)
;练习3.5
(define (random-in-range low high)
(let ((range (- high low)))
(+ low (random range))))
; (+ low (* (random range)))))
(random-in-range 10 20)
(define (estimate-integral P x1 x2 y1 y2 trials)
(define (experiment)
(P (random-in-range x1 x2) (random-in-range y1 y2)))
;矩形面积
(let ((rt-area (* (abs (- x1 x2))
(abs (- y1 y2)))))
(* rt-area (monte-carlo trials experiment)))); 矩形面积乘以比率