Jul 10
(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))
Jul 10
(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))
Jul 9
(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))
Jul 6
(define (sum term a next b)
  (define (iter a result)
    (if (> a b)
        result
        (iter (next a) (+ result (term a)))))
  (iter a 0))

 

Jul 6
(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更准确点

 

May 17

 

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

 

 

 

 

May 14

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

May 13

当输入的数大于2^32-1时,计算时间大量增加,就不符合增长的阶o(log(n))了,而之前得到的值一定小于m,显然在计算大数的时候,速度差异就会很大

May 12

用了一段时间的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的大牛队= =

 

 

May 12

这题没什么难度,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

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