(a)这明显是个疑问句,而不是什么可执行的东东
(b)分母不能为0
(define (filtered-accumulate filter combiner null-value term a next b)
(cond ((> a b) null-value)
((filter a)
(combiner (term a)
(filtered-accumulate filter combiner null-value term (next a) next b)))
(else (filtered-accumulate filter combiner null-value term (next a) next b))))
(define (square x) (* x x))
(define (fast-prime? n)
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder (square (expmod base (/ exp 2) m))
m))
(else
(remainder (* (remainder base m) (expmod base (- exp 1) m))
m))))
(define (fermat-test a)
(= (expmod a n n) a))
(define (test a)
(cond ((= a 0) #t)
((fermat-test a) (fermat-test (- a 1)))
(else #f)))
(test (- n 1)))
(define (indentity n) n)
(define (inc n) (+ n 1))
(define (prime-add a b)
(filtered-accumulate fast-prime? + 0 square a inc b))
(define (gcd a b)
(if (= b 0)
a
(gcd b (remainder a b))))
(define (check-add n)
(define (check m)
(if (= (gcd m n) 1)
#t
#f))
(filtered-accumulate check * 1 indentity 1 inc n))
(define (accumulate combiner null-value term a next b)
(if (> a b)
null-value
(combiner (term a)
(accumulate combiner null-value term (next a) next b))))
(define (sum term a next b)
(accumulate + 0 term a next b))
(define (product term a next b)
(accumulate * 1 term a next b))
(define (accumulate-iter combiner null-value term a next b)
(define (iter a result)
(if (> a b)
result
(iter (next a) (combiner (term a) result))))
(iter a null-value))
(define (sum-iter term a next b)
(accumulate-iter + 0 term a next b))
(define (product-iter term a next b)
(accumulate-iter * 1 term a next b))
(define (product term a next b)
(if (> a b)
1
(* (term a)
(product term (next a) next b))))
(define (inc n) (+ n 1))
(define (indentity x) x)
(define (factoral n)
(product indentity 1 inc n))
(define (pi-up n)
(cond ((= n 1) 2)
((even? n) (+ n 2))
(else (+ n 1))))
(define (pi-down n)
(cond ((even? n ) (+ n 1))
(else (+ n 2))))
(define (sequare x) (* x x))
(define (get-pi n)
(* 1.0 (/ (product pi-up 1 inc n)
(product pi-down 1 inc n))))
(define (pi-all n)
(cond ((= n 1) (/ 2 3))
((even? n) (/ (+ n 2) (+ n 1)))
(else (/ (+ n 1) (+ n 2)))))
(define (fast-pi n)
(* 1.0 (product pi-all 1 inc n)))
(define (product term a next b result)
(if (> a b)
result
(product term (next a) next b (* result (term a)))))
(define (product term a next b)
(define (iter a result)
(iter (next a)
(* (term a) result)))
(iter a 1))
(define (sum term a next b)
(define (iter a result)
(if (> a b)
result
(iter (next a) (+ result (term a)))))
(iter a 0))
(define (cube x) (* x x x))
(define (sum term a next b)
(if (> a b)
0
(+ (term a)
(sum term (next a) next b))))
(define (simpson f a b n)
(define (get-h) (/ (- b a) n))
(define (get-y k) (f (+ a (* k (get-h)))))
(define (simpson-term k)
(cond ((or (= k 0) (= k n)) (get-y k))
((even? k) (* 2.0 (get-y k)))
(else (* 4.0 (get-y k)))))
(define (simpson-next k) (+ k 1))
(* (sum simpson-term 0 simpson-next n) (/ (get-h) 3.0)))
(define (integral f a b dx)
(define (add-dx x) (+ x dx))
(* (sum f (+ a (/ dx 2.0)) add-dx b)
dx))
(simpson cube 0 1 10000)
(integral cube 0 1 0.0001)
输出
0.2500000000000011
0.24999999874993412
simpson更准确点
(define (square x) (* x x))
(define (check-one count a m)
(define (check-1? t)
(cond ((and (not (= a (- m 1))) (= count 1)) t)
((< count 1) t)
(else 0)))
(check-1? (remainder (square a) m)))
(define (check-two count a m)
(define (check-2? t)
(cond ((and (not (= a 1)) (not (= a (- m 1))) (= count 0)) t)
((not (= count 0)) t)
(else 0)))
(check-2? (remainder (square a) m)))
(define (expmod base exp m count)
(cond ((= exp 0) 1)
((and (even? exp) (even? (/ exp 2)))
(check-one count (expmod base (/ exp 2) m count) m))
((and (even? exp) (not (even? (/ exp 2))))
(check-two (- count 1) (expmod base (/ exp 2) m (- count 1)) m))
(else
(remainder (* base (expmod base (- exp 1) m count))
m))))
(define (prime-test m)
(define (try-it a)
(= (expmod a (- m 1) m 1) 0))
(try-it (+ 1 (random (- m 1)))))
(define (miller-test m times)
(cond ((= m 1) #f)
((= m 2) #t)
((even? m) #f)
((= times 0) #t)
((prime-test m) (miller-test m (- times 1)))
(else #f)))
这个是了解miller-test之后写的
http://en.wikipedia.org/wiki/Miller%E2%80%93Rabin_primality_test
然后这个是按照书上的要求写的
(define (square x) (* x x))
(define (check a n)
(define (check-1? t)
(if (and (> a 1)
(< a (- n 1))
(= t 1))
0 t))
(check-1? (remainder (square a) n)))
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(check (expmod base (/ exp 2) m) m))
(else
(remainder (* base (expmod base (- exp 1) m)) m))))
(define (prime-test m)
(define (try-it a)
(= (expmod a (- m 1) m) 1))
(try-it (+ 1 (random (- m 1)))))
(define (miller-test m times)
(cond ((= m 1) #f)
((= times 0) #t)
((prime-test m) (miller-test m (- times 1)))
(else #f)))
当输入的数大于2^32-1时,计算时间大量增加,就不符合增长的阶o(log(n))了,而之前得到的值一定小于m,显然在计算大数的时候,速度差异就会很大
用了一段时间的PLT和MIT,就测试数据方面来讲,MIT的效率可以说是龟爬了,开c语言编译都不用这么久,当然了,这是指大数方面
PLT的(current-milliseconds)远远没有MIT的(runtime)稳定,数据差距很大,但是测试多组数据的时候就很好,MIT适合单组数据
因为不知道PLT下面的random是什么,所以就MIT解决了
(define (fast-prime? n times)
(define (expmod base exp m)
(cond ((= exp 0) 1)
((even? exp)
(remainder (square (expmod base (/ exp 2) m))
m))
(else
(remainder (* (remainder base m) (expmod base (- exp 1) m))
m))))
(define (fermat-test n)
(define (try-it a)
(= (expmod a n n) a))
(try-it (+ 1 (random (- n 1)))))
(cond ((= times 0) true)
((fermat-test n) (fast-prime? n (- times 1)))
(else false)))
(define (start-prime-test n start-time times)
(if (fast-prime? n times)
(- (runtime) start-time)))
(define (time-prime-test n times)
(start-prime-test n (runtime) times))
(time-prime-test 100003 1000)
;Value: .1559999999999917
(time-prime-test 100019 1000)
;Value: .1880000000000024
(time-prime-test 100043 1000)
;Value: .17199999999999704
(time-prime-test 10000000019 1000)
;Value: .3430000000000035
(time-prime-test 10000000033 1000)
;Value: .3269999999999982
(time-prime-test 10000000061 1000)
;Value: .37400000000000944
还是比较接近两倍的,
http://eli.thegreenplace.net/2007/07/09/sicp-section-126/
http://oss.timedia.co.jp/show
这两个都说不是,eli的循环策略很不错,明天可以试试看,放大效果会比较好吧,那个日本人写的测试数据太小了,感觉不靠谱
http://www.cppblog.com/cuigang/archive/2008/03/31/45781.html
http://silvernightsky.ycool.com/post.2239761.html
以上两个是支持两倍的,各种不懂,测试数据大了,MIT就要求debug,具体不清楚
不死心,试了试10^4和10^8
(time-prime-test 10007 1000)
;Value: .125
(time-prime-test 10009 1000)
;Value: .125
(time-prime-test 100000007 1000)
;Value: .25
(time-prime-test 100000037 1000)
;Value: .25
这数据让我很无语,真够整的,但是
(time-prime-test 1000033 1000)
;Value: .15600000000000058
(time-prime-test 1000033 1000)
;Value: .17099999999999937
(time-prime-test 1000033 1000)
;Value: .15600000000000236
(time-prime-test 1000033 1000)
;Value: .1559999999999988
(time-prime-test 1000033 1000)
;Value: .15600000000000236
(time-prime-test 1000033 1000)
;Value: .14099999999999824
(time-prime-test 1000 1000)
;Unspecified return value
(time-prime-test 1009 1000)
;Value: .09299999999999997
(time-prime-test 1013 1000)
;Value: .10900000000000176
很明显10^3和10^6就有问题了,基本是1.5倍,这让我对(runtime)的精度很是怀疑,有个办法很好,去找本scheme的书看看,手册之类的,但是scheme的书真的很少,而且英文居多,这本sicp基本让我头大了,虽然立志看原版书,但是没有大牛的英语水平,基本靠着google dictionary过来的
其实我很怀疑eli和那个日本人的测试,理论上来讲,fermat的测试前面是o(log(n)),但是times必须是一定的,可能他们两个的times跟测试数据有关吧,我也不好猜测什么
明天把,先做个循环,看下测试数据,然后联系下eli,问问他是怎么看的,估计又要麻烦裘老师了,恨自己没考上PKU啊,学校差距真的很大,大到提前招的新生都能秒了HDU的大牛队= =