SICP 1.3.4 Procedures as Returned Values - プログラミング再入門

戻り値としての手続き。

ノート

Average dumpを使ったsqrtの定義は、

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

y ⇒ x/yと言う定義の形を崩してしまっている。(lambda (y) (/ x y))と言う形をそのまま使う為には、yを受け取ってf(y)を返す関数ではなく、この関数（f(y)）をyを引数に取りyとf(y)の中間点を返す関数に変換する関数を作り、それをfixed-pointに渡せば良い事になる。実行してみる。DrRacketではsqrtは再定義出来ない事になっているので、mysqrtとして定義した。

> (mysqrt 10)
3.162277660168379
> (mysqrt 5)
2.236067977499978
>

指摘の通り数学的な概念を極めて忠実にプログラムに変換出来ている。

cube-rootも実行してみる。

> (cube-root 8)
1.9999981824788517
> (cube-root 27)
2.9999972321057697
> (cube-root 64)
3.9999968907104337
> (cube-root 125)
4.999997738179904
>

Newton's Method

ニュートン法。
点(x₁, g(x₁))を通り、傾きg'(x₁)の直線は
$y - g(x_1) = g'(x_1)(x - x_1)$
と表せる。教科書ではg'(x)ではなくDg(x)と表しているので、
$y - g(x_1) = Dg(x_1)(x - x_1)$
この直線のx切片（x₂, 0)とすると
$0 - g(x_1) = Dg(x_1)(x_2 - x_1)$
$x_2 - x_1 = - g(x_1) / Dg(x_1)$
$x2 = x_1 - g(x_1) / Dg(x_1)$
x₂をx₁の関数と捉えて、f(x)で表すと
$f(x) = x - g(x) / Dg(x)$
このf(x)の不動点が近似解。

1.1.7ではDg(x)を数学的に求めてしまってプログラムを作っていた。ここではこの導関数を求める部分もプログラムで行う。本来概念的な極限値のdxを具体的に0.00001とか決めてしまって、Dg(x)を力技で求めてしまう。２点(x, g(x))と(x+dx, g(x + dx))を結ぶ直線の傾きを求める。x方向の増量はdx、y方向の増量はg(x + dx) - g(x)。
これをdxを大域の定数として関数gを引数に、xを与えるとDg(x)を返す関数に変換する関数が(deriv g)。

例をDrRacket上で実行してみる。f(x) = x³なので、f'(x)=3x²。x=5の時f'(5)=3×5²=75。

> (define (cube x) (* x x x))
> ((deriv cube) 5)
75.00014999664018
>

関数g(x)から関数g'(x)に変換する関数がnewton-transform。
関数g(x)とxの初期値guessを受け取り、newton-transformでg'(x)に変換して、fixed-pointを捜す関数がnewtons-method。
xの平方根を求めるにはg(y) = y² - xとした時のg(y)=0の解を求める事になる。

> (mysqrt2 5)
2.2360679775020436
> (sqrt 5)
2.23606797749979
>

Abstractions and first-class procedures

抽象と第一級手続き。
更にaverage-dumpとnewton-ftransforの二種類の変換も抽象化可能。
このfixed-point-of-transformを使ったsqrtの定義をmysqrt3とmysqrt4として定義して

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

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

実行してみる。

> (mysqrt3 5)
2.236067977499978
> (mysqrt4 5)
2.2360679775020436
> (sqrt 5)
2.23606797749979
>

手続きを引数として渡す、関数値として返す事で手続きが高度に抽象化出来る事を見て来た。必ずしも極限まで抽象化するべきと言う訳ではない。

より制限が少ない要素を第一級の状態を持つと表現する。それらは

変数として名前を付けられる
引数として関数に渡す事が出来る
関数値として返す事が出来る
データ構造の中に含める事が出来る

Lispは手続きを完全に第一級として扱っている。これはより効率的な実装環境の構築を難しくしているが、この事による表現力の増大は絶大である。

Exercise 1.40

x³+ax²+bx+cをそのまま関数として定義

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

ニュートン法でx³+ax²+bx+c=0となる点が求められる。検算としてx³-27=0の階を計算してみる。

> (newtons-method (cubic 0 0 -27) 1.0)
3.0000000000019176
>

Exercise 1.41

与えられた関数を２回適用する関数double。

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

例を実行してみる。

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

(double double)は2回適用する関数を2回適用するので、都合引数の関数を４回適用する。

ここまでは良い。ここからがややこしい。

(double (double double))はついincを４回適用して、更にincをもう４回適用する様なイメージを持ってしまうが、実はこれをincに適用する前にもう１回doubleを適用している。doubleは引数として与えられた関数を２回適用する関数を返すので(double double)に(double double)を適用した結果の関数を返す事になり、これにincを適用している。(double double)は引数として与えられた関数（ここではそれも(double double)）を４回適用する関数だから、与えられた引数を４回適用する関数を更に４回適用する関数を返して、それにincを適用している。つまりincを合計で16回適用するので5+16=21が答え。

勘違いしている方は４回適用する関数を２回適用すると勘違いしているのでプログラムにすると

> ((double ((double double) inc)) 5)
13
>

をイメージしてしまっている。括弧の場所が少しだけ違う。

Exercise 1.42

関数を合成する手続き。ひとつ引数を取る関数fとgを取り、gを先に適用して、その結果にfを適用する関数を生成する。

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

実行してみる。

> ((compose square inc) 6)
49
>

Exercise 1.43

関数fをn回適用する関数に変換する。nが1ならfだけ、それ以外ならrepeatedを再帰呼び出しして、Exercise 1.42のcomposeを使ってfを合成する。

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

実行してみる。

> ((repeated square 2) 5)
625
>

Exercise 1.44

関数fを取り、f(x)に対してf(x - dx)、f(x)、f(x + dx)の平均を値とする関数を生成する関数smoothを定義する。

(define (smooth f dx)
  (define (average3 a b c)
    (/ (+ a b c) 3))
  (lambda (x) (average3 (f (- x dx)) (f x) (f (+ x dx)))))

smoothとrepeatedを使ってn回smoothを掛けるn-fold smoothed functionを定義する。

repeatedするのはfではなくsmooth
repeatedに渡す関数はcomposeを使っているので、１引数である必要がある。
smoothの二つの引数のうちどちらかは部分適用しておく必要がある。
smooth自身は１引数の関数を返す。

(define (n-folded-smooth f dx n)
  (repeated (smooth f dx) n))

検証が難しい。

> (define pi 3.1415926535)
> ((n-folded-smooth sin (/ pi 4) 1) (/ pi 4))
0.569035593723558
> (sin (/ pi 4))
0.7071067811706742
> (/ (+ (sin 0) (sin (/ pi 4)) (sin (/ pi 2))) 3)
0.569035593723558
>

一応、意図した通りの平均値にはなっている。が、smoothを10回繰り返すと

> ((n-folded-smooth sin (/ pi 4) 10) (/ pi 4))
0.07031759412420362
>

smoothを繰り返した時にどうなるべきなのか今ひとつ理解出来ていない。

Exercise 1.45

プログラムは

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

averageの回数１回のとき

> (nth-root 2 4 1)
2.000000000000002
> (nth-root 3 8 1)
1.9999981824788517
> (nth-root 4 16 1)
. . user break
>

３乗根までは収束する。
averageの回数２回なら

> (nth-root 4 16 2)
2.0000000000021965
> (nth-root 5 32 2)
2.0000015129957607
> (nth-root 6 64 2)
2.0000029334662086
> (nth-root 7 128 2)
2.0000035538623377
> (nth-root 8 256 2)
. . user break
>

７乗根まで収束する。
average回数３回なら

> (nth-root 8 256 3)
2.0000000000039666
> (nth-root 9 512 3)
1.9999997106840102
> (nth-root 10 1024 3)
2.000001183010332
> (nth-root 11 2048 3)
1.999997600654736
> (nth-root 12 4096 3)
1.9999976914703093
> (nth-root 13 8192 3)
2.0000029085658984
> (nth-root 14 16384 3)
1.9999963265447058
> (nth-root 15 32768 3)
2.0000040951543028
> (nth-root 16 65536 3)
. . user break
>

15乗根まで。
３→７→１５は7=3×2+1、15=7×2+2と来ているので、次は15×2+1=31。average回数４回の場合31乗まで収束すると予想。

> (nth-root 31 (expt 2 31) 4)
1.9999951809750391
> (nth-root 32 (expt 2 32) 4)
. . user break
>

予想通り。次はaverage回数５回の場合は63乗まで収束すると予想。

> (power 2 63)
9223372036854775808
> (power 2 64)
18446744073709551616
> (nth-root 63 (expt 2 63) 5)
2.0000048698155544
> (nth-root 64 (expt 2 64) 5)
. . user break
>

予想通り。
averageの回数をa、計算出来る最大乗数をnとするとn=2^a+1-1の関係がある。
これをaについて解くと
$n=2^{a+1}-1$
$2^{a+1}=n+1$
$a+1=\log_2(n+1)$
$a=\frac{\log(n+1)}{\log2}-1$
これでaverageを適用する回数num-averageを決めるバージョンを書く。

> (ceiling (- (/ (log (+ 63 1)) (log 2)) 1))
5.0
> (ceiling (- (/ (log (+ 64 1)) (log 2)) 1))
6.0
>

63乗根はaverageを５回で良いが、64乗根で６回必要。これを使ってnth-root2を定義してみる。

(define (nth-root2 n x)
  (let ((num-average (ceiling (- (/ (log (+ n 1)) (log 2)) 1))))
    (fixed-point ((repeated average-damp num-average)
                  (lambda (y) (/ x (expt y (- n 1)))))
                 1.0)))

実行してみる。

> (nth-root2 2 4)
2.000000000000002
> (nth-root2 3 8)
1.9999981824788517
> (nth-root2 4 16)
2.0000000000021965
> (nth-root2 8 256)
2.0000000000039666
> (nth-root2 16 65536)
2.000000000076957
> (nth-root2 32 (expt 2 32))
2.000000000000006
> (nth-root2 64 (expt 2 64))
2.0000000000000853
>

Exercise 1.46

究極の抽象手続きiterative-improve。1.1.7のsqrt-iterも原型に出来るが、再帰の際に二つの値を渡している所が1.3.3と異なるので、1.3.3のfixed-pointを原型とする。close-enough?を引数で渡す様にし、iterative-improve自信は関数を生成して返す様に変形する。

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

close-enough?に渡す関数を作る関数を定義する。

(define (make-close-enough tolerance)
  (lambda (x1 x2) (< (abs (- x1 x2)) tolerance)))

1.7.7のsqrtを再定義すると

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

(define (sqrt-1.7.7 x)
  ((iterative-improve (make-close-enough 0.001)
                      (lambda (guess) (average guess (/ x guess))))
   1.0))

実行してみる。

> (sqrt-1.7.7 2)
1.4142135623746899
>

次に、1.3.3のfixed-pointを再定義する。fixed-pointを使ったsqrtの定義はそのまま。

(define (fixed-point f guess)
  ((iterative-improve (make-close-enough 0.001) f) guess))

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

実行してみる。

> (sqrt-1.3.3 2)
1.4142135623746899
>

ついでにaverage-dampを使った1.3.4のsqrtも再定義してみる。

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

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

実行してみる。

> (sqrt-1.3.4 2)
1.4142135623746899
>