secje

(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)))

这个之前为了在PLT下面运行已经写过一个了，详情请看another-prime-test

当输入的数大于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的大牛队= =

这题没什么难度，just modify a little bit from the previos procedure (prime? n)

(define (runtime) (current-milliseconds))
(define (square x) (* x x))
(define (next-divisor n)
  (if (= n 2)
      3
      (+ n 2)))
(define (smallest-divisor n)
  (find-divisor n 2))
(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (next-divisor test-divisor)))))
(define (divides? a b)
  (= (remainder b a) 0))
(define (prime? n)
  (= n (smallest-divisor n)))
(define (timed-prime-test n)
  (start-prime-test n (runtime)))
(define (start-prime-test n start-time)
  (if (prime? n)
      (- (runtime) start-time)))
(timed-prime-test 10000000019)
(timed-prime-test 10000000033)
(timed-prime-test 10000000061)

然后输出

594
523
600

previous version

(define (runtime) (current-milliseconds))
(define (square x) (* x x))
(define (smallest-divisor n)
  (find-divisor n 2))
(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (+ test-divisor 1)))))
(define (divides? a b)
  (= (remainder b a) 0))
(define (prime? n)
  (= n (smallest-divisor n)))
(define (timed-prime-test n)
  (start-prime-test n (runtime)))
(define (start-prime-test n start-time)
  (if (prime? n)
      (- (runtime) start-time)))
(timed-prime-test 10000000019)
(timed-prime-test 10000000033)
(timed-prime-test 10000000061)

输出

1507
1029
1203

两倍的关系还是比较相近的