SICP 1.3.3 Procedures as General Methods - プログラミング再入門

ノート

汎用メソッドとしての手続き。

Finding roots of equations by the half-interval method

二分方による方程式の解（根）。

f(a) < 0 < f(b)の時に、aとbの間にあるf(x)=0となるxを求める。aとbの中間点xに対してf(x)が正ならbをxと置き換え、負ならaをxと置き換えて繰り返す。いずれf(x)が十分に小さくなったらxを答えと見なす。

紹介されている定義を実際に動かしてみる。追加で必要な関数はaverageのみ。

(define (average a b)
  (/ (+ a b) 2))

実行結果

> (half-interval-method sin 2.0 4.0)
3.14111328125
> (half-interval-method (lambda (x) (- (* x x x) (* 2 x) 3))
                        1.0
                        2.0)
1.89306640625
>

Finding fixed points of functions

関数の不動点を探す。
f(x)=xを満たすxが不動点。例を実行してみる。

> (fixed-point cos 1.0)
0.7390822985224023
> (fixed-point (lambda (y) (+ (sin y) (cos y)))
               1.0)
1.2587315962971173
>

ある値xの平方根を求めるには $y^2=x$ を満たす $y$ を探す事になる。仮に $Y$ がこの条件を満たしている場合 $Y=\frac{x}{Y}$ となる筈で、これがちょうど $y=f(y)$ の形になっているので $f(y)=\frac{x}{y}$ の不動点を見つけると、その $y$ は $Y$ である事になる。計算の過程では $y^2$ を $y_1\times y_2$ に分解して $y_1\times y_2=x$ から $y_2=\frac{x}{y_1}$ と言う式でyの不動点を求めると考える。

ところが、 $y$ の初期値を $y_1$ として順に $y_2$ 、 $y_3$ と求めて行こうとすると $y_2=\frac{x}{y_1}$ 、 $y_3=\frac{x}{y_2}=\frac{x}{\frac{x}{y_1}}=y_1$ となり $y_1$ と $y_2$ の間を行き来するだけで収束しない。

$y_2=\frac{x}{y_1}$ で $y_1$ と $y_2$ が離れていると仮定して、次の $y$ を $y_2$ ではなく $y_1$ と $y_2$ の中間点とする事で行ったり来たりを回避する。 $y$ の移動が半分になるので戻ろうとしても元の $y$ には戻って来ない。

組み込みのsqrtは再定義出来ないので、mysqrtとして定義して

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

実行してみる。

> (mysqrt 5)
2.236067977499978
> (sqrt 5)
2.23606797749979
>

このテクニックをavarage dumping（平均緩和法）と呼ぶらしい。

Exercise 1.35

1.2.2に示された黄金比の定義は
$\phi^2=\phi+1$
両辺を $\phi$ で割って
$\phi=\frac{\phi}{\phi} + \frac{1}{\phi} = 1 + \frac{1}{\phi}$
定義の通りfixed-pointで計算する。大体1.6と分かっているので1.0を初期値とする。同じく定義から(1 + √5) / 2と比較してみる。

> (fixed-point (lambda (φ) (+ 1 (/ 1 φ))) 1.0)
1.6180327868852458
> (/ (+ 1 (sqrt 5)) 2)
1.618033988749895
>

Exercise 1.36

tryにdisplayとnewlineを挿入。

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

一応式の変形をおさらいしておく。
$x^x=1000$
両辺の対数を取って
$\log{x^x}=\log{1000}$
$x\log{x}=\log{1000}$
$x=\frac{\log{1000}}{\log{x}}$
普通に計算する。1を初期値にすると0除算となるとの事なので、2から始める。

> (fixed-point (lambda (x) (/ (log 1000) (log x))) 2.0)
2.0
9.965784284662087
3.004472209841214
6.279195757507157
3.7598507024015393
5.2158437849258945
4.182207192401398
4.82776509834459
4.387593384662677
4.671250085763899
4.481403616895052
4.6053657460929
4.5230849678718865
4.577114682047341
4.541382480151454
4.564903245230833
4.549372679303342
4.559606491913287
4.552853875788271
4.557305529748263
4.554369064436181
4.556305311532999
4.555028263573554
4.555870396702851
4.555315001192079
4.5556812635433275
4.555439715736846
4.555599009998291
4.555493957531389
4.555563237292884
4.555517548417651
4.555547679306398
4.555527808516254
4.555540912917957
4.555532270803653
>

average dumpingを使う。

> (fixed-point (lambda (x) (average x (/ (log 1000) (log x)))) 2.0)
2.0
5.9828921423310435
4.922168721308343
4.628224318195455
4.568346513136243
4.5577305909237005
4.555909809045131
4.555599411610624
4.5555465521473675
4.555537551999825
>

圧倒的に速く収束する。

Exercise 1.37

a.
continued fractionとは連分数の事。連分数とは分母に更に分数が含まれている分数の事。で、無限連分数からk項の有限連分数。
最初の単純な回答

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

試してみる。

> (/ 1 (cont-frac (lambda (i) 1.0)
                  (lambda (i) 1.0)
                  10))
1.6181818181818184
> (/ 1 (cont-frac (lambda (i) 1.0)
                  (lambda (i) 1.0)
                  100))
1.618033988749895
>

それなりに合っていそうだが、定義に(n k) (d k)が２回ずつ登場するのがちょっと気に入らない。定義の式をそのまま関数にしてみる。

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

これはこれで、+ 0ってのが無駄ではある。。。ま、実行してみる。

> (/ 1 (cont-frac (lambda (i) 1.0)
                  (lambda (i) 1.0)
                  10))
1.6181818181818184
> (/ 1 (cont-frac (lambda (i) 1.0)
                  (lambda (i) 1.0)
                  100))
1.618033988749895
>

一応、合っていそう。
問題は４桁合わせるのにkはいくつ必要なのか。場合によっては誤差は非常に小さくても５桁目を切り捨てた時に４桁目が合っているとは限らないので、「4桁合っている」を「誤差が0.001未満」と解釈する。
プログラムで求めてみる。繰り返し回数を引数に取る関数f、正確な値ans、最大誤差toleranceを引数に取る関数iteration

(define (iteration f ans tolerance)
  (define (iter i)
    (let* ((candidate (f i))
           (error (abs (- candidate ans))))
      (display candidate)
      (display " (")
      (display error)
      (display ")")
      (newline)
      (if (<  error tolerance)
          i
          (iter (+ i 1)))))
  (iter 1))

途中経過を表示する様にした。実行結果

> (iteration (lambda (k) (/ 1 (cont-frac (lambda (i) 1.0)
                                         (lambda (i) 1.0)
                                         k)))
             (/ (+ 1 (sqrt 5)) 2)
             0.001)
1.0 (0.6180339887498949)
2.0 (0.3819660112501051)
1.5 (0.1180339887498949)
1.6666666666666665 (0.048632677916771616)
1.6 (0.018033988749894814)
1.625 (0.0069660112501050975)
1.6153846153846154 (0.0026493733652794837)
1.619047619047619 (0.0010136302977241662)
1.6176470588235294 (0.00038692992636546464)
9
> (iteration (lambda (k) (/ 1 (cont-frac (lambda (i) 1.0)
                                         (lambda (i) 1.0)
                                         k)))
             (/ (+ 1 (sqrt 5)) 2)
             0.0001)
1.0 (0.6180339887498949)
2.0 (0.3819660112501051)
1.5 (0.1180339887498949)
1.6666666666666665 (0.048632677916771616)
1.6 (0.018033988749894814)
1.625 (0.0069660112501050975)
1.6153846153846154 (0.0026493733652794837)
1.619047619047619 (0.0010136302977241662)
1.6176470588235294 (0.00038692992636546464)
1.6181818181818184 (0.0001478294319234852)
1.6179775280898876 (5.6460660007306984e-05)
11
>

11項必要。

b.
aで作成したcont-fracは再帰的（再帰呼び出しの戻り道で計算する）なので、繰り返し的（再帰呼び出しの行きで計算する）に書き換える。
が、繰り返し計算が必要な部分が分母の足し算なので（掛け算であれば簡単だが）、添字の1からkに向かって計算するのではなく、逆にkから1に向かって計算する事にする。
最初にN_k/D_kを計算し、それに対して、N_k-1/(D_k-1 + N_k/D_k)を計算する。これを計算し最後にN₁/(D₁+…)で計算完了。プログラムにすると

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

実行結果

> (/ 1 (cont-frac (lambda (i) 1.0)
                  (lambda (i) 1.0)
                  10))
1.6181818181818184
> (/ 1 (cont-frac (lambda (i) 1.0)
                  (lambda (i) 1.0)
                  100))
1.618033988749895
>

recursiveの時と同じ結果が出ている。

Exercise 1.38

オイラーが示したネイピア数の連分数表現との事。Diが意外と厄介。３つずつの規則性があるので３で割った商と余りを考えて：

i	1	2	3	4	5	6	7	8	9
D_i	1	2	1	1	4	1	1	6	1
i/3	0	0	1	1	1	2	2	2	3
i%3	1	2	0	1	2	0	1	2	0

単純には、iを3で割った余りが2の時、D_i=i - i / 3。それ以外は1。
プログラムにすると

> (+ 2 (cont-frac (lambda (i) 1.0)
                  (lambda (i) (if (= 2 (modulo i 3))
                                  (- i (quotient i 3))
                                  1))
                  100))
2.7182818284590455
>

ネイピア数は
e = 2.71828 18284 59045 23536 02874 71352 …
との事なのでまあまあ近似出来ていそう。

Exercise 1.39

ランバートによるtanの連分数展開。
N_iはi=1の時のみx、あとは-x²。マイナスに注意。
D_iは1、3、5…なので、2i - 1。
プログラムにすると

(define (tan-cf x k)
  (cont-frac (lambda (i) (if (= i 1)
                             x
                             (* -1 (* x x))))
             (lambda (i) (- (* 2 i) 1))
             k))

実行してみる。

> (tan-cf (/ 3.1415926535 4) 100)
0.9999999999551035
> (tan (/ 3.1415926535 4))
0.9999999999551034
> (tan-cf (/ 3.1415926535 6) 100)
0.5773502691696717
> (tan (/ 3.1415926535 6))
0.5773502691696718
>

そこそこ近似出来ていそう。