1.3 高階手続きによる抽象

高階手続き（手続きを返したり受け取ったりする手続き）による強力な抽象の力を見ていく。

1.3.1 引数としての手続き

問題1.29

(define (inc n) (+ n 1))
(define (cube n) (* n n n))
(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 h (/ (- b a) n))
  (define (y k)
    (f (+ a (* k h))))
  (define (term k)
    (cond ((= k 0) (y 0))
          ((= k n) (y n))
          ((even? k) (* 2 (y k)))
          (else (* 4 (y k)))))
  (* (/ h 3)
     (sum term 0 inc n)))

gosh> (simpson cube 0 1 100)
1/4
gosh> (simpson cube 0 1 1000)
1/4

有理数になってる.ばっちり正確なわけですが.

gosh> (simpson cube 0.0 1 100)
0.24999999999999992
gosh> (simpson cube 0.0 1 1000)
0.2500000000000003

integralとの比較

(define (integral f a b dx)
  (define (add-dx x) (+ x dx))
  (* (sum f (+ a (/ dx 2.0)) add-dx b)
     dx))

gosh> (integral cube 0 1 0.01)
0.24998750000000042
gosh> (integral cube 0 1 0.001)
0.249999875000001

simpsonの方が正確ですな。
実行時間を計測してみよう.

(time (dotimes (i 5000) (integral cube 0 1 0.001)))
;(time (dotimes (i 5000) (integral cube 0 1 0.001)))
; real   4.935
; user   4.880
; sys    0.030
(time (dotimes (i 5000) (simpson cube 0.0 1 1000)))
;(time (dotimes (i 5000) (simpson cube 0.0 1 1000)))
; real   7.369
; user   7.340
; sys    0.020

ふむ。simpsonの方が遅い.ちょっと、実装を工夫してみよう.condの条件分岐をなくしてみる.

(define (simpson2 f a b n)
  (define h (/ (- b a) n))
  (define (y k)
    (f (+ a (* k h))))
  (define (inc2 n) (+ n 2))
  (define (term n)
    (lambda (k)
      (* n (y k))))
  (* (/ h 3)
     (+ (y 0)
        (sum (term 4) 1 inc2 (- n 1))
        (sum (term 2) 2 inc2 (- n 1))
        (y n))))

これで、integralとほとんど同じ実行時間になった。

問題1.30

(define (sum term a next b)
  (define (iter a result)
    (if (> a b)
        result
        (iter (next a) (+ result (term a)))))
  (iter a 0))
  
gosh> (time (sum identity 0 inc 1000000))
;(time (sum identity 0 inc 1000000))
; real   0.684
; user   0.690
; sys    0.000
500000500000

gosh> (time (sum2 identity 0 inc 1000000))
;(time (sum2 identity 0 inc 1000000))
; real   0.328
; user   0.330
; sys    0.000
500000500000

sum2は今回定義したもの。倍ぐらい速くなっている。

問題1.31

;; 反復的
(define (product term a next b)
  (define (iter a result)
    (if (> a b)
        result
        (iter (next a) (* result (term a)))))
  (iter a 1))
;; 再帰的
(define (product term a next b)
  (if (> a b)
      1
      (* (term a)
         (product term (next a) next b))))

(define (factorial n)
  (product identity 1 inc n))

(define (pi/4 n)
  (define (inc n) (+ n 1))
  (define (term n)
    (cond ((even? n) (/ (+ n 2) (+ n 1)))
          (else (/ (+ n 1) (+ n 2)))))
  (product term 1.0 inc n))

    
  
gosh> (* 4 (pi/4 10.0))
3.2751010413348065
gosh> (* 4 (pi/4 100))
3.1570301764551654
gosh> (* 4 (pi/4 1000))
3.1431607055322552
gosh> (* 4 (pi/4 10000))
3.141749705738009
gosh> (* 4 (pi/4 100000))
3.1416083612779406

円周率の収束は遅いな。

問題1.32

sum, productはそれぞれ、accumerateを使って以下のように定義できる.

(define (accumerate combiner nulll-avlue term a next b)
  (if (> a b)
      null-value
      (combiner (term a)
                (accumerate combiner null-value term (next a) next b))))
(define (sum term a next b)
  (accumerate + 0 term a next b))

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

;; 反復的なaccumerate
(define (accumerate 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))

でも、これだと、null-valueは項がなくなったときに使うという条件を満たしていない. *, + は交換則が成り立つので、(誤差を除けば）値は代わらないはずだが.

;; null-valueを最後に使う
(define (accumerate2 combiner null-value term a next b)
  (define (iter a result)
    (if (> a b)
        (combiner result null-value)
        (iter (next a) (combiner result (term a)))))
  (if (> a b)
      null-value
      (iter (next a) (term a))))

問題1.33

filtered-accumulateの定義で明らかに再帰呼び出しの部分が重複したコードになっているけど，これは後に出てくるletを使うことで重複をなくすことができる。内部定義でもできるけどあまりいいスタイルではないと思う。

(define (filtered-accumulate combiner null-value term a next b filter)
  (cond ((> a b) null-value)
        ((filter a)
         (combiner (term a)
                   (filtered-accumulate combiner
                                        null-value
                                        term
                                        (next a)
                                        next
                                        b
                                        filter)))
        (else 
         (filtered-accumulate combiner
                              null-value
                              term
                              (next a)
                              next
                              b
                              filter))))
  
a. 区間a, bの素数の2乗の和
(define (sum-of-square-prime-range a b)
  (define (inc n) (+ n 1))
  (filtered-accumulate + 0 square a inc b prime?))

b. nと互いに素で,nより小さい正の整数の積
(define (sum-of-relatively-prime n)
  (define (inc n) (+ n 1))
  (define (gcd a b)
    (if (= b 0) a (gcd b (remainder a b))))
  (define (filter a)
    (= 1 (gcd n a)))
  (filtered-accumulate + 0 identity 1 inc n filter))

1.3.2 lambdaを使う手続きの構築

ついに、lambdaが出てきた.
内部定義は、手続きのためだけにしたほうがいいってある。でも、インデントが深くなるから手続き以外でも使いたいよなぁ.

問題1.34

(f f)
(f 2)
(2 2) ==> エラー. 2は手続きではない！

1.3.3 一般方法としての手続き

fixed-pointって、fや与える値によっては、不動点があったとしても無限ループになるのでは？その後の説明できっと出てくるだろうけどやってみる。

(define tolerance 0.0001)
(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))
;; f(x) = -x + 2
(define (f x) (+ (- x) 2))
(fixed-point f 1.0) ==> 1.0
(fixed-point f 0.5) ==> 無限ループ(0.5と1.5を繰り返す)

問題1.35

Wikipediaによると、黄金比とは a : b = b : (a + b) が成り立つように分割したときの a : bのことである.とある。a = 1とすると、 b^2 = 1 + b <=> b = 1 + 1/bとなるので、黄金比φは x -> 1 + 1/x の不動点である.

平均緩和法を使わずに収束するのだろうか？たぶんしない。

(fixed-point (lambda (x) (+ 1 (/ 1 x))) -10.0)
=> 1.618060695780903

いや、収束した.

(fixed-point (lambda (x) (/ 1 x)) -10.0)

は収束しない（プロセスが止まらない）。y = x に対して対称かどうかの差だろう.

問題1.36

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

;; 平均緩和法を使わない
(fixed-point (lambda (x) (/ (log 1000) (log x))) 10.0)
gosh> 10.0
2.9999999999999996
6.2877098228681545
3.7570797902002955
5.218748919675316
4.1807977460633134
4.828902657081293
4.386936895811029
4.671722808746095
4.481109436117821
4.605567315585735
4.522955348093164
4.577201597629606
4.541325786357399
4.564940905198754
4.549347961475409
4.5596228442307565
4.552843114094703
4.55731263660315
4.554364381825887
4.556308401465587
4.555026226620339
4.55587174038325
4.555314115211184
4.555681847896976
4.555439330395129
4.555599264136406
4.555493789937456
4.555563347820309
-> 29ステップ

;; 平均緩和法を使う
(fixed-point (lambda (x) (/ (+ x (/ (log 1000) (log x))) 2)) 10.0)
gosh> 10.0
6.5
5.095215099176933
4.668760681281611
4.57585730576714
4.559030116711325
4.55613168520593
4.555637206157649
4.55555298754564
-> 9ステップ

平均緩和法を使わない方は、解の回りで振動している.それに対し平均緩和法を使った方は振動しておらず、収束が早い.

問題1.37

(連分数は英語でcontinued fractionっていうのか)

a.

(define (cont-frac n d k)
  (define (f i)
    (if (= i k)
        (/ (n i) (d i))
        (/ (n i) (+ (d i)
                    (f (+ i 1))))))
  (f 1))

φ = 0.6 ... だから、4桁までの精度ということは黄金比との差が5e-5に納めればってことかな。

(define (check a b)
  (define (close-enougth? n)
    (let ((phi (/ (+ 1 (sqrt 5)) 2)))
      (< (abs (- n (/ 1 phi)))
         5e-5)))
  (if (<= a b)
      (let ((n (cont-frac (lambda (i) 1.0)
                          (lambda (i) 1.0)
                          a)))
        (display a)
        (display " ")
        (display n)
        (if (close-enougth? n)
            (display " ok")
            #f)
        (newline)
        (check (+ a 1) b))
      #f))

-- 実行結果
gosh> (check 1 30)
1 1.0
2 0.5
3 0.6666666666666666
4 0.6000000000000001
5 0.625
6 0.6153846153846154
7 0.6190476190476191
8 0.6176470588235294
9 0.6181818181818182
10 0.6179775280898876
11 0.6180555555555556 ok
12 0.6180257510729613 ok
13 0.6180371352785146 ok
14 0.6180327868852459 ok
15 0.6180344478216819 ok
16 0.6180338134001252 ok
17 0.6180340557275542 ok
18 0.6180339631667064 ok
19 0.6180339985218034 ok
20 0.6180339850173578 ok
21 0.6180339901755971 ok
22 0.6180339882053251 ok
23 0.618033988957902 ok
24 0.6180339886704432 ok
25 0.6180339887802426 ok
26 0.6180339887383031 ok
27 0.6180339887543225 ok
28 0.6180339887482037 ok
29 0.6180339887505407 ok
30 0.6180339887496482 ok
#f

というわけで、k>=11で4桁の精度がえられる.

b.
先ほどのは再帰的なので反復的なcont-frac

(define (cont-frac n d k)
  (define (iter k a)
    (if (= k 0)
        a
        (iter (- k 1) (/ (n k) (+ (d k) a)))))
  (iter k 0))

問題1.38

うまくDの一般項を見つける必要がある。`

(define (euler-e k)
  (+ 2
     (cont-frac (lambda (i) 1.0)
             (lambda (i)
               (let ((n (remainder i 3)))
                 (cond ((= n 2) (/ (+ i i 2.0) 3.0))
                       (else 1.0))))
             k)))

gosh> (euler-e 10)
2.7182817182817183
gosh> (euler-e 100)
2.7182818284590455
gosh> (euler-e 1000)
2.7182818284590455
gosh> (euler-e 10000)
2.7182818284590455

実際eの値は
2.718281828459045
なので k = 100で良い近似値がえられている.

問題1.39

(define (tan-cf x k)
  (cont-frac (lambda (k)
                 (if (= k 1)
                       x
                       (- (* x x))))
               (lambda (k)
                 (+ k k -1.0))
               k))

テスト

(define (test x a b d)
  (define (loop a)
    (if (> a b)
        #f
        (begin (format #t "~4a ~20a ~a\n"
                       a
                       (tan-cf x a)
                       (- (tan x) (tan-cf x a)))
               (loop (+ a d)))))
  (print "tan(x) = " (tan x))
  (format #t "~4a ~20a ~a\n" "k" "tan-cf" "diff")
  (loop a))

tan(x) = 1.5574077246549023
k    tan-cf               diff
1    1.0                  0.5574077246549023
2    1.4999999999999998   0.057407724654902514
3    1.5555555555555558   0.00185216909934649
4    1.5573770491803278   3.067547457447084e-5
5    1.5574074074074076   3.1724749471884195e-7
6    1.557407722401769    2.253133235541327e-9
7    1.5574077246432194   1.1682876888130522e-11
8    1.557407724654856    4.618527782440651e-14
9    1.557407724654902    2.220446049250313e-16
10   1.557407724654902    2.220446049250313e-16
11   1.557407724654902    2.220446049250313e-16
12   1.557407724654902    2.220446049250313e-16
13   1.557407724654902    2.220446049250313e-16
14   1.557407724654902    2.220446049250313e-16
15   1.557407724654902    2.220446049250313e-16
16   1.557407724654902    2.220446049250313e-16
17   1.557407724654902    2.220446049250313e-16
18   1.557407724654902    2.220446049250313e-16
19   1.557407724654902    2.220446049250313e-16
20   1.557407724654902    2.220446049250313e-16
#f

x=1.0の時は、k = 10からは値の変化が無い. xの絶対値が大きくなる程に収束が遅くなる.x=100の時はk=133で収束する.xがどんな値でも、tan(x)があまり大きく（小さく）ならない範囲ではほぼ正しい値に値に収束する.

1.3.4 値として返される手続き

数列の収束を不動点として考えるっていうのは初めてだ.
X(n+1) = Xn / 2
という数列の収束する値は
f(x) = x / 2
の不動点ってことか。

Newton法は図に書いてみると一目瞭然だ.というか高校でこれやった気がする.

第一級身分という概念が出てくる.今までにも聞いたことはあるけど明確な定義は初めて。

問題1.40

(define dx 0.00001)
(define (deriv g)
  (lambda (x)
    (/ (- (g (+ x dx)) (g x))
       dx)))
(define (newton-transform g)
  (lambda (x)
    (- x (/ (g x) ((deriv g) x)))))
(define (newtons-method g guess)
  (fixed-point (newton-transform g) guess))

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

(newtons-method (cubic 1 2 3) 1) ==> -1.2756822036498454
((cubic 1 2 3) -1.2756822036498454) ==> 4.935607478273596e-12

問題1.41

doubleの定義は

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

(((double (double double)) inc) 5)
がどういう値を返すか、実行する前に考えてみる.

(double (double double))

まず、(double double)は手続きを取り、4回作用させる手続きを返す手続きだ.それを、doubleに適用すると手続きを取り、それに4回の作用をを4回作用させる手続きになる。4回の作用を4回行うということは4*4=4^2=16回作用させるということだ.なので、問題の式の値は21になる。

式の置き換えをしてみる。

(((double (double double)) inc) 5)
==> (((double
       (lambda (x) (double (double x))))
      inc)
     5)
==> (((lambda (x)
        ((lambda (x) (double (double x)))
         ((lambda (x) (double (double x))) x)))
      inc)
     5)
==> ((((lambda (x) (double (double x)))
         ((lambda (x) (double (double x))) inc))
      )
     5)
==> ((((lambda (x) (double (double x)))
         (double (double inc)))
      )
     5)
==> (((double (double (double (double inc)))))
     5)

ok. doubleが4個なので2^4=16。
実際に評価すると

gosh> (((double (double double)) inc) 5)
21

問題1.42

composeはgaucheに既にあるので、名前をcompose2とする。

(define (compose2 f g) (lambda (x) (f (g x))))

問題1.43

;; 再帰的
(define (repeated f n)
  (lambda (x)
    (define (iter n)
      (if (= n 0)
          x
          (f (iter (- n 1)))))
    (iter n)))

;; 再帰的2
(define (repeated2 f n)
  (lambda (x)
    (if (= n 0)
        x
        (f ((repeated2 f (- n 1)) x)))))

;; 反復的
(define (repeated3 f n)
  (lambda (x)
    (define (iter cnt result)
      (if (= cnt 0)
          result
          (iter (- cnt 1) (f result))))
    (iter n x)))

repeated2は内部で新たな関数を使っていないが、毎回手続きの生成とその評価をしているので遅い.
あ、ヒントを見逃してた.composeを使うと

(define (repeated4 f n)
  (if (= n 0)
      identity
      (compose f (repeated4 f (- n 1)))))

だ。repeated2と考え方は同じだ.

問題1.44

(define (smooth f)
  (define dx 0.1)
  (lambda (x)
    (/ (+ (f (- x dx)) (f x) (f (+ x dx))) 3)))
(define (fold-smoothed f n)
  ((repeated smooth n) f))

問題1.45

収束しない場合、おそらく同じ値を繰り返すことになるだろうという予想の元、確かめる手続きを作ってみる.

(define (average-damp f)
  (lambda (x) (/ (+ x (f x)) 2)))

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


(define (root x n k)
  (fixed-point ((repeated average-damp k)
              (lambda (y)
                (/ x (expt y (- n 1)))))
             1.0))

このrootを使って、nを増やしていったときに収束する最小のkを調べた. x には 2を使った.結果は以下の表.

n	k
1	0
2	1
3	1
4	2
5	2
6	2
7	2
8	3
9	3
15	3
16	4
31	4
32	5
63	5
64	6
65	6

どうやら必要な平均緩和の回数は
$\lfloor \log_2{x} \rfloor$
になる。で、n乗根を求める手続きはこうなる.

(define (root n x)
  (fixed-point ((repeated average-damp (inexact->exact (round (/ (log n)
                                                                 (log 2)))))
                (lambda (y) (/ x (expt y (- n 1)))))
               1.0))

問題1.46

(define (iterative-improve good-enough? improve)
  (define (iter guess)
    (if (good-enough? guess)
        guess
        (iter (improve guess))))
  iter)

(define (sqrt x)
  (define (good-enough? guess)
    (< (abs (- (* guess guess) x)) 0.001))
  (define (improve guess)
    (/ (+ guess (/ x guess)) 2))
  ((iterative-improve good-enough? improve) 1.0))

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