计算机程序的构造和解释(SICP) 第1章 习题解析 Part2

这次回顾第一章的部分习题。

学习资料：

https://ocw.mit.edu/courses/electrical-engineering-and-computer-science/6-001-structure-and-interpretation-of-computer-programs-spring-2005/index.htm

https://github.com/DeathKing/Learning-SICP

https://mitpress.mit.edu/sites/default/files/sicp/index.html

https://www.bilibili.com/video/BV1Xx41117tr?from=search&seid=14983483066585274454

参考资料：

https://sicp.readthedocs.io/en/latest

1.16(p30)

利用如下事实即可：

$$b^n = (b^{2})^{n/2}= (b^{n/2})^2, n=2k$$

代码如下：

(define (expt-iter b counter product)
  (cond ((= counter 0) product)
        ((even? counter) (expt-iter (* b b) (/ counter 2) product))
        (else (expt-iter b (- counter 1) (* product b)))))

(define (expt b n)
    (expt-iter b n 1))

(define (even? n)
    (= (remainder n 2) 0))

1.17(p31)

原始递归方法：

$$\mathrm{prod}(a, b)=\begin{cases} 0 & b=0\\ \mathrm{prod}(a, b - 1) + a & b>0 \end{cases}$$

快速求法：

$$\mathrm{prod}(a, b)=\begin{cases} 0 & b=0\\ \mathrm{prod}(a, b/2) + \mathrm{prod}(a, b/2) & b=2k, k\ge 1\\ \mathrm{prod}(a, b -1) + a& b=2k-1, k\ge 1 \end{cases}$$

代码如下：

(define (double a)
    (+ a a))

(define (halve n)
    (/ n 2))

(define (even? n)
    (= (remainder n 2) 0))

(define (* a b)
    (cond ((= b 0) 0)
          ((even? b) (double (* a (halve b))))
          (else (+ a (* a (- b 1))))))

1.18(p31)

将前一个过程改写成迭代的形式：

(define (double a)
    (+ a a))

(define (halve n)
    (/ n 2))

(define (even? n)
    (= (remainder n 2) 0))

(define (prod-iter a b product)
    (cond ((= b 0) product)
          ((even? b) (+ product (double (prod-iter a (halve b) 0))))
          (else (prod-iter a (- b 1) (+ product a)))))

(define (* a b)
    (prod-iter a b 0))

1.19(p31)

证明：

$$\begin{aligned} T_{pq}^2 (a, b) &= T_{pq}(bq + aq + ap, bp+aq)\\ &=((bp+aq) q + (bq + aq + ap)q + (bq + aq + ap)p, ( bp+aq)p + (bq + aq + ap)q)\\ &= (b(q^2+2pq) + a(q^2+2pq) + a(p^2 + q^2), b(p^2+ q^2) + a(q^2 + 2pq))\\ &= T_{p'q'}(a,b)\\ p'&= p^2 + q^2\\ q'&= q^2 + 2pq \end{aligned}$$

转换成代码即为：

(define (fib n)
  (fib-iter 1 0 0 1 n))

(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                   (+ (* p p) (* q q)) ; compute p'
                   (+ (* q q) (* 2 (* p q))) ; compute q'
                   (/ count 2)))
        (else (fib-iter (+ (* b q) (* a q) (* a p))
                        (+ (* b p) (* a q))
                        p
                        q
                        (- count 1)))))

; for test
(define (fib1 n)
  (fib1-iter 1 0 n))

(define (fib1-iter a b count)
  (if (= count 0)
      b
      (fib1-iter (+ a b) a (- count 1))))

(define (test f1 f2 n)
  (cond ((< (f1 n) (f2 n)) 
          (display "error")
          (newline))
        ((< (f1 n) (f2 n))
          (display "error")
          (newline))))

(test fib fib1 100)

一点思考：

题目中的变换是如何想到的，斐波那契数列为何能用对数时间复杂度计算？

将题目中递推公式写成矩阵形式可得：

$$\begin{aligned} \left( \begin{matrix} bq + aq + ap\\ bp+aq \end{matrix} \right) &= \left( \begin{matrix} p+q & q\\ q & p \end{matrix} \right) \left( \begin{matrix} a\\ b \end{matrix} \right)\\ &\triangleq A_{pq} \left( \begin{matrix} a\\ b \end{matrix} \right) \end{aligned}$$

特别的：

$$\begin{aligned} A_{pq} ^2 &= \left( \begin{matrix} p+q & q\\ q & p \end{matrix} \right)^2\\ &= \left( \begin{matrix} q^2 +2pq & q^2 +2pq\\ q^2 +2pq & q^2 + pq \end{matrix} \right)\\ &=A_{p'q'}\\ p'&= p^2 + q^2\\ q'&= q^2 + 2pq \end{aligned}$$

所以上式用矩阵形式理解更方便。

不过上述形式还是比较难凑出的，这里给出一个一般的方法，对每个线性递推式，必然存在一个对应的矩阵$A$，递推实际上是计算$A^n$，如果$A$可对角化，即

$$A=P\Lambda P^{-1}$$

那么

$$A^n= P \Lambda^n P^{-1}$$

由于$\Lambda$是对角阵，所以只要计算对角元即可，可以利用快速幂的方法在对数时间复杂度内计算。

1.20(p33)

参考资料：

https://sicp.readthedocs.io/en/latest/chp1/20.html

应用序：remainder调用次数为5次。

(gcd 206 40) ; 1
(gcd 40 6) ; 2
(gcd 6 4) ; 3
(gcd 4 2) ; 4
(gcd 2 0)
2

正则序：

(gcd 206 40)
(gcd 40 (remainder 206 40)) ; t1 = (remainder 206 40) = 6, 1
(if (= t1 0)) ; #f
(gcd t1 (remainder 40 t1)) ; t2 = (remainder 40 t1) = (remainder 40 6) = 4, 2
(if (= t2 0)) ; #f
(gcd t2 (remainder t1 t2)) ; t3 = (remainder t1 t2) = (remainder 6 4) = 2, 3
(if (= t3 0)) ; #f
(gcd t3 (remainder t2 t3)) ; t4 = (remainder t2 t3) = (remainder 4 2) = 0, 4
(if (= t4 0)) ; #t
t3
(remainder t1 t2)
(remainder t1 (remainder 40 t1))
(remainder t1 (remainder 40 (remainder 206 40))) ; 5
(remainder t1 (remainder 40 6)) ; 6
(remainder t1 4)
(remainder (remainder 206 40) 4) ; 7
(remainder 6 4) ; 8 
2

remainder调用次数为8次。

1.21(p35)

(display (smallest-divisor 199))
(newline)
(display (smallest-divisor 1999))
(newline)
(display (smallest-divisor 19999))
(newline)

得到结果：

199
1999
7

1.22(p35)

参考资料：https://www.scheme.com/csug8/system.html#./system:h10

Chez Scheme不包含runtime函数，这里使用如下方式计算时间：

(define (runtime)
    (real-time))

完整的代码如下：

(define (runtime)
    (real-time))

(define (square a)
  (* a a))

(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 (report-prime n)
  (display n)
  (display " is prime.")
  (newline))

(define (report-time elapsed-time)
  (display "time is :")
  (display elapsed-time)
  (newline))

(define (search-for-primes n)
    (define start-time (runtime))
    (search-for-primes-iter (+ n 1) start-time 3))

(define (search-for-primes-iter n start-time num)
    (cond ((= num 0) (report-time (- (runtime) start-time)))
          ((prime? n) 
           (report-prime n)
           (search-for-primes-iter (+ n 1) start-time (- num 1)))
          (else 
           (search-for-primes-iter (+ n 1) start-time num))))

; for test
(search-for-primes 1000000)
(search-for-primes 10000000)
(search-for-primes 100000000)
(search-for-primes 1000000000)
(search-for-primes 10000000000)

结果如下：

1000003 is prime.
1000033 is prime.
1000037 is prime.
time is :6
10000019 is prime.
10000079 is prime.
10000103 is prime.
time is :1
100000007 is prime.
100000037 is prime.
100000039 is prime.
time is :1
1000000007 is prime.
1000000009 is prime.
1000000021 is prime.
time is :2
10000000019 is prime.
10000000033 is prime.
10000000061 is prime.
time is :9

结果并不是按照$\sqrt{10}$成倍增长，原因应该是$n$的数量级还不够大。

1.23(p36)

代码如下：

(define (runtime)
    (real-time))

(define (square a)
  (* a a))

(define (next test-divisor)
    (if (= test-divisor 2)
        3
        (+ test-divisor 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 test-divisor)))))

(define (divides? a b)
  (= (remainder b a) 0))

(define (prime? n)
  (= n (smallest-divisor n)))

;;
(define (report-prime n)
  (display n)
  (display " is prime.")
  (newline))

(define (report-time elapsed-time)
  (display "time is :")
  (display elapsed-time)
  (newline))

(define (search-for-primes n)
    (define start-time (runtime))
    (search-for-primes-iter (+ n 1) start-time 3))

(define (search-for-primes-iter n start-time num)
    (cond ((= num 0) (report-time (- (runtime) start-time)))
          ((prime? n) 
           (report-prime n)
           (search-for-primes-iter (+ n 1) start-time (- num 1)))
          (else 
           (search-for-primes-iter (+ n 1) start-time num))))

(search-for-primes 1000000)
(search-for-primes 10000000)
(search-for-primes 100000000)
(search-for-primes 1000000000)
(search-for-primes 10000000000)

结果：

1000003 is prime.  
1000033 is prime.  
1000037 is prime.  
time is :0
10000019 is prime. 
10000079 is prime. 
10000103 is prime. 
time is :1
100000007 is prime.
100000037 is prime.
100000039 is prime.
time is :4
1000000007 is prime.
1000000009 is prime.
1000000021 is prime.
time is :2
10000000019 is prime.
10000000033 is prime.
10000000061 is prime.
time is :3

结果和预期不同，原因应该是$n$的数量级还不够大。

1.24(p36)

代码按如下：

(define (runtime)
    (real-time))

(define (square a)
  (* a a))

(define true #t)

(define false #f)

;; fast-prime?
(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 (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) true)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (else false)))

;;
(define (report-prime n)
  (display n)
  (display " is prime.")
  (newline))

(define (report-time elapsed-time)
  (display "time is :")
  (display elapsed-time)
  (newline))

(define (search-for-primes n)
    (define start-time (runtime))
    (search-for-primes-iter (+ n 1) start-time 12))

(define (search-for-primes-iter n start-time num)
    (cond ((= num 0) (report-time (- (runtime) start-time)))
          ((fast-prime? n 10) 
           (report-prime n)
           (search-for-primes-iter (+ n 1) start-time (- num 1)))
          (else 
           (search-for-primes-iter (+ n 1) start-time num))))

(search-for-primes 1000000)
(search-for-primes 10000000)
(search-for-primes 100000000)
(search-for-primes 1000000000)
(search-for-primes 10000000000)

实验结果：

1000003 is prime.
1000033 is prime.
1000037 is prime.
1000039 is prime.
1000081 is prime.
1000099 is prime.
1000117 is prime.
1000121 is prime.
1000133 is prime.
1000151 is prime.
1000159 is prime.
1000171 is prime.
time is :1
10000019 is prime.
10000079 is prime.
10000103 is prime.
10000121 is prime.
10000139 is prime.
10000141 is prime.
10000169 is prime.
10000189 is prime.
10000223 is prime.
10000229 is prime.
10000247 is prime.
10000253 is prime.
time is :3
100000007 is prime.
100000037 is prime.
100000039 is prime.
100000049 is prime.
100000073 is prime.
100000081 is prime.
100000123 is prime.
100000127 is prime.
100000193 is prime.
100000213 is prime.
100000217 is prime.
100000223 is prime.
time is :3
1000000007 is prime.
1000000009 is prime.
1000000021 is prime.
1000000033 is prime.
1000000087 is prime.
1000000093 is prime.
1000000097 is prime.
1000000103 is prime.
1000000123 is prime.
1000000181 is prime.
1000000207 is prime.
1000000223 is prime.
time is :2
10000000019 is prime.
10000000033 is prime.
10000000061 is prime.
10000000069 is prime.
10000000097 is prime.
10000000103 is prime.
10000000121 is prime.
10000000141 is prime.
10000000147 is prime.
10000000207 is prime.
10000000259 is prime.
10000000277 is prime.
time is :13

这次增加了$n$，但是和期望的结果仍然有差距，个人认为仍然是数量级的原因。

1.25(p36)

不正确，因为$a^n$可能溢出。

1.26(p36)

(* (expmod base (/ exp 2) m)
   (expmod base (/ exp 2) m))

根据应用序，上述代码会计算两次(expmod base (/ exp 2) m)，所以时间复杂度的关系为

$$T(n) =2T(n/ 2)$$

所以是线性时间复杂度。

1.27(p36)

(define (square a)
  (* a a))

(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 (fermat-test n m)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it m))

(define (Carmichael-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (define (iter m)
    (cond ((= m n) #t)
          ((try-it m) (iter (+ m 1)))
          (else #f)))
  (iter 1))

; for test
(define (test n)
    (display n)
    (display " ")
    (display (Carmichael-test n))
    (newline))

(test 561)
(test 1105)
(test 1729)
(test 2465)
(test 2821)
(test 6601)

测试结果：

561 #t
1105 #t
1729 #t
2465 #t
2821 #t
6601 #t

1.28(p37)

参考资料：

https://sicp.readthedocs.io/en/latest/chp1/28.html

代码如下：

(define (square a)
  (* a a))

(define true #t)

(define false #f)

; 判断是不是1取模n的非平凡平方根
; 不等于1, n - 1, 平方模n不等于1
(define (square-mod-test a n)
    (and (not (= a 1))
         (not (= a (- n 1)))
         (= (remainder (square a) n) 1)))

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

(define (Miller-Rabin-test n)
  (define (try-it a)
    (= (expmod a (- n 1) n) 1))
  (try-it (+ 1 (random (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) true)
        ((Miller-Rabin-test n) (fast-prime? n (- times 1)))
        (else false)))

; for test
(define (report-prime n)
  (display n)
  (display " is prime.")
  (newline))

(define (test n)
    (define times 10)
    (define (test-iter m)
        (cond ((= m n) (display "End.") (newline))
              ((fast-prime? m 10)
               (report-prime m)
               (test-iter (+ m 1)))
              (else (test-iter (+ m 1)))))
    (test-iter 2))

(test 100)

测试结果：

2 is prime.
3 is prime.
5 is prime.
7 is prime.
11 is prime.
13 is prime.
17 is prime.
19 is prime.
23 is prime.
29 is prime.
31 is prime.
37 is prime.
41 is prime.
43 is prime.
47 is prime.
53 is prime.
59 is prime.
61 is prime.
67 is prime.
71 is prime.
73 is prime.
79 is prime.
83 is prime.
89 is prime.
97 is prime.
End.

值的回顾的习题

对数时间复杂度计算斐波那契数列：1.19。

Chez Scheme中如何计算程序运算时间：1.22。

随机算法判断质数：1.24, 1.28.