Jade Dungeon

ch01 构造过程抽象 part03

用高阶函数做抽象

过程作为参数

下面三个过程有很大的相似性:

求\(a\)到\(b\)的整数和:

\[ \begin{equation} 1 + 2 + 3 + \cdots \end{equation} \]
(define (sum-integers a b)
  (if (> a b)
      0
      (+ a (sum-integers (+ a 1) b))))

求\(a\)到\(b\)的整数的立方的和:

\[ \begin{equation} 1^3 + 2^3 + 3^3 + \cdots \end{equation} \]
(define (sum-cubes a b)
  (if (> a b)
      0
      (+ (cube a) (sum-cubes (+ a 1) b))))

还有求圆周率的方法,(非常缓慢地)收敛到\(\pi / 8\):

(define (pi-sum a b)
  (if (> a b)
      0
      (+ (/ 1.0 (* a (+ a 2))) (pi-sum (+ a 4) b))))

用数学表示就是:

\[ \begin{equation} \frac{\pi}{8} = \frac{1}{1 \times 3} + \frac{1}{5 \times 7} + \frac{1}{9 \times 11} + \cdots \end{equation} \]

说个题外话,莱布尼茨有更高效的方法求圆周率:

\[ \begin{equation} \frac{\pi}{4} = 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \cdots \end{equation} \]

以上三个方法的主体可以概括为:

(define (<name> a b)
  (if (> a b)
      0
      (+ (<term> a)
			   (<name> (<next> a) b))))

相当于数学中的求和:

\[ \sum_{n=a}^{b}f(n) = f(a) + \cdots + f(b) \]

我们把共用的代码抽象出来,把差异部分作为两个过程参数:

  • term表示每个项的算法。
  • next取下一个数。
(define (sum term a next b)
  (if (> a b)
      0
      (+ (term a)
         (sum term (next a) next b))))

于是,上面的三个方法只要抽象出不同的部分就可以了:

求立方和:

(define (inc n) (+ n 1))    ; 取下一个整数

(define (sum-cubes a b)
  (sum cube a inc b))

;: (sum-cubes 1 10)         ;3025

求整数和:

(define (identity x) x)     ; 每项的值就是整数自身

(define (sum-integers a b)
  (sum identity a inc b))

;: (sum-integers 1 10)     ;55

求PI:

(define (pi-sum a b)
  (define (pi-term x)          ; 计算每一项
    (/ 1.0 (* x (+ x 2))))
  (define (pi-next x)          ; 求下一项的数
    (+ x 4))
  (sum pi-term a pi-next b))

;: (* 8 (pi-sum 1 1000))    ; 3.139592655589783

还可以用更加复杂的例子,求函数\(f\)在范围\([a, b]\)之间定积分的近似值:

\[ \int_{a}^{b} = \big\lceil f\big( a + \frac{dx}{2} \big) + f\big( a + dx + \frac{dx}{2} \big) + f\big( a + 2dx + \frac{dx}{2} \big) + \cdots \big\rceil dx \]
(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2)) add-dx b)
     dx))

(integral cube 0 1 0.01)       ; .24998750000000042
(integral cube 0 1 0.001)      ; .249999875000001

练习 1.29

辛普森公式的积分近似值比上面的方法更加精确:

\[ \frac{h}{3}[ y_0 + 4y_1 + 2y_2 + 4y_3 + 2y_4 + \dotsi + 2y_{n-2} + 4y_{n-1} + y_n ] \]

我们先将题目给出的条件翻译成表达式,将公式主体的分析放到后面。

变量 h 的定义由条件\(h = (b - a)/n\)给出,将它翻译成表达式:

(define h (/ (- b a) n))

变量 a 和 b 是自由的,因为等会 h 就要放到辛普森函数主体内,所以,不用担心这些 自由变量。

继续将下一个条件\(y_k = f(a + kh)\)翻译成函数:

(define (y k)
    (f (+ a (* k h))))

函数\(f\),变量\(a\)和\(h\)都是自由变量。

现在将注意力放到近似值公式上:

\[ \frac{h}{3}[ y_0 + 4y_1 + 2y_2 + 4y_3 + 2y_4 + \dotsi + 2y_{n-2} + 4y_{n-1} + y_n ] \]

在公式的最外边,大括号之外,乘上了\(\frac{h}{3}\),我们同样可以在自己的函数上也用 (* (/ h 3) ...)来达到同样的效果。

在公式的内部,大括号包围的部分,计算多个\(y_k\)之和,根据下标\(k\)的不同,\(y_k\)有 三个不同的乘法因子:

  • 当 \(k = 0\) 或 \(k = n\) 时,乘法因子为 \(1\)
  • 当 \(k\) 为奇数时,乘法因子为 \(4\)
  • 当 \(k\) 为偶数时,乘法因子为 \(2\)

根据这三条规则,可以给出函数factor,它接受一个参数 \(k\) ,返回相应的 \(y_k\) 的乘法因子:

(define (factor k)
    (cond ((or (= k 0) (= k n))  1)
          ((odd? k)              4)
          (else                  2)))

和之前的几个定义一样,factor内部也有一个自由变量 \(n\) 。

有了factor函数,我们就可以使用表达式(* (factor k) (y k))计算出大括号内的各个加法项。

综合以上这些条件,现在可以写出完整的辛普森函数了:

;;; 29-simpson.scm

(load "p38-sum.scm")

(define (simpson f a b n)
    
    (define h (/ (- b a) n))

    (define (y k)
        (f (+ a (* k h))))

    (define (factor k)
        (cond ((or (= k 0) (= k n)) 1)
              ((odd? k)             4)
              (else                 2)))
    
    (define (term k)
        (* (factor k)
           (y k)))

    (define (next k)
        (+ k 1))

    (if (not (even? n))
        (error "n can't be odd")
        (* (/ h 3)
           (sum term (exact->inexact 0) next n))))

sum的定义照抄书本 38 页的代码:

;;; p38-sum.scm

(define (sum term a next b)
    (if (> a b)
        0
        (+ (term a)
           (sum term (next a) next b))))

辛普森函数检查了输入参数n的情况,确保它是一个偶数,否则引发一个错误。

测试

分别将 n 设为 100 和 1000 ,用simpson函数求出 cube 在 0 和 1 的积分:

1 ]=> (load "29-simpson.scm")

;Loading "29-simpson.scm"...
;  Loading "p38-sum.scm"... done
;... done
;Value: simpson

1 ]=> (define (cube x) (* x x x)) ;Value: cube
1 ]=> (simpson cube 0 1 100)      ;Value: .24999999999999992
1 ]=> (simpson cube 0 1 1000)     ;Value: .2500000000000003

可以看到,和书本 39 页integral函数计算出的积分相比,simpson的计算结果精度 更高。

练习 1.30

sum改为迭代版本:

;;; 30-iter-sum.scm

(define (sum term a next b)
    (define (iter a result)
        (if (> a b)
            result
            (iter (next a)
                  (+ (term a) result))))
    (iter a 0))


测试:

1 ]=> (load "30-iter-sum.scm")

;Loading "30-iter-sum.scm"... done
;Value: sum

1 ]=> (sum (lambda (x) x)
           1
           (lambda (i) (+ 1 i))
           10)

;Value: 55

练习 1.31

定义product计算给定范围中各点的某个函数值的乘积。

product(递归版本)

product计算给定范围中各点的某个函数值的乘积,它的递归版本和书本 38 页的 递归版本 sum 非常相似:

;;; 31-rec-product.scm

(define (product term a next b)
    (if (> a b)
        1
        (* (term a)
           (product term (next a) next b))))

测试:

1 ]=> (load "31-rec-product.scm")

;Loading "31-rec-product.scm"... done
;Value: product

1 ]=> (product (lambda (x) x)
               1
               (lambda (i) (+ i 1))
               10)

;Value: 3628800
product(迭代版本)

迭代版本的 product 也和 练习 1.30 的迭代版本 sum 非常相似:

;;; 31-iter-product.scm

(define (product term a next b)
    (define (iter a result)
        (if (> a b)
            result
            (iter (+ a 1)
                  (* (term a) result))))
    (iter a 1))

测试:

1 ]=> (load "31-iter-product.scm")

;Loading "31-iter-product.scm"... done
;Value: product

1 ]=> (product (lambda (x) x)
               1
               (lambda (i) (+ i 1))
               10)

;Value: 3628800
使用 product 定义 factorial

factorial 也就是阶乘,书本 21 页介绍过阶乘的概念。

以下是使用 product 重定义的 factorial :

;;; 31-factorial.scm

(load "31-iter-product.scm")

(define (factorial n)
    (product (lambda (x) x)
             1
             (lambda (i) (+ i 1))
             n))

测试:

1 ]=> (load "31-factorial.scm")

;Loading "31-factorial.scm"...
;  Loading "31-iter-product.scm"... done
;... done
;Value: factorial

1 ]=> (factorial 10)

;Value: 3628800
计算 pi 的近似值

可以将题目给出的公式:

\[ \begin{equation} \begin{split} \frac{\pi}{4} = \frac{ 2 \cdot 4 \cdot 4 \cdot 6 \cdot 6 \cdot 8 \dotsi}{ 3 \cdot 3 \cdot 5 \cdot 5 \cdot 7 \cdot 7 \dotsi} \end{split} \end{equation} \]

转换成:

\[ \begin{equation} \begin{split} \pi = 4 \cdot \bigg (\frac{2 \cdot 4 \cdot 4 \cdot 6 \cdot 6 \cdot 8 \dotsi}{3 \cdot 3 \cdot 5 \cdot 5 \cdot 7 \cdot 7 \dotsi} \bigg) \end{split} \end{equation} \]

很明显,公式括号里面的分子和分母,可以分别用两个乘法序列来表示。

分子的乘法序列\(2, 4, 4, 6, 6, 8, \dotsi\)可以用函数numer-term生成:

;;; 31-numer-term.scm

(define (numer-term i)
    (cond ((= i 1)
            2)
          ((even? i)
            (+ i 2))
          (else
            (+ i 1))))

分母的乘法序列\(3, 3, 5, 5, 7, 7, \dotsi\)可以用函数denom-term生成:

;;; 31-denom-term.scm

(define (denom-term i)
    (if (odd? i)
        (+ i 2)
        (+ i 1)))

组合起productnumer-termdenom-term,完整的\(\pi\)值计算程序的定义如下:

;;; 31-pi.scm

(load "31-iter-product.scm")
(load "31-numer-term.scm")
(load "31-denom-term.scm")

(define (pi n)
    (* 4
        (exact->inexact
            (/ (product numer-term
                        1
                        (lambda (i) (+ i 1))
                        n)
               (product denom-term 
                        1
                        (lambda (i) (+ i 1))
                        n)))))

程序使用了exact->inexact 函数转换除法的商,确保计算所得的结果为浮点数格式(而不是份数格式)。

最后,对程序进行测试:

1 ]=> (load "31-pi.scm")
;Loading "31-pi.scm"...
;  Loading "31-iter-product.scm"... done
;  Loading "31-numer-term.scm"... done
;  Loading "31-denom-term.scm"... done
;... done
;Value: pi

1 ]=> (pi 1)
;Value: 2.6666666666666665

1 ]=> (pi 10)
;Value: 3.2751010413348074

1 ]=> (pi 100)
;Value: 3.1570301764551676

1 ]=> (pi 1000)
;Value: 3.1431607055322663

1 ]=> (pi 10000)
;Value: 3.1417497057380523

1 ]=> (pi 100000)
;Value: 3.1416083612781764

增大和减少输入参数n可以控制计算出的\(\pi\)值的精度。

练习 1.32

sumproduct,我们可以进一步抽象出「累积」的概念,并将它写成相应的 accumulate函数,它和之前的sumproduct函数都非常相似:

;;; 32-rec-accumulate.scm

(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))))
重定义 sum 和 product

接着,使用accumulate函数重定义sum

;;; 32-sum.scm
(load "32-rec-accumulate.scm")

(define (sum term a next b)
	(accumulate + 0 term a next b))

还有product函数:

;;; 32-product.scm
(load "32-rec-accumulate.scm")

(define (product term a next b)
	(accumulate * 1 term a next b))

测试新的sum

1 ]=> (load "32-sum.scm")

;Loading "32-sum.scm"...
;  Loading "32-rec-accumulate.scm"... done
;... done
;Value: sum

1 ]=> (sum (lambda (x) x) 1 (lambda (i) (+ i 1)) 10)

;Value: 55

测试新的product

1 ]=> (load "32-product.scm")

;Loading "32-product.scm"...
;  Loading "32-rec-accumulate.scm"... done
;... done
;Value: product

1 ]=> (product (lambda (x) x) 1 (lambda (i) (+ i 1)) 10)

;Value: 3628800
迭代版 accumulate

迭代版的accumulate和之前迭代版的sumproduct都非常相似:

;;; 32-iter-accumulate.scm

(define (accumulate combiner null-value term a next b)
	(define (iter a result)
		(if (> a b) result
			(iter (next a) (combiner result (term a)))))
	(iter a null-value))

测试:

1 ]=> (load "32-iter-accumulate.scm")

;Loading "32-iter-accumulate.scm"... done
;Value: accumulate

1 ]=> (accumulate +
				  0
				  (lambda (x) x)
				  1
				  (lambda (i) (+ i 1))
				  10)

;Value: 55

练习 1.33

结合过滤器的概念,生成一个更加泛用的过程filtered-accumulate。能够过滤不在 范围内的项。

函数比 练习 1.33 的accumulate函数多增加一个谓词函数参数valid?,用于确保在 计算时只对符合条件的项进行累积:

;;; 33-filtered-accumulate.scm

(define (filtered-accumulate combine null-value term a next b valid?)
	(if (> a b)
		null-value
		(let ((rest-terms (filtered-accumulate 
									 combine null-value term (next a) next b valid?)))
			(if (valid? a)
				(combine (term a) rest-terms)
				rest-terms))))
a) 计算素数之和

使用filtered-accumulate定义primes-sum,它加起给定范围内的所有素数,其中 素数检测函数prime?来自书本 33 页:

;;; 33-primes-sum.scm
(load "33-filtered-accumulate.scm")
(load "p33-prime.scm")

(define (primes-sum a b)
	(filtered-accumulate 
		+ 0 (lambda (x) x)
		a (lambda (i) (+ i 1)) b prime?))

测试:

1 ]=> (load "33-primes-sum.scm")
;Loading "33-primes-sum.scm"...
;  Loading "33-filtered-accumulate.scm"... done
;  Loading "p33-prime.scm"...
;    Loading "p33-smallest-divisor.scm"...
;      Loading "p33-divides.scm"... done
;      Loading "p33-find-divisor.scm"... done
;    ... done
;  ... done
;... done
;Value: primes-sum

1 ]=> (primes-sum 1 10)     ; 1 + 2 + 3 + 5 + 7 = 18
;Value: 18
b) 计算互素正整数之乘积

根据题目给出的互素性质,使用gcd函数写出coprime?谓词,用于检查两个整数是否 互素:

;;; 33-coprime.scm

(define (coprime? i n)
	(and (< i n)
		(= 1 (gcd i n))))

测试:

1 ]=> (load "33-coprime.scm")
;Loading "33-coprime.scm"... done
;Value: coprime?

1 ]=> (coprime? 2 7)
;Value: #t

1 ]=> (coprime? 2 4)
;Value: #f

然后将这个coprime?filtered-accumulate组合起来,写成product-of-coprimes

;;; 33-product-of-coprimes.scm

(load "33-coprime.scm")
(load "33-filtered-accumulate.scm")

(define (product-of-coprimes n)
	(filtered-accumulate 
		* 1 (lambda (x) x) 1 (lambda (i) (+ i 1))
		n (lambda (x) (coprime? x n))))

注意,因为coprime?函数需要两个参数,所以在accumulate函数的最后部分,使用了 一个lambda表达式将参数n闭包进去,作为coprime?的第二个参数。

测试:

1 ]=> (load "33-product-of-coprimes.scm")
;Loading "33-product-of-coprimes.scm"...
;  Loading "33-filtered-accumulate.scm"... done
;  Loading "33-coprime.scm"... done
;... done
;Value: product-of-coprimes

1 ]=> (product-of-coprimes 10)      ; 1 * 3 * 7 * 9 = 189
;Value: 189
迭代版 filtered-accumulate
;;; 33-iter-filtered-accumulate.scm

(define (filtered-accumulate combine null-value term a next b valid?)
	(define (iter i result)
		(cond 
			((> i b) result)
			((valid? i) (iter (next i) (combine (term i) result)))
			(else (iter (next i) result))))
	(iter a null-value))

测试:

1 ]=> (load "33-iter-filtered-accumulate.scm")
;Loading "33-iter-filtered-accumulate.scm"... done
;Value: filtered-accumulate

1 ]=> (filtered-accumulate +                ; 2 + 4 + 6 + 8 + 10 = 30
               0
               (lambda (x) x)
               1
               (lambda (i) (+ i 1))
               10
               even?)

;Value: 30

用lambda构造过程

定义lambda

定义格式与define过程一样,只不过没有过程名称:

(lambda (<formal-parameters>) <body>)

例:

(lambda (a b) (+ a b))

过程定义就是比lambda多一个名称关联,比如:

(define (plus4 x) (+ x 4))

等价于:

(define plus4 (lambda (x) (+ x 4)))

调用lambda

调用,exps相当于实参列表:

((lambda (<formal-parameters>) <body>) <exps>)

例:

((lambda (a b) (+ a b)) 1 2)     ; a和b实参数为1和2,Value = 3

let限定变量作用域

还有let可以定义一系列变量,并把作用域限制在<body>内:

(let ((<var1> <exp1>)
      (<var2> <exp2>)
      (<var3> <exp3>)
      ...
      (<varn> <expn>))
  <body>)

其实这和lambda等价:

((lambda (<var1> ... <varn>) <body>)
	<exp1> <exp2> ...  <expn>)

例:

((lambda (a b) (+ a b)) 1 2)     ; Value = 3

相当于:

(let ((a 1)
      (b 2))
  (+ a b))

因为lambda表达式也可以理解为一个只有指定作用域中有效和过程。let表达式只是 基础lambda表达式的语法外衣,在写法上形参与实参一一对应。

使用let的优点

  • let在尽可能接近使用的地方建立局部变量约束。比如当\(x = 5\)时,下面第二行的 \(x=3\)而第三行的\(x=5\)页:
(+ (let ((x 3))       ; x = 3
      (+ x (* x 10))) ; (+ 3 (* 3 10)) = 33
  x)                  ; x = 5
                      ; (+ 33 5)    = 38
  • 变量的值是在let之外计算的。比如当\(x=2\)时,下面第二行的\(x=2\)而第三行\(x=3\):
(let ((x 3)           ; x = 3
      (y (+ x 2)))    ; y = 2 + 2 = 4
  (* x y))            ; (* 3 4) = 12

如果要使用define达到同样的效果:

(define (f x y)
  (define a (+ 1 (* x y)))
  (define b (- 1 y))
  (+ (* x (square a))
	   (* y b)
		 (* a b)))

相比之下,用let定义变量更加简洁,而define更适合定义内部方法。

练习 1.34

对于定义

;;; 34-f.scm

(define (f g) (g 2))

执行表达式(f f)会造成错误:

1 ]=> (load "34-f.scm")
;Loading "34-f.scm"... done
;Value: f

1 ]=> (f f)
;The object 2 is not applicable.
;To continue, call RESTART with an option number:
; (RESTART 2) => Specify a procedure to use in its place.
; (RESTART 1) => Return to read-eval-print level 1.

要了解出错的原因,我们展开(f f)的执行过程:

(f f)

(f (lambda (g) (g 2)))

((lambda (g) (g 2))
	(lambda (g) (g 2)))

((lambda (g) (g 2))
	2)

(2 2)

执行到(2 2)这一步时,解释器试图以2作为函数进行调用,但是2并不是一个函数, 所以执行出错并打印信息: The object 2 is not applicable 。

过程作为一般性的方法

以过程作为参数的过程可以抽象出多个不同逻辑中的共同部分。以下两个例子都有着相似 的计算过程:

  • 猜测一个值作为参数。
  • 检查是否达到要求:
    • 达到要求就作为返回值。
    • 没有达到要求就改进这个值再作为参数。

通过区间折半方法求方程的平方根

  1. 如果\(f(a) \lt 0 \lt f(b)\)那么\(f()\)在\(a\)与\(b\)之间必然有一个零点。
  2. 取\(a\)与\(b\)的平均值\(x\),看\(f(x)\)是否大于0。
  3. 重复直到误差(\(T\))足够小。

如果初始长度是\(L\),那么需要的步数的增长的为\(\Theta(\log(L/T))\)。

实现方法:

定义过程search,参数为函数f以及使f的值为正与负的两个点。先求出两点的 中点,再检查区间是否足够小:

  • 足够小就返回。
  • 不够小,再取f在这个点的中值与0点,继续判断。
(define (search f neg-point pos-point)
	(let (
			(midpoint (average neg-point pos-point)))
	  (if (close-enough? neg-point pos-point)
			midpoint
	  	(let (
					(test-value (f midpoint)))
	  	  (cond (
					(positive? test-value) (search f neg-point midpoint))
	  	  	((negative? test-value) (search f midpoint pos-point))
	  	  	(else midpoint))))))

判断两个端点足够近的方法:

(define (close-enough? x y)
	(< (abs (- x y)) 0.001))

因为调用search方法时给出的两个点不一定正好是正负两端有值,所以要再加一个 half-interval-method过程检查是否是一正一负两个值,不符合就抛错:

(define (half-interval-method f a b)
	(let (
			(a-value (f a))
			(b-value (f b)))
		(cond 
			((and (negative? a-value) (positive? b-value))
				(search f a b))
			((and (negative? b-value) (positive? a-value))
				(search f b a))
			(else
				(error "Values are not of opposite sign" a b)))))

例:求\(\pi\)近似值,它正好是在\(\sin x = 0\)在2和4之间的根:

(half-interval-method sin 2.0 4.0)

例:求\(x^3 - 2x -3 = 0\)在1和2之间的根:

(half-interval-method 
	(lambda (x) (- (* x x x) (* 2 x) 3)) 1.0 2.0)

找出函数的不动点

如果\(x\)满足\(f(x)=x\),那么\(x\)就是方程\(f()\)的不动点。对于某些函数,可以通过重复 调用\(f(x)\)的方法:

\[ \begin{equation} \begin{split} f(x) & \\ f(f(x)) & \\ f(f(f(x))) & \\ f(f(f(f(x)))) & \\ \cdots{} \end{split} \end{equation} \]

赶到值的变化不大时,可以找到一个不动点。以这个算法实现过程fixed-point

(define tolerance 0.00001)

(define (fixed-point f first-guess)
	(define (close-enough? v1 v2)
		(< (abs (- v1 v2)) tolerance))
	(define (try guess)
		(let (
				(next (f guess)))
			(if (close-enough? guess next)
				next
				(try next))))
	(try first-guess))

例:以1为初始值求余弦函数的不动点

(fixed-point cos 1.0)

例:求\(y=\sin y + \cos y\)的一个解:

(fixed-point (lambda (y) (+ (sin y) (cos y))) 1.0)
lambda的数学符号表示

可以看到求不动点和求平方根根很像:都是从一个值开始不断猜测接近结果。比如求\(x\) 的平方根就是要找到一个\(y\)使得\(y=x/y\):

  • 用程序可以表示为:(lambda (y) (/ x y))
  • 用数学符号表示为:\(y \mapsto x/y\)。

但是用不动点的方式来求平方根的问题在于它是不收敛的,答案在两边振荡:

  1. 以\(y_1\)开始
  2. 下一个数是\(y_2=x/y_1\)
  3. 下一个数是\(y_3=x/y_2=x/(x/y_1)=y_1\)。进入无限循环。
平均阻尼技术

因为答案是在\(y\)和\(x/y\)之间的,所以为了把振荡的范围从\(x/y\)减小,可以用\(y\)与\(x/y\) 的平均值。这样\(y\)的下一个猜测为\((1/2)(y+x/y)\)而不是\(x/y\)。原理是:

\[ \begin{equation} \begin{split} y &= x/y \\ y + y &= x/y + y & \quad \quad \text{两边加上$y$} \\ y &= (1/2)(x/y + y) & \quad \quad \text{两边除以$2$} \end{split} \end{equation} \]

这种取逼近一个解的一系列值的平均值的方法,被称为「平均阻尼技术」。常常在搜寻 不动点的过程中作为收敛手段。

这样计算过程就成了求\(y \mapsto (1/2)(y+x/y)\)的不动点:

(define (sqrt x)
	(fixed-point
		(lambda (y) (average y (/ x y)))
		1.0))

练习 1.35

证明黄金分割率\(\phi\)是变换\(x \mapsto 1 + 1/x\)的不动点。 按定义计算出黄金分割率的值:

;;; 35-golden-ratio.scm
(load "p46-fixed-point.scm")

(define golden-ratio
	(fixed-point
		(lambda (x) (+ 1 (/ 1 x)))
		1.0))

并通过fixed-point算出\(\phi\)的值。测试:

1 ]=> (load "35-golden-ratio.scm")
;Loading "35-golden-ratio.scm"...
;  Loading "p46-fixed-point.scm"... done
;... done
;Value: golden-ratio

1 ]=> golden-ratio
;Value: 1.6180327868852458

练习 1.36

首先,用练习11.22中的newlinedisplay过程把46页的fixed-point函数修改为 一个打印每一步中的猜测值:

;;; 36-fixed-point.scm
(define tolerance 0.000001)

(define (display-info guess step)
	(display "Step: ")
	(display step)
	(display " ")
	
	(display "Guess: ")
	(display guess)
	(newline))

(define (fixed-point f first-guess)
          
	(define (close-enough? v1 v2)
		(< (abs (- v1 v2)) tolerance))

	(define (try guess step)
		(display-info guess step)             ; 每次进入测试时打印一次猜测
		(let (
				(next (f guess)))
			(if (close-enough? next guess)
				(begin                            ; 如果猜测完成
					(display-info next (+ 1 step))  ; 记得算上最后一次计算 next 的猜测
						next)
					(try next (+ 1 step)))))

	(try first-guess 1))

为了在 if 形式中执行多条表达式,fixed-point函数内部使用了begin形式。

函数display-info用于打印猜测值guessdisplay-info的另一个参数step用于 进行步数计算。

然后要求找出\(x \mapsto \log(1000)/\log(x)\)的不动点。

除了不动点函数外,我们还需要一个平均阻尼函数(在书本 48 页定义):

;;; p48-average-damp.scm
(load "p15-average.scm")

(define (average-damp f)
	(lambda (x) (average x (f x))))

最后,将给定的函数映射转换成前序表达式:

;;; 36-formula.scm

(define formula 
    (lambda (x) (/ (log 1000) (log x))))

现在,可以进行求值测试了,先来试试不带平均阻尼的计算:

1 ]=> (load "36-fixed-point.scm")
;Loading "36-fixed-point.scm"... done
;Value: display-info

1 ]=> (load "36-formula.scm")
;Loading "36-formula.scm"... done
;Value: formula

1 ]=> (fixed-point formula 2.0)
Step: 1 Guess: 2.
Step: 2 Guess: 9.965784284662087
Step: 3 Guess: 3.004472209841214
Step: 4 Guess: 6.279195757507157
Step: 5 Guess: 3.759850702401539
Step: 6 Guess: 5.215843784925895
Step: 7 Guess: 4.182207192401397
Step: 8 Guess: 4.8277650983445906
Step: 9 Guess: 4.387593384662677
Step: 10 Guess: 4.671250085763899
Step: 11 Guess: 4.481403616895052
Step: 12 Guess: 4.6053657460929
Step: 13 Guess: 4.5230849678718865
Step: 14 Guess: 4.577114682047341
Step: 15 Guess: 4.541382480151454
Step: 16 Guess: 4.564903245230833
Step: 17 Guess: 4.549372679303342
Step: 18 Guess: 4.559606491913287
Step: 19 Guess: 4.552853875788271
Step: 20 Guess: 4.557305529748263
Step: 21 Guess: 4.554369064436181
Step: 22 Guess: 4.556305311532999
Step: 23 Guess: 4.555028263573554
Step: 24 Guess: 4.555870396702851
Step: 25 Guess: 4.555315001192079
Step: 26 Guess: 4.5556812635433275
Step: 27 Guess: 4.555439715736846
Step: 28 Guess: 4.555599009998291
Step: 29 Guess: 4.555493957531389
Step: 30 Guess: 4.555563237292884
Step: 31 Guess: 4.555517548417651
Step: 32 Guess: 4.555547679306398
Step: 33 Guess: 4.555527808516254
Step: 34 Guess: 4.555540912917957
Step: 35 Guess: 4.555532270803653
Step: 36 Guess: 4.555537970114198
Step: 37 Guess: 4.555534211524127
Step: 38 Guess: 4.555536690243655
Step: 39 Guess: 4.555535055574168
Step: 40 Guess: 4.5555361336081
Step: 41 Guess: 4.555535422664798
;Value: 4.555535422664798

接着,试试使用平均阻尼(别忘了载入 48 页的average-damp函数):

1 ]=> (load "p48-average-damp.scm")
;Loading "p48-average-damp.scm"...
;  Loading "p15-average.scm"... done
;... done
;Value: average-damp

1 ]=> (fixed-point (average-damp formula) 2.0)
Step: 1 Guess: 2.
Step: 2 Guess: 5.9828921423310435
Step: 3 Guess: 4.922168721308343
Step: 4 Guess: 4.628224318195455
Step: 5 Guess: 4.568346513136242
Step: 6 Guess: 4.5577305909237005
Step: 7 Guess: 4.555909809045131
Step: 8 Guess: 4.555599411610624
Step: 9 Guess: 4.5555465521473675
Step: 10 Guess: 4.555537551999825
Step: 11 Guess: 4.555536019631145
Step: 12 Guess: 4.555535758730802
;Value: 4.555535758730802

对比发现,不带平均阻尼的计算使用了 41 步,另一方面,使用平均阻尼的计算只使用了 12 步,说明平均阻尼有助于计算快速收敛。

练习 1.37

如下形式的表达式叫作「无穷连分式」:

\[ f = \cfrac{N_1}{D_1 + \cfrac{N_2}{D_2 + \cfrac{N_3}{D_3 + \cdots}}} \]

当所有的\(N_i\)与\(D_i\)都为1时,值为\(1/\phi\)。如果只取前\(k\)项,那么就称为:\(k\)项 有限边分式,形式为:

\[ \cfrac{N_1}{D_1 + \cfrac{N_2}{D_2 + \cfrac{N_3}{\ddots + \cfrac{N_k}{D_k}}}} \]

过程nd都是只有一个参数i的过程,返回每一项的\(N_i\)与\(D_i\)。定义过程 (cont-frac n d k)求值。查检以下过程是否逼近\(1/\phi\):

(cont-frac
	(lambda (i) 1.0)
	(lambda (i) 1.0)
	k)

k为多大时能有十进制的4位精度?

因为连分式本质上就是一个除法计算序列,所以题目给出\(k\)项连分式:

\[ \begin{equation} \begin{split} \cfrac{N_1}{D_1 + \cfrac{N_2}{\ddots + \cfrac{N_k}{D_k}}} \end{split} \end{equation} \]

可以转换成以下等价的除法计算序列:

\[ \begin{equation} \begin{split} N_1 / (D_1 + (N_2 / (D_2 + \dotsi + (N_k / D_k)))) \end{split} \end{equation} \]

而这个除法序列又可以用一个递归表达式来表示:

\[ \begin{equation} \begin{split} cf(1) \\ N_1 / (D_1 + cf(2)) \\ N_1 / (D_1 + (N_2 / (D_2 + cf(3)))) \\ N_1 / (D_1 + (N_2 / (D_2 + (N_3 / (D_3 + cf(4)))))) \\ \dotsi \\ N_1 / (D_1 + (N_2 / (D_2 + (N_3 / (D_3 + \dotsi + (N_k / D_k)))))) \end{split} \end{equation} \]

其中函数\(cf(i)\)表示连分式的第\(i\)个项。根据这个递归表达式,我们可以给出(递归 计算版本)连分式过程的定义:

;;; 37-rec-cont-frac.scm

(define (cont-frac N D k)

	(define (cf i)
		(if (= k i)
			(/ (N k) (D k))
			(/ (N i) (+ (D i) (cf (+ i 1))))))
	
	(cf 1))
迭代版本的连分式过程

前面说过,一个\(k\)项连分式等价于除法计算序列:

\[ \begin{equation} \begin{split} N_1 / (D_1 + (N_2 / (D_2 + \dotsi + (N_k / D_k)))) \end{split} \end{equation} \]

如果要迭代地计算这个除法计算序列,我们必须倒转公式中各个除法项计算的先后顺序, 先计算高次项,然后才是低次项。也即是说,我们必须先计算:

\[ \begin{equation} \begin{split} cf(k) = N_k / D_k \end{split} \end{equation} \]

然后计算:

\[ \begin{equation} \begin{split} cf(k-1) = N_{k-1} / (D_{k-1} + cf(k)) \end{split} \end{equation} \]

再然后计算:

\[ \begin{equation} \begin{split} cf(k-2) = N_{k-2} / (D_{k-2} + cf(k-1)) \end{split} \end{equation} \]

一直这样反方向回溯下去,直到到达:

\[ \begin{equation} \begin{split} N_1 / (D_1 + cf(2)) \end{split} \end{equation} \]

这时我们就得出了整个连分式的解,而且整个计算过程是迭代地进行的。

这个迭代计算连分式过程的定义如下:

;;; 37-iter-cont-frac.scm

(define (cont-frac N D k)

	(define (iter i result)
		(if (= i 0)
			result
			(iter (- i 1)
				(/ (N i)
					(+ (D i) result)))))

	(iter (- k 1)
		(/ (N k) (D k))))
连分式定义的黄金分割率函数

使用连分式定义的黄金分割率函数的定义如下:

;;; 37-golden-ratio.scm

(load "37-rec-cont-frac.scm")

(define (golden-ratio k)
	(+ 1 (cont-frac 
					(lambda (i) 1.0)
					(lambda (i) 1.0) k)))

测试:

1 ]=> (load "37-golden-ratio.scm")
;Loading "37-golden-ratio.scm"...
;  Loading "37-rec-cont-frac.scm"... done
;... done
;Value: golden-ratio

1 ]=> (golden-ratio 1)
;Value: 2.

1 ]=> (golden-ratio 10)
;Value: 1.6179775280898876

1 ]=> (golden-ratio 11)
;Value: 1.6180555555555556

从测试结果可以看出,只要\(k\)的值大于等于\(11\),就可以保证计算所得的黄金分割率的 精度到达前四位:\(1.618\)。

测试迭代连分式过程

前面的黄金分割率函数使用的是递归版本的连分式过程,现在来试试迭代版本的连分式 程序,确保它运作良好:

;;; 37-golden-ratio-using-iter-cont-frac.scm
(load "37-iter-cont-frac.scm")

(define (golden-ratio k)
	(+ 1
		(cont-frac 
			(lambda (i) 1.0)
			(lambda (i) 1.0)
			k)))

测试:

1 ]=> (load "37-golden-ratio-using-iter-cont-frac.scm")
;Loading "37-golden-ratio-using-iter-cont-frac.scm"...
;  Loading "37-iter-cont-frac.scm"... done
;... done
;Value: golden-ratio

1 ]=> (golden-ratio 1)
;Value: 2.

1 ]=> (golden-ratio 10)
;Value: 1.6179775280898876

1 ]=> (golden-ratio 11)
;Value: 1.6180555555555556

练习 1.38

欧拉关于自然对数底\(e\)提出的连分展开式\(e-2\),在这个式子中\(N_i\)永远\(1\),而\(D_i\) 依次为:

i 1 2 3 4 5 6 7 8 9 10 11 ...
\(D_i\) 1 2 1 1 4 1 1 6 1 1 8 ...

使用练习1.37中的cont-frac过程基于欧拉展开式求出\(e\)的近似值。

观察以上序列,可以发现函数 Di 的规律:

  • 当\((i+1)\)取模 3 等于 0 时,\(D_i\)等于\((i+1)/3 \times 2\)
  • 其他情况下,\(D_i\)返回 1

根据以上规律,可以写出完整的求\(e\)函数了:

;;; 38-e.scm
(load "37-iter-cont-frac.scm")

(define (e k)

	(define (N i) 1)

	(define (D i)
		(if (= 0 (remainder (+ i 1) 3))
			(* 2 (/ (+ i 1) 3))
			1))

	(+ 2.0 (cont-frac N D k)))

测试:

1 ]=> (load "38-e.scm")
;Loading "38-e.scm"...
;  Loading "37-iter-cont-frac.scm"... done
;... done
;Value: e

1 ]=> (e 1)
;Value: 3.

1 ]=> (e 2)
;Value: 2.6666666666666665

1 ]=> (e 3)
;Value: 2.75

1 ]=> (e 4)
;Value: 2.7142857142857144

1 ]=> (e 5)
;Value: 2.71875

1 ]=> (e 10)
;Value: 2.7182817182817183

1 ]=> (e 100)
;Value: 2.718281828459045

练习 1.39

根据正切函数的连分式:

\[ \begin{equation} \begin{split} \tan x = \frac{x}{1- \frac{x^2}{3- \frac{x^3}{5- \ddots }}} \end{split} \end{equation} \]

定义过程(tan-cf x k)

这题直接将相应的\(N\)函数和\(D\)函数给出就行了。

其中:

  • \(N_i\)只有当\(i\)为 1 时返回\(x\),
  • 其他情况下都返回\(−(x2)\)(因为连分式内\(D_i\)进行的是减法操作)。

\(D_i\)则返回\(i\)所指定的位置的奇数,也即是:

  • 当\(i\)为\(1\)时,返回第一个奇数\(1\),
  • \(i\)为\(2\)时,返回第二个奇数\(3\)
  • 以此类推。
;;; 39-tan-cf.scm
(load "37-iter-cont-frac.scm")

(define (tan-cf x k)
    
	(define (N i)
		(if (= i 1)
			x
			(- (square x))))
	
	(define (D i)
		(- (* i 2) 1))
	
	(exact->inexact (cont-frac N D k)))

可以使用 MIT Scheme 内置的tan函数和tan-cf进行对比测试,验证tan-cf的正确性 :

1 ]=> (load "39-tan-cf.scm")
;Loading "39-tan-cf.scm"...
;  Loading "37-iter-cont-frac.scm"... done
;... done
;Value: tan-cf

1 ]=> (tan 10)
;Value: .6483608274590866

1 ]=> (tan-cf 10 100)
;Value: .6483608274590866

1 ]=> (tan 25)
;Value: -.13352640702153587

1 ]=> (tan-cf 25 100)
;Value: -.13352640702153587

过程作为返回值

结合不动点搜寻、平均阻尼、和\(y \mapsto x/y\)的例子:

之前求平均阻尼时所用的求平均数过程,所求的是\(y\)与\(x/y\)的平均值:

(lambda (y) (average y (/ x y)))

以后再写一个求\(y\)与\(y^2\)的平均数:

(lambda (y) (average y (square y)))

这里「求平均」是共同的,可以抽象成方法average-damp

(define (average-damp f)
	(lambda (y) (average y (f y))))
(average-damp (lambda (y) (/ x y))     ;; y与 x/y的平均值
(average-damp (lambda (y) (square y))  ;; y与 y^2的平均值

所以求平方根的过程与求立方根过程可以重写为:

(define (sqrt x)
	(fixed-point
		(average-damp (lambda (y) (/ x y)))
		1.0))

求立方根的过程可以重写为:

(define (cube-root x)
	(fixed-point 
		(average-damp (lambda (y) (/ x (square y))))
    1.0))

牛顿方法

如果\(x \mapsto g(x)\)是一个可微函数,则方程\(g(x)=0\)的一个解就是\(x \mapsto g(x)\) 的一个不动点,其中:

\[ \begin{equation} \begin{split} \label{eqnewton01} f(x) = x - \frac{g(x)}{Dg(x)} \end{split} \end{equation} \]

\(Dg(x)\)是\(g\)对\(x\)的导数,所以牛顿方法就是使用用找不动点的方法逼近方程的解。

函数到函数的变换

导数是从函数到函数的变换。例如:\(x \mapsto x^3\)的导数是\(x \mapsto 3x^2\)。

一般来说,对于函数\(g\)和一个很小的数\(dx\),那么函数\(g\)的导数由以下函数 (\(dx\)的极限)给出:

\[ \begin{equation} \begin{split} Dg(x) = \frac{g(x+dx)-g(x)}{dx} \end{split} \end{equation} \]

用过程描述导数的概念,它以函数为参数,返回值类型也是函数:

(define (deriv g)
	(lambda (x)
		(/ (- (g (+ x dx)) (g x)) dx)))

因为\(dx\)代表一个很小的数,以这里取值为\(0.00001\):

(define dx 0.00001)

例子,求\(x \mapsto x^3\)在\(5\)的导数的近似值:

(define (cube x) (* x x x))

((deriv cube) 5)      ;; result is 75.00014999664018

借助deriv把牛顿法表示为一个求不动点的过程,公式\ref{eqnewton01}可以描述为:

(define (newton-transform g)
	(lambda (x)
		(- x (/ (g x) ((deriv g) x)))))

然后再定义牛顿法:

(define (newtons-method g guess)
	(fixed-point (newton-transform g) guess))

它的参数是一个过程,然后希望找到零点的函数。这里还需要一个初始的猜测。

以求\(x\)的平方根为例,可以初始猜测值为\(1\),然后用牛顿法找\(y \mapsto y^2 -x \)的 零点。这样得到了求平方根函数的另一种形式:

(define (sqrt x) 
	(newtons-method (lambda (y) (- (square y) x)) 1.0))

抽象与第一级过程

之前介绍的平方根计算方式、求不动点、牛顿法都是类似的形式:

多一个函数出发,找到这个函数在某个变换下的不动点,可以用过程表示为:

(define (fixed-point-of-transform g transform guess)
  (fixed-point (transform g) guess))

三个参数分别为:过程gg的变换同方式、初始猜测值。

这样可以重新定义平方根计算(寻找\(y \mapsto x/y\)在平均阻尼下的不动点):

(define (sqrt x)
	(fixed-point-of-transform
		(lambda (y) (/ x y))
		average-damp
		1.0))

牛顿平方根(寻找\(y \mapsto y^2 -x\)牛顿变换的实例)也可以重新定义为:

(define (sqrt x)
	(fixed-point-of-transform
		(lambda (y) (- (square y) x))
		newton-transform
		1.0))

一般而言,程序语言会对计算元素可以使用的方式作出限制,而限制最少的状态被称为 具有 第一级 的状态。它的权限包括:

  • 可以使用变量名。
  • 可以把过程作为参数。
  • 可以把过程作为结果。
  • 可以包含在数据结构中。

实现一线过程的主要代价是,为了能把过程作为返回值,就需要为过程里的自由变量保留 空间,即使这一过程并不执行。在4.1节有关Scheme实现的研究中,这些变量都被存储在 过程的环境变量中。

练习 1.40

定义一个过程cubic,它可以作为newtons-method参数,像是这样:

(newtons-method (cubic a b c) 1)

可以用来逼近\(x^3 + ax^2 + bx + c\)的零点。

首先根据公式写出相应的cubic过程,它的返回值是一个接受参数x的过程:

;;; 40-cubic.scm
(load "8-cube.scm")

(define (cubic a b c) 
	(lambda (x) (
		+
		(cube x) 
		(* a (square x)) 
		(* b x) c)))

然后将cubic过程传给newtons-method就可以进行测试了:

1 ]=> (load "p50-newtons-method.scm")
;Loading "p50-newtons-method.scm"...
;  Loading "p46-fixed-point.scm"... done
;  Loading "p49-newton-transform.scm"...
;    Loading "p49-deriv.scm"... done
;  ... done
;... done
;Value: newtons-method

1 ]=> (load "40-cubic.scm")
;Loading "40-cubic.scm"...
;  Loading "8-cube.scm"... done
;... done
;Value: cubic

1 ]=> (newtons-method (cubic 3 2 1) 1.0)
;Value: -2.3247179572447267

1 ]=> (newtons-method (cubic 2 5 5) 1.0)
;Value: -1.2332293376819243

练习 1.41

定义一个过程double,以一个过程为参数,返回是另一个把参数过程执行两次的过程。

对于一个接受单个参数x的函数f来说,要将它应用多一次(总共两次)的办法是执行 以下表达式:

(f (f x))

由此可以给出相应的double函数,它接受一个函数f,并且返回一个能将f应用两次 的过程:

;;; 41-double.scm

(define (double f)
    (lambda (x)
        (f (f x))))

接着,用这个double函数来执行书本给出的表达式:

1 ]=> (load "41-double.scm")
;Loading "41-double.scm"... done
;Value: double

1 ]=> (((double (double double)) 1+) 5)
;Value: 21

练习 1.42

定义过程compose它的作用是复合两个函数数,如x \mapsto f(x)x \mapsto g(x) 组合就是\(x \mapsto f(g(x))\)。

例子,inc过程把参数的值加1,square把参数平方,那么:

((compose square inc) 6)   ;; 结果为(6 +1)^2 = 49

对于表达式:

((compose square inc) 6)

要计算出值 49 ,这个表达式应该执行以下计算序列:

((compose square inc) 6)

((lambda (x) (square (inc x))) 6)

(square (inc 6))

(square 7) 

49

也即是, (compose f g) 的展开式应该为:

(lambda (x) (f (g x)))

为这一模式赋予一个名字,得出的就是复合函数的定义了:

;;; 42-compose.scm

(define (compose f g) 
	(lambda (x) (f (g x))))

测试:

1 ]=> (load "42-compose.scm")
;Loading "42-compose.scm"... done
;Value: compose

1 ]=> ((compose square 1+) 6)
;Value: 49

练习 1.43

定义过程,把一个方法f重复调用n次,比如\(f(f(\cdots f(x)\cdots))\)。

对于表达式:

(repeated f n)

repeated 函数返回一个接受单个参数 x 的函数 g ,在这个函数 g 的函数体内, f 被调用了 n 次。

比如说,表达式(repeated square 4)的展开式应该为:

(repeated square 4)

(lambda (x) 
	(square ((repeated square 3) x)))

(lambda (x) 
	(square 
		((lambda (x) (square ((repeated square 2) x)))
			x)))

(lambda (x) 
	(square 
		((lambda (x) (square 
				((lambda (x) (square ((repeated square 1) x))) 
					x))) 
			x)))

(lambda (x) 
	(square 
		((lambda (x) (square 
				((lambda (x) (square (square x)))
					x))) 
			x)))
repeated

根据前面列出的规则,可以给出相应的repeated函数(递归计算):

;;; 43-repeated.scm
(define (repeated f n)
	(if (= n 1)
		f
		(lambda (x) 
			(let
				((fs (repeated f (- n 1))))
				(f (fs x))))))

无 let 版本

(define (repeated f n)
	(if (= n 1)
		f
		(lambda (x)
			(f ((repeated f (- n 1)) x)))))

repeated 也可以迭代地计算:

;;; 43-iter-repeated.scm
(define (repeated f n)
	(define (iter i repeated-f)
		(if (= i 1)
			repeated-f
			(iter (- i 1)
				(lambda (x) (f (repeated-f x))))))
	(iter n f))
带 compose 的 repeated

事实上,我们完全不必使用 lambda 来显式地组合起两个函数(因为这样容易出错), 使用 练习 1.42 的 compose 就可以更好地完成这件事。

以下是使用 compose 定义的递归计算的 repeated :

;;; 43-rec-repeated-using-compose.scm
(load "42-compose.scm")

(define (repeated f n)
	(if (= n 1)
		f
		(compose f (repeated f (- n 1)))))

迭代计算的 repeated 也可以用 compose 来定义:

;;; 43-iter-repeated-using-compose.scm
(load "42-compose.scm")

(define (repeated f n)
	(define (iter i repeated-f)
		(if (= i 1)
			repeated-f
			(iter (- i 1)
				(compose f repeated-f))))
	(iter n f))
测试

无 compose ,递归计算的 repeated :

1 ]=> (load "43-repeated.scm")
;Loading "43-repeated.scm"... done
;Value: repeated

1 ]=> ((repeated square 2) 5)
;Value: 625

无 compose ,迭代计算的 repeated :

1 ]=> (load "43-iter-repeated.scm")
;Loading "43-iter-repeated.scm"... done
;Value: repeated

1 ]=> ((repeated square 2) 5)
;Value: 625

带 compose ,递归计算的 repeated :

1 ]=> (load "43-rec-repeated-using-compose.scm")
;Loading "43-rec-repeated-using-compose.scm"...
;  Loading "42-compose.scm"... done
;... done
;Value: repeated

1 ]=> ((repeated square 2) 5)
;Value: 625

带 compose ,迭代计算的 repeated :

1 ]=> (load "43-iter-repeated-using-compose.scm")
;Loading "43-iter-repeated-using-compose.scm"...
;  Loading "42-compose.scm"... done
;... done
;Value: repeated

1 ]=> ((repeated square 2) 5)
;Value: 625
评论

参考快速幂乘的实现:

(define (repeated f n)
	(define (itr f n result)
		(cond ((= n 0) result)
			((even? n)
				(itr (compose f f) (/ n 2) result))
			(else
				(itr f (- n 1) (compose f result)))))
	(itr f n identity))
评论

我的racket代码

(define (repeated f n)
	(lambda (x)
		(if (<= n 0)
			x
			((compose f (repeated f (- n 1))) x))))
评论

这样就可以了…

(define (repeated f n)
	(if (= n 0)
		(lambda (x) x)
		(lambda (x)
			(f ((repeated f (- n 1)) x)))))
评论

想n等于0的时候(lambda (x) (x))的,发现不行…用了(lambda (x) (+ 0 x)), 求不挫的玩法

评论

好像没考虑n等于0的情况…就是为了知道lz怎么处理才过来看的…

评论

lz好像没考虑n等于0的情况。

练习 1.44

smooth过程实现平滑的功能,平滑是信号处理中的一个概念:

函数fdx是一个很小的值,那么f的平滑也是一个函数,它在点x的值就是 \(f(x - dx)\)、\(f(x)\)和\(f(x + dx)\)的平均值。

有时候需要进行多次平滑,请结合smoothrepeated实现多次平滑。

根据题目给出的定义,写出平滑函数 smooth :

;;; 44-smooth.scm
(define dx 0.00001)

(define (smooth f)
	(lambda (x)
		(/ (+ (f (- x dx))
				(f x)
				(f (+ x dx)))
			3)))

测试:

1 ]=> ((smooth square) 5)
;Value: 25.000000000066663
smooth-n-times

题目还要求我们写出一个可以连续进行 n 次平滑的函数(叫它 smooth-n-times 好了)。

根据描述,表达式(smooth-n-times f 3)的展开式应该为:

(smooth-n-times f 3)

(smooth (smooth-n-times f 2))

(smooth (smooth (smooth-n-times f 1)))

(smooth (smooth (smooth f)))

根据这个展开式,最简单直观的smooth-n-times可以这样定义:

;;; 44-smooth-n-times.scm
(load "44-smooth.scm")

(define (smooth-n-times f n)
	(if (= n 0)
		f
		(smooth (smooth-n-times f (- n 1)))))

上面的smooth-n-times是一个递归计算过程,smooth-n-times的迭代计算过程定义如下:

;;; 44-iter-smooth-n-times.scm
(load "44-smooth.scm")

(define (smooth-n-times f n)
	(define (iter i smoothed-f)
		(if (= i 0)
			smoothed-f
			(iter (- i 1)
				(smooth smoothed-f))))
	(iter n f))
使用 repeated 定义 smooth-n-times

事实上,在前面的smooth-n-times定义里,我们做了重复的工作,因为在 练习 1.43 里 ,我们就写过生成连续调用序列的函数repeated,而smooth-n-times里对smooth的 连续调用完全可以使用repeated来完成,从而避免显式的递归,以及off-by-one错误:

(repeated smooth 3)

(lambda (f)
	(smooth (repeated smooth 2)))

(lambda (f)
	(smooth (smooth (repeated smooth 1))))

(lambda (f)
	(smooth (smooth (smooth f))))

以下是使用 repeated 实现的 smooth-n-times 的定义:

;;; 44-smooth-n-times-using-repeated.scm
(load "44-smooth.scm")
(load "43-repeated.scm")

(define (smooth-n-times f n)
	(let ((n-times-smooth (repeated smooth n)))
		(n-times-smooth f)))
测试

无 repeated ,递归计算的 smooth-n-times :

1 ]=> (load "44-smooth-n-times.scm")
;Loading "44-smooth-n-times.scm"...
;  Loading "44-smooth.scm"... done
;... done
;Value: smooth-n-times

1 ]=> ((smooth-n-times square 10) 5)
;Value: 25.000000000666663

无 repeated ,迭代计算的 smooth-n-times :

1 ]=> (load "44-iter-smooth-n-times.scm")
;Loading "44-iter-smooth-n-times.scm"...
;  Loading "44-smooth.scm"... done
;... done
;Value: smooth-n-times

1 ]=> ((smooth-n-times square 10) 5)
;Value: 25.000000000666663

带 repeated 的 smooth-n-times (是递归计算还是迭代计算取决于所使用的 repeated 函数):

1 ]=> (load "44-smooth-n-times-using-repeated.scm")
;Loading "44-smooth-n-times-using-repeated.scm"...
;  Loading "44-smooth.scm"... done
;  Loading "43-repeated.scm"... done
;... done
;Value: smooth-n-times

1 ]=> ((smooth-n-times square 10) 5)
;Value: 25.000000000666663
评论

直接用

(((repeated smooth 10) square) 5)

调用似乎就可以了吧?

练习 1.45

之前已经研究过了:

  • 如何找\(y \mapsto x/y\)的不动点,结合平均阻尼求平方根。
  • 根据\(y \mapsto x/y^2\)的不动点,结合平均阻尼求立方根。

但是到了四次方就行不通了。一次平均阻尼在\(y \mapsto x/y^3\)不会收敛,要两次平均 阻尼(平均阻尼的平均阻尼)。

解决需要求\(n\)次方根时基于\(y \mapsto x/y^{n-1}\)需要反复做多少次平均阻尼。可以 使用fixed-pointaverage-damprepeated过程计算n次方根。

写出乘幂函数 expt

乘幂函数expt用于计算公式的y^{n−1}部分,它接受两个参数:basen,并计算出 base^n

乘幂函数的定义早在书本的29页就已经给出了,但是这里不妨用一种新的方式来实现它。

根据定义,乘幂其实就是对 base 进行 n 次自乘,比如说,\(2^5\)可以展开成计算序列:

\(2 \times 2 \times 2 \times 2 \times 2\)

而这个自乘序列又可以用以下的repeated调用表示(假设自由变量base已经被闭包为 数字2):

((repeated (lambda (x) (* base x)) 5) 1)

(* 2 ((repeated (lambda (x) (* base x)) 4) 1))

(* 2 (* 2 ((repeated (lambda (x) (* base x)) 3) 1)))

(* 2 (* 2 (* 2 ((repeated (lambda (x) (* base x)) 2) 1))))

(* 2 (* 2 (* 2 (* 2 ((repeated (lambda (x) (* base x)) 2) 1)))))

(* 2 (* 2 (* 2 (* 2 (* 2 1)))))

(* 2 (* 2 (* 2 (* 2 2))))

(* 2 (* 2 (* 2 4)))

(* 2 (* 2 8))

(* 2 16)

32

展开式的乘法序列还做了一点小变换,它将原本的

\(2 \times 2 \times 2 \times 2 \times 2\)

改成了

\(2 \times 2 \times 2 \times 2 \times 2 \times 1\)

这种变换并不影响计算结果,只是让repeated处理的条件简单一些,避免off-by-one 错误而已。

将这一展开模式写成相应的过程,我们就得出了使用repeated定义的expt函数:

;;; 45-expt.scm
(load "43-repeated.scm")

(define (expt base n)
	(if (= n 0)
		1
		((repeated (lambda (x) (* base x)) n) 1)))

测试:

1 ]=> (load "45-expt.scm")
;Loading "45-expt.scm"...
;  Loading "43-repeated.scm"... done
;... done
;Value: expt

1 ]=> (expt 2 0)
;Value: 1

1 ]=> (expt 2 1)
;Value: 2

1 ]=> (expt 2 5)
;Value: 32

1 ]=> (expt 2 10)
;Value: 1024
average-damp-n-times

因为需要对输入的公式进行不定数量的average-damp以确保不动点收敛,为了保持代码 的可读性,我们可以写一个辅助函数来做这件事:

;;; 45-average-damp-n-times.scm
(load "43-repeated.scm")
(load "p48-average-damp.scm")

(define (average-damp-n-times f n)
	((repeated average-damp n) f))

average-damp-n-times接受两个参数fn,并返回一个经过了naverage-damp 变换的f作为返回值。

测试:

1 ]=> (load "45-average-damp-n-times.scm")
;Loading "45-average-damp-n-times.scm"...
;  Loading "43-repeated.scm"... done
;  Loading "p48-average-damp.scm"...
;    Loading "p15-average.scm"... done
;  ... done
;... done
;Value: average-damp-n-times

1 ]=> ((average-damp-n-times square 10) 10.0)
;Value: 10.087890625

1 ]=> ((average-damp-n-times square 100) 10.0)
;Value: 10.
damped-nth-root

既然已经有了求幂函数expt,以及可以进行任意次平均阻尼变换的 average-damp-n-times,那么组合起这两个函数,再加上fixed-point,就可以写出求 n次方根的函数damped-nth-root了:

;;; 45-damped-nth-root.scm
(load "p46-fixed-point.scm")
(load "45-expt.scm")
(load "45-average-damp-n-times.scm")

(define (damped-nth-root n damp-times)
	(lambda (x)
		(fixed-point 
			(average-damp-n-times 
				(lambda (y) 
					(/ x (expt y (- n 1)))) 
				damp-times)
			1.0)))

函数damped-nth-root接受两个参数ndamp-timesn表示要计算的方根次数, damp-times指定要对公式进行多少次平均阻尼变换。

damped-nth-root的返回值是一个过程,它接受参数x,并计算xn次方根。

可以通过定义平方根、立方根和四次方根来测试damped-nth-root(因为暂时只知道 这三个方根需要多少次平均阻尼):

1 ]=> (load "45-damped-nth-root.scm")

;Loading "45-damped-nth-root.scm"...
;  Loading "p46-fixed-point.scm"... done
;  Loading "45-expt.scm"...
;    Loading "43-repeated.scm"... done
;  ... done
;  Loading "45-average-damp-n-times.scm"...
;    Loading "43-repeated.scm"... done
;    Loading "p48-average-damp.scm"...
;      Loading "p15-average.scm"... done
;    ... done
;  ... done
;... done
;Value: damped-nth-root

1 ]=> (define sqrt (damped-nth-root 2 1))
;Value: sqrt

1 ]=> (sqrt (* 3 3))
;Value: 3.

1 ]=> (define cube-root (damped-nth-root 3 1))
;Value: cube-root

1 ]=> (cube-root (* 3 3 3))
;Value: 2.9999972321057697

1 ]=> (define 4th-root (damped-nth-root 4 2))
;Value: 4th-root

1 ]=> (4th-root (* 3 3 3 3))
;Value: 3.000000000000033
收敛条件

接着要解决的问题是,找出计算 n 次方根和收敛计算所需的平均阻尼次数之间的关系, 以下是一些实验数据:

n次方根 1 2 3 4 5 6 7 8 ... 15 16 ... 31 32 ...
收敛所需的平均阻尼次数 1 1 1 2 2 2 2 3 ... 3 4 ... 4 5 ....

可以看出,要使得计算\(n\)次方根的不动点收敛,最少需要\(\lfloor \lg n \rfloor\)次平均阻尼。

其中\(lg\)可以用函数lg计算得出:

;;; 45-lg.scm
(define (lg n)
	(cond 
		((> (/ n 2) 1)
			(+ 1 (lg (/ n 2))))
		((< (/ n 2) 1) 0)
		(else 1)))

另外,lg函数已经完成了取整的工作,因此计算完(lg n)之后,就不必再使用其他函数 进行向下取整工作了。

测试:

1 ]=> (load "45-lg.scm")
;Loading "45-lg.scm"... done
;Value: lg

1 ]=> (lg 0)
;Value: 0

1 ]=> (lg 1)
;Value: 0

1 ]=> (lg 2)
;Value: 1

1 ]=> (lg 3)
;Value: 1

1 ]=> (lg 1024)
;Value: 10
nth-root

现在可以给出函数nth-root的定义了,它调用damped-nth-root函数计算n次方根, 并使用lg函数计算足以令不动点收敛的平均阻尼次数:

;;; 45-nth-root.scm
(load "45-damped-nth-root.scm")
(load "45-lg.scm")

(define (nth-root n)
	(damped-nth-root n (lg n)))

最终得出的nth-root函数非常简单,因为所有工作都已经在子函数里完成了。

可以通过定义一些次方根来测试nth-root

1 ]=> (load "45-nth-root.scm")
;Loading "45-nth-root.scm"...
;  Loading "45-damped-nth-root.scm"...
;    Loading "p46-fixed-point.scm"... done
;    Loading "45-expt.scm"...
;      Loading "43-repeated.scm"... done
;    ... done
;    Loading "45-average-damp-n-times.scm"...
;      Loading "43-repeated.scm"... done
;      Loading "p48-average-damp.scm"...
;        Loading "p15-average.scm"... done
;      ... done
;    ... done
;  ... done
;  Loading "45-lg.scm"... done
;... done
;Value: nth-root

1 ]=> (define sqrt (nth-root 2))
;Value: sqrt

1 ]=> (sqrt (* 3 3))
;Value: 3.

1 ]=> (define cube-root (nth-root 3))
;Value: cube-root

1 ]=> (cube-root (* 3 3 3))
;Value: 2.9999972321057697

1 ]=> (define 4th-root (nth-root 4))
;Value: 4th-root

1 ]=> (4th-root (* 3 3 3 3))
;Value: 3.000000000000033

1 ]=> (define 100th-root (nth-root 100))
;Value: 100th-root

1 ]=> (100th-root 100)
;Value: 1.047130656622199

练习 1.46

本章很多算法都是「迭代式改进」形式的:先猜一个值,检查是不是够好,不够好就改进。

把这种模式抽象成一个过程interative-improve,它有两个参数类型都是过程:

  1. 检查结果是否足够好
  2. 改进猜测的方法

interative-improve返回的结果也是一个过程,以某一个猜测为参数,通过不断改进, 得到结果。

interative-improve重写1.1.7节的sqrt过程和1.3.3节的fixed-point过程。

题目要求我们给出这样一个iterative-improve函数:它接受一个用于检测猜测值是否 足够好的函数(close-enough?),以及一个用于改进猜测值的函数(improve),并返回 一个接受单个初始猜测值(first-guess)的过程,这个过程可以一直改进猜测值,直到 猜测足够好。

根据描述,可以给出以下形式的函数:

(define (iterative-improve close-enough? improve)
    (lambda (first-guess)
        ; ...
    ))

这个过程和 46 页的fixed-point非常的相似:

;;; p46-fixed-point.scm

(define tolerance 0.00001)

(define (fixed-point f first-guess)
	(define (close-enough? v1 v2)
		(< (abs (- v1 v2)) tolerance))
	(define (try guess)
		(let ((next (f guess)))
			(if (close-enough? guess next)
				next
				(try next))))
	(try first-guess))

fixed-point 函数不仅仅和iterative-improve非常相似,事实上,iterative-improve 就隐藏在fixed-point当中!

fixed-point中,close-enough?负责检查猜测是否足够好,而函数f负责改进猜测 ,如果我们将close-enough?函数抽取出来,成为额外的参数,那么fixed-point 的定义就变成了:

(define (fixed-point f first-guess close-enough?)
	(define (try guess)
		(let ((next (f guess)))
			(if (close-enough? guess next)
				next
				(try next))))
	(try first-guess))

(define (close-enough? v1 v2)
	(< (abs (- v1 v2)) tolerance))

(define tolerance 0.00001)

如果再将first-guessfixed-point函数的参数列表中移走,变成另一个包裹 try函数的lambda的参数,fixed-point函数的定义就变成了这样:

(define (fixed-point f close-enough?)
	(lambda (first-guess)
		(define (try guess)
			(let ((next (f guess)))
				(if (close-enough? guess next)
					next
					(try next))))
		(try first-guess)))

(define (close-enough? v1 v2)
	(< (abs (- v1 v2)) tolerance))

(define tolerance 0.00001)

可以看到,现在的first-point定义已经和前面给出的iterative-improve形式一样了, 如果将fixed-point函数改名成iterative-improve,将参数f改名成improve, 并且删除close-enough函数和dx变量,题目要求的iterative-improve就(神奇地) 显出庐山真面目了:

;;; 46-iterative-improve.scm
(define (iterative-improve close-enough? improve)
	(lambda (first-guess)
		(define (try guess)
			(let ((next (improve guess)))
				(if (close-enough? guess next)
					next
					(try next))))
		(try first-guess)))

注意我们是如何一步步地从fixed-point函数中抽象出iterative-improve函数的,将 这两个函数放在一起对比将是一个非常有趣的练习。

iterative-improve重定义fixed-point
;;; 46-fixed-point.scm
(load "46-iterative-improve.scm")

(define (fixed-point f first-guess)
	(define tolerance 0.00001)
	(define (close-enough? v1 v2)
		(< (abs (- v1 v2)) tolerance))
	(define (improve guess)
		(f guess))
	((iterative-improve close-enough? improve) first-guess))

测试:

1 ]=> (load "46-fixed-point.scm")
;Loading "46-fixed-point.scm"...
;  Loading "46-iterative-improve.scm"... done
;... done
;Value: fixed-point

1 ]=> (fixed-point cos 1.0)
;Value: .7390822985224024
iterative-improve重定义 sqrt
;;; 46-sqrt.scm
(load "46-iterative-improve.scm")

(define (sqrt x)
	(define dx 0.00001)
	(define (close-enough? v1 v2)
		(< (abs (- v1 v2)) dx))
	(define (improve guess)
		(average guess (/ x guess)))
	(define (average x y)
		(/ (+ x y) 2))
	((iterative-improve close-enough? improve) 1.0))

测试:

1 ]=> (load "46-sqrt.scm")
;Loading "46-sqrt.scm"...
;  Loading "46-iterative-improve.scm"... done
;... done
;Value: sqrt

1 ]=> (sqrt 9)
;Value: 3.
评论

这样是不是表述更清晰些……

(define (iterative-improve close-enough? improve)
	(define (try guess)
		(let ((next (improve guess)))
			(if (close-enough? guess next)
				next
				(try next))))
	(lambda (first-guess) (try first-guess)))
评论

是的。楼主在iterative-improve里用lambda多此一举。

(define (iterative-improve good-enough? improve)
	(define (iter x)
		(display x)
		(newline)
		(if (good-enough? x)
			x
			(iter (improve x))))
	iter)

(define (fixed-point f first-guess)
	(define (good-enough? x)
		(< (abs (- x (f x))) tolerance))
	((iterative-improve good-enough? (average-damp f)) first-guess))