プログラミング再入門

プログラミングをもう一度ちゃんと勉強する読書ノート

SICP 1.3 Formulating Abstractions with Higher-Order Procedures

ノート

そのまま訳すと『高階手続きによる抽象の定式化』か。単に定式化とするか抽象化とする方がピンときそう。

引数を使って具体的な数値あるいは変数から切り離して値の操作を手続きとして纏める事もひとつの抽象化であると。普通の言語は一連の操作のパターンに名前をつけて、以降はそれを使って計算出来る機構を備えているが、仮にその手続きが数値にしか適用出来ない(数値しか引数に出来ない)と、それは大きな制限となる。しばしば複数の手続きが共通のパターンを持つ事があり、これらを概念として抽象化する為には、引数としてあるいは関数値として手続き(あるいは関数)を受け渡しする必要がある。手続きを捜査する手続きを高階手続きと呼ぶ。

1.3.1 Procedures as Arguments

『引数としての手続き』

数学の\Bigsumの概念を手続きとして表す。ここではそれに加えてnがどう変化するかも、単に整数として1ずつ増えるだけではなく、それを制御する手続きも引数として取っている。

Σに渡すべき関数fと言う程ではなく単にaからbまでの数を足す場合には\Bigsum_{n=a}^{b}nで、この場合はf(n)=nとなる様な関数fを定義する。

nを整数ではなく小さな値(dx)で変化させると考えれば\Bigsum\Bigintに変形出来る。ここで紹介されている公式は長方形近似だが、最も単純なf(x)\times dxの足し算ではなくf(x+\frac{dx}{2})\times dxとなっている。精度としてどう変わるのかは良く分からない。

Exercise 1.29

長方形近似、台形近似よりも精度が高いシンプソンの公式による近似。正確には合成シンプソン公式と呼ばれるものらしい。

(define (simpson f a b n)
  (define (sim h)
    (define (c k)
      (cond ((or (= k 0) (= k n)) 1)
            ((odd? k) 4)
            (else 2)))
    (define (y k)
      (f (+ a (* k h))))
    (define (term k)
      (* (c k) (y k)))
    (* (/ h 3 ) (sum term 0 inc n)))
  (sim (/ (- b a) n)))

実行して長方形近似と比較してみる。

> (integral cube 0 1 0.01)
0.24998750000000042
> (simpson cube 0 1 100)
1/4
> (integral cube 0 1 0.001)
0.249999875000001
> (simpson cube 0 1 1000)
1/4
> 

シンプソン法の方は小数を扱っていないので分数の結果となってしまい、あまり上手く比較出来ている様な気はしない。
実際nが小さくてもそれなりの答えが出て来る。

> (simpson cube 0.0 1.0 2.0)
0.25
> 

a=0, b=1, n=2の時、h=\frac{1}{2}となるので
\frac{1}{2}\times\frac{1}{3}\times \left{1 \times (0 + 0 \times \frac{1}{2})^3 + 4 \times (0 + 1 \times \frac{1}{2})^3 + 1 \times (0 + 2 \times \frac{1}{2})^3\right}
=\frac{1}{6} \times \left(0 + 4 \times \frac{1}{8} + 1\right)
=\frac{1}{6} \times \frac{3}{2}
=\frac{1}{4}
なので、この条件ではそもそも近似ではなく正確な答えが出てしまう模様。

Exercise 1.30

1.2.1で出て来たlinear recursionとiterationの定義。再帰の帰り道で計算されるのをrecursionと呼び、末尾再帰になっていて再帰の行きで計算する事をiterationと呼んでいる。確かにsumの定義は+演算が残るので末尾再帰にはなっていない。
穴埋め問題を解くと

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

引数を特に何も加工せずに返す関数identicalを定義する。

(define (identical x) x)

簡単な検算。

> (sum identical 1 inc 10)
55
> 

更に検算のためExercise 1.29と同じ計算を実行してみる。

> (integral cube 0 1 0.01)
0.24998750000000042
> (simpson cube 0 1 100)
1/4
> (integral cube 0 1 0.001)
0.24999987500000073
> (simpson cube 0 1 1000)
1/4
> 

b=0.001のケースの計算結果がExercise 1.29の時と異なるのは足し算の順序が逆になるからだと思われる。再帰バージョンのsumの場合、再帰の巻き戻しの際に0.99953+0.99853+...+0.00153+0.00053と計算して行くが、繰り返しバージョンでは0.00053+0.00153+...+0.99853+0.99953と計算される。前者の方が結果の値が最初から大きいので、最後の方に足されるごく小さな値の端数が桁落ちする割合がより大きいのではないか。

Exercise 1.31

a.
sumと同様に掛け算をするproductを定義する。繰り返しタイプで作ると:

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

sumとは演算子がひとつ違うだけなので、当然『纏めよう』と言う話に繋がると思われる。

これを使ってfactorialを定義する。

(define (identical x) x)
(define (inc x) (+ x 1))

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

この本ではまだ登場していないが無名関数を使うと

(define (factorial2 n)
  (product (lambda (x) x) 1 (lambda (x) (+ x 1)) n))

検算

> (factorial 10)
3628800
> (factorial2 10)
3628800
> 

問題に示された漸化式を用いてπを計算する。
ウォリス積と呼ばれる公式らしい。各項は2/3, 4/3, 4/5…と言う数列。nextで2/3から4/3、4/3から4/5を計算する事も可能だが(法則は「分子か分母の小さい方に2を足す」)、終値が決められない。なので整数引数に対してf(1)=\frac{2}{3}f(2)=\frac{4}{3}...となる様な一般項を定義してtermに与える事にする。
分子はnに対して、nが奇数であれば+1、偶数であれば+2、分母は逆にnが奇数であれば+2、偶数であれば+1。

n 1 2 3 4 5 6
分子 2 4 4 6 6 8
分母 3 3 5 5 7 7

関数にしてみると

(define (term n)
  (if (odd? n)
      (/ (+ 1 n) (+ 2 n))
      (/ (+ 2 n) (+ 1 n))))

簡単に検算。

> (term 1)
2/3
> (term 2)
1 1/3
> (term 3)
4/5
> (term 4)
1 1/5
> 

(帯分数で表示されている)

これを使ってpiを計算する。incは普通に1足す物として

(define (inc n) (+ 1 n))

公式からpiを以下の様に定義する。掛け算する項の数を引数に取る。

(define (pi n) (* (product term 1 inc n) 4.0))

最後の4倍は整数の掛け算にしてしまうと分数で答えが表示されてしまうので、敢えて4.0として小数で表示される様にする。
検算

> (pi 100)
3.1570301764551676
> (pi 1000)
3.1431607055322663
> (pi 10000)
3.1417497057380523
> (pi 100000)
3.1416083612781764
> 

最後の100000項の計算は少し時間がかかる。

b.
最初に繰り返し版を書いたので、ここでは再帰版を定義する。

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

検算

> (pi 100)
3.1570301764551676
> (pi 1000)
3.1431607055322663
> (pi 10000)
3.1417497057380523
> (pi 100000)
. . user break
> 

結果は同じだが、こちらの方が計算が遅く100000では待てど暮らせど答えは返って来ない。

Exercise 1.32

a.
Exercise 1.31で予想した通り、sumとproductを一般化せよと。演算子+と*、それからそれぞれの演算の初期値、加算に対しては0、乗算に対しては1を引数に取り、sum、productと同様の演算をする。Schemeでは簡単に定義出来て:

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

これを使うとsumとproductはこうなる。

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

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

簡単な検算。

> (sum identical 1 inc 10)
55
> (product identical 1 inc 6)
720
> 

これでsimpsonとpiを計算してみる。

> (simpson cube 0 1 100)
1/4
> (pi 1000)
3.1431607055322663
> 

ちゃんと動いていそう。

b.
a.では繰り返し版を書いたので再帰版を書くと:

(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 identical 1 inc 10)
55
> (product identical 1 inc 6)
720
> (simpson cube 0 1 1000)
1/4
> (pi 1000)
3.1431607055322663
> 

動いた。

Exercise 1.33

範囲内の項について無条件に項を加算あるいは乗算するのではなく条件を満たすもののみを演算の対象とするfiltered-accumulate。計算の効率は良くないが構造が単純なrecursiveな定義を用いる。

2種類実装してみた。条件を満たした時と満たさなかった時のfiltered-accumulateを再帰呼び出しする部分が全く重複するのでnext-valueとして定義する。

(define (filtered-accumulate combiner null-value filter term a next b)
  (define (next-value) (filtered-accumulate combiner null-value filter term (next a) next b))
  (cond ((> a b) null-value)
        ((filter a) (combiner (term a) (next-value)))
        (else (next-value))))

ifも式である事を利用して、条件を満たさなかった場合には(term a)の代わりにnull-valueを演算に組み込む事でfiltered-accumulateの呼び出しを1回にした。

(define (filtered-accumulate combiner null-value filter term a next b)
  (if (> a b)
      null-value
      (combiner (if (filter a)
                    (term a)
                    null-value)
                (filtered-accumulate combiner null-value filter term (next a) next b))))

検算用に無条件に#tを返す関数trueを用意する。

(define (true a) #t)

(define (product term a next b)
  (filtered-accumulate * 1 true term a next b))

(define (sum term a next b)
  (filtered-accumulate + 0 true term a next b))

(define (identical x) x)
(define (inc n) (+ 1 n))

簡単に検算する。

> (sum identical 1 inc 10)
55
> (product identical 1 inc 6)
720
> 

a.
範囲内の素数の自乗和。1.2.6のprime?を使用。squareを定義して、1は素数に含めないので2以上で計算。

> (filtered-accumulate + 0 prime? square 2 inc 20)
1028
> 

2^2+3^2+5^2+7^2+11^2+13^2+17^2+19^2
=4+9+25+49+121+169+289+361
=1027
合ってそうかな。

b.
relatively primeは「互いに素」の事。1.2.5のgcdを使用。今度は色々と複雑なので関数として定義する。

(define (product-of-all-positive-intergers-relatively-prime-to n)
  (define (relatively-prime-to-n? a)
    (= 1 (gcd a n)))
  (define (gcd a b)
    (if (= b 0)
        a
        (gcd b (remainder a b))))
  (filtered-accumulate * 1 relatively-prime-to-n? identical 1 inc (- n 1)))

計算結果

> (product-of-all-positive-intergers-relatively-prime-to 2)
1
> (product-of-all-positive-intergers-relatively-prime-to 3)
2
> (product-of-all-positive-intergers-relatively-prime-to 4)
3
> (product-of-all-positive-intergers-relatively-prime-to 5)
24
> (product-of-all-positive-intergers-relatively-prime-to 6)
5
> (product-of-all-positive-intergers-relatively-prime-to 7)
720
> (product-of-all-positive-intergers-relatively-prime-to 8)
105
> (product-of-all-positive-intergers-relatively-prime-to 9)
2240
> (product-of-all-positive-intergers-relatively-prime-to 10)
189
> (product-of-all-positive-intergers-relatively-prime-to 11)
3628800
> 

1の扱いが難しいが、「GCDが1である」と言う定義には反していないし、答えは掛け算で特に影響が無いので1も対象とする事にすると、
2に対して互いに素なのは1だけなので、1。以降1は自明として…
3に対して互いに素なのは3だけなので、2。
4に対して互いに素なのは3だけなので、3。
5に対して互いに素なのは2、3、4なので、24。
6に対して互いに素なのは5だけなので、5。
7に対して互いに素なのは2、3、4、5、6なので、6!で720。
8に対して互いに素なのは3、5、7なので、105。
9に対して互いに素なのは2、4、5、7、8なので、2240。
10に対して互いに素なのは3、7、9なので189。
11に対して互いに素なのは10!なので3628800。
合ってそうかな。