SICP学习笔记 1.2.6 实例:素数检测

    练习1.22

;; runtime函数在stk、racket中都不支持，而GNU的Mit-Scheme十分难用，而且在Fedora16上编译安装后启动就报错，
;; 后来总算是在selinux的提示下让它能够正常启动了
;; grep scheme /var/log/audit/audit.log | audit2allow -M mypol
;; semodule -i mypol.pp

 (define (start-prime-test n start-time)
  (and (prime? n)
	  (report-prime n (- (runtime) start-time))))
	  
(define (prime? n)
  (= n (smallest-divisor n)))
  
(define (report-prime n elapsed-time)
  (display n)
  (display " *** ")
  (display elapsed-time)
  (newline))
  
(define (search-for-primes n)
  (if (even? n)
      (search-for-primes-iter (+ n 1) 3 (runtime))
      (search-for-primes-iter n 3 (runtime))))
      
(define (search-for-primes-iter n count start-time)
  (if (= count 0)
      ((newline) (display "search over"))
      (if (start-prime-test n start-time)
          (search-for-primes-iter (+ n 2) (- count 1) (runtime))
          (search-for-primes-iter (+ n 2) count (runtime)))))
          
;; 数据太小的话根本不能进行比较
(search-for-primes 10000000000)
10000000019 *** .23
10000000033 *** .23
10000000061 *** .24
search over

(search-for-primes 100000000000)
100000000003 *** .75
100000000019 *** .73
100000000057 *** .73
search over

(search-for-primes 1000000000000)
1000000000039 *** 2.34
1000000000061 *** 2.33
1000000000063 *** 2.33
search over

;; 耗时大概呈√10倍增长

 

    练习1.23

;; 修改如下
(define (next test)
  (if (= n 2)
      3
      (+ 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 test-divisor)))))
		
(search-for-primes 10000000000)
10000000019 *** .13
10000000033 *** .14
10000000061 *** .14
search over

(search-for-primes 100000000000)
100000000003 *** .45
100000000019 *** .45
100000000057 *** .44
search over

(search-for-primes 1000000000000)
1000000000039 *** 1.41
1000000000061 *** 1.39
1000000000063 *** 1.41
search over

;; 随着数据的增大，耗时大致呈2倍增长

 

    练习1.24

;; 1009 1013 1019
;; 1000003 1000033 1000037

timed-prime-test 1009)
1009 *** .39

(timed-prime-test 1013)
1013 *** .42

(timed-prime-test 1019)
1019 *** .44

(timed-prime-test 1000003)
1000003 *** .77

(timed-prime-test 1000033)
1000033 *** .76

(timed-prime-test 1000037)
1000037 *** .79

;; lg n^2 = 2 * lg n,对比以上结果，比较接近，预计测试数据增大后将更加接近

 

    练习1.25

    理论上可行，但是直接求base^exp的话可能会因为结果太大而出现溢出

 

    练习1.26

    修改后将对于base^2n进行如下方式的计算过程  -->  (base*base*base*...)*(base*base*base*...)

    而原来的方式将进行如下方式的计算过程   -->  (base^n)^2 => (base^(n/2)^2)^2 => ...
    因此原有方式为O(log n)而修改后为O(n)

 

    练习1.27

;; 561 1105 1729 2465 2821 6601

(define (expmod base exp m)
  (cond ((= exp 0) 1)
		((even? exp)
		 (remainder (square (expmod base (/ exp 2) m)) m))
		(else
		 (remainder (* base (expmod base (- exp 1) m)) m))))

(define (carmichael n)
  (carmichael-iter n 1))

(define (carmichael-iter n a)
  (cond ((= a n) true)
		((= (expmod a n n) a) (carmichael-iter n (+ a 1)))
		(else false)))
		
(carmichael 561)
;Value: #t

(carmichael 1105)
;Value: #t

(carmichael 1729)
;Value: #t

(carmichael 2465)
;Value: #t

(carmichael 2821)
;Value: #t

(carmichael 6601)
;Value: #t

 

    练习1.28

费马小定理:
如果n是素数, 其中1<a<n, 则有 a^n ≡ a (mod n)

变形过程为:
已知 a^n ≡ a (mod n)
-->     a^n = k*n + a  
-->     a^(n-1) = (k/a)*n + 1
-->     a^(n-1) = k'*n + 1
-->     a^(n-1) ≡ 1 (mod n)

设n是素数, 其中1<a<n-1, 假设 a^2 ≡ 1 (mod n)
-->     a^2 - 1 = k*n
-->     (a+1)*(a-1) = k*n
-->     (a+1)*(a-1) ≡ 0 (mod n)
-->     a+1 ≡ 0 (mod n)
        a-1 ≡ 0 (mod n)
-->     a = n-1
        a = 1
        而1<a<n-1，则n不为素数，或不存在1<a<n-1, a^2 ≡ 1 (mod n)
       
;; 561 1105 1729 2465 2821 6601

(define (miller-rabin n)
  (miller-rabin-iter n 1))

(define (miller-rabin-iter n a)
  (cond ((= a n) true)
		((= (expmod a (- n 1) n) 1) (miller-rabin-iter n (+ a 1)))
		(else false)))
		
(define (expmod base exp m)
  (cond ((= exp 0) 1)
		((even? exp)  
		 (nontrivial-square-root (expmod base (/ exp 2) m) m))
		(else
		 (remainder (* base (expmod base (- exp 1) m)) m))))
		 
(define (nontrivial-square-root a n)
  (define (try-it value)
    (if (and (> a 1) (< a (- n 1)) (= value 1))
        0
        value))
  (try-it (remainder (square a) n)))
  
(miller-rabin 561)
;Value: #f

(miller-rabin 1105)
;Value: #f

(miller-rabin 1729)
;Value: #f

(miller-rabin 2465)
;Value: #f

(miller-rabin 2821)
;Value: #f

(miller-rabin 6601)
;Value: #f