[SICP] 問題 1.37 : 連分数

a. 無限の連分数(continued fraction)は
の形の式である. 例えばNiとDiがすべて1の無限連分数展開が1/φになることが示せる. φは(1.2.2節で示した)黄金比. 無限連分数の近似値のとり方の一つは, 与えられた項数で展開を中断することで, そういう中断--- k項有限連分数(k-term finite continued fraction)という---は
の形である. nとdを一引数(項の添字i)で連分数の項のNiとDiを返す手続きとする. (cont-frac n d k)がk項有限連分数を計算するような手続きcont-fracを定義せよ.

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

のkの順次の値で1/φの近似をとり, 手続きを調べよ. 4桁の精度の近似を得るのに, kはどのくらい大きくしなければならないか.

b. cont-fracが再帰的プロセスを生成するなら, 反復的プロセスを生成するものを書け. 反復的プロセスを生成するなら, 再帰的プロセスを生成するものを書け.

連分数を求める手続きを定義します.

;; 再帰的プロセスを生成
(define (cont-frac-r n d k)
  (define (iter i)
    (if (> i k)
        0
        (/ (n i)
           (+ (d i) (iter (+ i 1))))))
  (iter 1))

;; 反復的プロセスを生成
(define (cont-frac-i n d k)
  (define (iter i result)
    (if (= i 0)
        result
        (iter (- i 1)
              (/ (n i)
                 (+ (d i) result)))))
  (iter k 0))

(define (phi-r k)
  (/ 1 (cont-frac-r (lambda (i) 1.0)
                    (lambda (i) 1.0)
                    k)))

(define (phi-i k)
  (/ 1 (cont-frac-i (lambda (i) 1.0)
                    (lambda (i) 1.0)
                    k)))

再帰的プロセスを生成する手続きは, 上の方から計算しています. 反復的プロセスを生成する手続きは, 下の方から計算するようにしています.

黄金比Φの計算をしてみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (phi-r 15)
1.6180327868852458
> (phi-i 15)
1.6180327868852458
> 

kの値を15にすれば良さそうです.

[SICP] 問題 1.36 : x^x=1000

問題1.22で示した基本のnewlineとdisplayを使い, 生成する近似値を順に印字するようfixed-pointを修正せよ. 次にx log(1000)/log(x)の不動点を探索することで, xx = 1000の解を見つけよ. (自然対数を計算するSchemeの基本log手続きを使う.) 平均緩和を使った時と使わない時のステップ数を比べよ. ({ fixed-pointの予測値を1にして始めてはいけない. log(1)=0による除算を惹き起すからだ.)

近似値を表示するための処理をfixed-pointに追加しなくてはなりません. 追加する場所は, 内部手続tryが適切でしょう. 修正した手続きは次のとおりです.

(define (average x y)
  (/ (+ x y) 
     2))

(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)))
      (display guess)
      (newline)
      (if (close-enough? guess next)
          next
          (try next))))
  (try first-guess))

x^x=1000を満たすxを求めてみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fixed-point (lambda (x) (/ (log 1000) (log x))) 
               2.0)
2.0
9.965784284662087
3.004472209841214
6.279195757507157
3.759850702401539
5.215843784925895
4.182207192401397
4.8277650983445906
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
> (expt 2 3)
8
> (expt 4.555532270803653 4.555532270803653)
999.9913579312362
> 

期待した結果が得られているようです.

平均緩和を使ってみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fixed-point (lambda (x) (average x (/ (log 1000) (log x)))) 
               2.0)
2.0
5.9828921423310435
4.922168721308343
4.628224318195455
4.568346513136242
4.5577305909237005
4.555909809045131
4.555599411610624
4.5555465521473675
4.555537551999825
> 

ステップ数が1/3以下になっています.

[SICP] 問題 1.35 : 不動点の探索による黄金比の計算

(1.2.2節の)黄金比φが変換 x → 1 + 1/x の不動点であることを示し, この事実を使いfixed-point手続きによりφを計算せよ.

黄金比の定義は次のとおりです.
　　φ² = φ + 1
この両辺をφで割ります.
　　φ = 1 + 1/φ
この結果から,
　　変換 x → 1 + 1/x
が得られます.

不動点を求める手続きは次のとおりです.

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

実行してみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fixed-point (lambda (x) (+ 1 (/ 1 x)))
               1.0)
1.6180327868852458
> 

p.21から黄金比の値は約1.6180ですから, 求める結果が得られています.

[SICP] 問題 1.34 : (f f)の評価

手続き

(define (f g)
  (g 2))

を定義したとする. その時

(f square)
4

(f (lambda (z) (* z (+ z 1))))
6

解釈系に組合せ(f f)を(意地悪く)評価させるとどうなるか. 説明せよ.

問題文の手続きを定義して実行してみます. 手続きは次の通り.

(define (f g)
  (g 2))

(define (square x) (* x x))

実行結果は次のとおりです.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (f square)
4
> (f (lambda (z) (* z (+ z 1))))
6
> (f f)
. . application: not a procedure;
 expected a procedure that can be applied to arguments
  given: 2
  arguments...:
   2
> 

(f f)を評価するとエラーが生じています. その理由を考えます.

手続きfの定義から(f f)を評価した場合の置き換えは次のようになります.

  (f f)
= (f 2)
= (2 2)

「２」は手続きではないのでエラーが生じます.

[SICP] 1.2.6 : expmodでbaseのexp乗をmで割った余りが求まる理由

baseのexp乗をmで割った余りを求める手続きexpmodは次のような定義になります.

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))

法演算について馴染みがないと, なぜ, この手続で余りが求まるか分からないと思います. ここで簡単に計算の仕組みを説明します.

まず, 整数aをmで割った余りをr, 整数bをmで割った余りをtとします. すると整数aとbは次のように表現できます.

   a = qm + r
   b = sm + t

このaとbの和をmで割った余りを求めます.

   a + b = (qm + r) + (sm + t)
         = (q + s)m + (r + t)

(q + s)mはmで割り切れるので, a + bをmで割った余りはr+tをmで割った余りと等しくなります.

aにある整数cを掛けた値をmで割った余りも求めてます.

   ca = cqm + cr

caをmで割った余りはcrをmで割った余りと等しくなります.

続けて, aとbの積をmで割った余りを求めます.

   ab = (qm + r)(sm + t)
      = qsm^2 + (qt + rs)m + rt

abをmで割った余りはrtをmで割った余りと等しくなります.

積についての結果から,　aのx乗をmで割った余りはrのx乗をmで割った余りと等しくなります.

ここで, 手続きの定義に戻り, baseのexp乗をmで割った余りを考えます.

１）expが偶数のときを考えます.

baseの(exp/2)乗をmで割った余りをrとします. 上で計算した結果から, rの2乗をmで割った余りは, baseのexp乗をmで割った余りと等しくなります.

2) expが奇数のときを考えます.

baseの(exp-1)乗をmで割った余りをrとします. 上で計算した結果から, base*rをmで割った余りは, baseのexp乗をmで割った余りと等しくなります.

これらの結果から, 手続きexpmodによりbaseのexp乗をmで割った余りを求められることが分かります.

[SICP] 問題 1.27 : カーマイケル数

脚注47のCarmichael数が本当にFermatテストをだますことを示せ. それには整数nをとり, a < nなるすべてのaで, aⁿがnを法としてaの合同になるかどうか見る手続きを書き, Carmichael数でその手続きを使ってみる.

a < nを満たす正の整数aについて, expmodを用いてaⁿがnを法としてaと合同になるか調べます. 手続きは次のとおりです.

(define (square x) (* x x))

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))

(define (pass-fermat-test? n)
  (define (loop a)
    (or (= a n)
        (and (= (expmod a n n) a)
             (loop (+ a 1)))))
  (loop 1))

評価してみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (pass-fermat-test? 561)
#t
> (pass-fermat-test? 1105)
#t
> (pass-fermat-test? 1729)
#t
> (pass-fermat-test? 2465)
#t
> (pass-fermat-test? 2821)
#t
> (pass-fermat-test? 6601)
#t
> (pass-fermat-test? 6605)
#f
> (pass-fermat-test? 6607)
#t
> 

脚注のカーマイケル数は, フェルマー・テストをパスしています. 間違えて, 必ず#tを返す手続きを実装していては困りますので, 合成数6605では#fが返ってくることを確認しています.

[SICP] 問題 1.26 : expmodの計算に必要なステップ数

Louis Reasonerは問題1.24がなかなかうまくいかなかった. 彼の fast-prime?はprime?よりずっと遅いようであった. Louisは友人のEva Lu Atorに助けを求めた. Louisのプログラムを見ていると, squareを呼び出す代りに乗算を陽に使うように変っていた.

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (* (expmod base (/ exp 2) m)
                       (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))

「違いが分らん」とLouis, 「分る」とEvaがいった. 「手続きをこのように書き替えたので, Θ(log n)のプロセスをΘ(n)のプロセスにしてしまった.」 説明せよ.

引数expの値に注目して, 手続きが生成するプロセスを図示してみます. まず, squareを用いた場合.

13 -- 12 -- 6 -- 3 -- 2 -- 1 -- 0

まず, squareを用いず, 乗算を用いた場合.

13 -- 12 -- 6 -- 3 -- 2 -- 1 -- 0
       |    |         |
       |    |         + -- 1 -- 0
       |    |
       |    + -- 3 -- 2 -- 1 -- 0
       |              |
       |              + -- 1 -- 0
       |
       + -- 6 -- 3 -- 2 -- 1 -- 0
            |         |
            |         + -- 1 -- 0
            |
            + -- 3 -- 2 -- 1 -- 0
                      |
                      + -- 1 -- 0

このように, squareを用いずに乗算を用いた場合, expの値が偶数の場合, expmodを２回呼び出しています. そのため, Θ(log n)のプロセスをΘ(n)のプロセスにしてしまい, 計算に時間がかかるようになります.

[SICP] 問題 1.25 : expmodの必要性

Alyssa P. Hackerはexpmodを書くのに多くの余分なことをやったと不満で, 結局べき乗の計算法は知っているから

(define (expmod base exp m)
  (remainder (fast-expt base exp) m))

と単純に書ける筈だといった. これは正しいか. これも高速素数テストと同じに使えるか, 説明せよ.

テキストp.29の脚注48にあるように, ここでは200桁の整数が素数であるかどうかを問題としています. expmodの引数であるbaseやexpに200桁の整数が与えられるため, fast-exptで計算することは困難であり現実的ではありません. そのため, p.28のexpmodの定義が必要となります.

[SICP] 問題 1.24 : フェルマー・テスト

問題1.22のtimed-prime-test手続きをfast-prime?(Fermat法参照) を使うように変え, その問題で見つけた12個の素数をテストせよ. FermatテストはΘ(log n)の増加なので, 1,000,000近くの素数をテストする時間を1000近くの素数をテストする時間と比べ, どの位と予想するか. データはその予想を支持するか. 違いがあれば説明出来るか.

フェルマーの小定理を用いた素数性のテストをする手続きを定義します. 手続きrandomの引数の範囲に制限があるため, 少し小細工をしています.

(define (square x) (* x x))

(define (expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (expmod base (- exp 1) m))
                    m))))  

(define (rnd n)
  (remainder (* (random (remainder n 4294967087))
                (random (remainder n 4294967087)))
             n))

(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (rnd (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) #t)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (else #f)))

;; 素数を求めるために必要な時間を計測する
(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (current-milliseconds)))

(define (start-prime-test n start-time)
  (if (fast-prime? n 10000)
      (report-prime (- (current-milliseconds) start-time))))

(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

(define (search-for-primes start end)
  (cond ((< start end)
         (timed-prime-test start)
         (search-for-primes (+ start 1) end))
        (else 
         (newline)
         (display "end"))))
 

素数を求めてみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (search-for-primes 100000000000000 100000000000100)
100000000000031 *** 325
100000000000067 *** 344
100000000000097 *** 312
end
> (search-for-primes 1000000000000000 1000000000000200)
1000000000000037 *** 342
1000000000000091 *** 351
1000000000000159 *** 340
end
> (search-for-primes 10000000000000000 10000000000000100)
10000000000000061 *** 390
10000000000000069 *** 375
10000000000000079 *** 366
end
> 

この程度のデータでは, 何も言えそうにありません.

[SICP] 問題 1.23 : 素数性を判定する数を奇数に限定

本節の初めのsmallest-divisorは多くの不要な計算をする: 数が2で割り切れるか調べた後は, より大きい偶数で割り切れるか調べるのは無意味である. test-divisorが使う値は, 2, 3, 4, 5, 6, ...ではなく, 2, 3, 5, 7, 9, ...であるべきだ. この変更を実装するため, 入力が2なら3 を返し, そうでなければ入力に2足したものを返す手続きnextを定義せよ. smallest-divisorを(+ test-divisor 1)ではなく, (next test-divisor)を使うように修正せよ. このように修正した smallest-divisorにしたtimed-prime-testで, 問題1.22で見つけた12 個の素数をテストせよ. この修正はテストのステップを半分にしたので, 二倍速く走ることが期待される. この期待は確められるか. そうでなければ, 得られた二つのアルゴリズムの速度の比はどうであったか. それが2と違う事情を説明せよ.

問題文にあるように手続きを修正します.

(define (square x) (* x x))

(define (divides? a b)
  (= (remainder b a) 0))

;; 問題文により追加
(define (next n)
  (if (= n 2)
      3
      (+ n 2)))

(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (next test-divisor)))))

(define (smallest-divisor n)
  (find-divisor n 2))

(define (prime? n)
  (= n (smallest-divisor n)))

;; 素数を求めるために必要な時間を計測する
(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (current-milliseconds)))

(define (start-prime-test n start-time)
  (if (prime? n)
      (report-prime (- (current-milliseconds) start-time))))

(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

(define (search-for-primes start end)
  (cond ((< start end)
         (timed-prime-test start)
         (search-for-primes (+ start 1) end))
        (else 
         (newline)
         (display "end"))))
 

問題 1.22で求めた12個の素数を判定するための時間を計測します.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (search-for-primes 100000000000000 100000000000100)
100000000000031 *** 624  (939)
100000000000067 *** 571  (1397)
100000000000097 *** 598  (949)
end
> (search-for-primes 1000000000000000 1000000000000200)
1000000000000037 *** 1907  (2891)
1000000000000091 *** 1965  (2928)
1000000000000159 *** 1877  (2919)
end
> (search-for-primes 10000000000000000 10000000000000100)
10000000000000061 *** 5925  (9149)
10000000000000069 *** 5887  (9166)
10000000000000079 *** 5887  (9174)
end
> 

カッコ内に問題 1.22で計測した時間を書いています. 問題 1.22の結果と見比べると, 時間は短くなっていますが, 半分とまでは言えません. 手続きnextでnが2であるかどうかを判定しているからでしょう.

[SICP] 問題 1.22 : 素数を求めるために必要な計算量

殆んどのLispの実装には基本手続きruntimeがあり, システムが走った時間を(例えばマイクロ秒で)示す整数を返す. 次のtimed-prime-test手続きは整数nで呼び出されると, nを印字し, nが素数かどうかを調べ, nが素数なら手続きは三つの星印とテストに要した時間を印字する.

(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (runtime)))

(define (start-prime-test n start-time)
  (if (prime? n)
      (report-prime (- (runtime) start-time))))

(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

この手続きを使い, 指定範囲の連続する奇整数について素数性を調べる手続きsearch-for-primesを書け. その手続きで, それぞれ1,000, 10,000, 100,000, 1,000,000より大きい, 小さい方から三つの素数を見つけよ. 素数をテストする時間に注意せよ. テストのアルゴリズムはΘ(√n)の増加の程度だから, 10,000前後の素数のテストは1,000前後の素数のテストの√10倍かかると考えよ. 時間のデータはこれを支持するか. 100,000, 1,000,000のデータは√n予想のとおりだろうか. 結果はプログラムが計算に必要なステップ数に比例した時間で走るという考えに合っているか.

手続きruntimeの代わりにcurrent-millisecondsを使います. 手続きは次のようになります.

(define (square x) (* x x))

(define (divides? a b)
  (= (remainder b a) 0))

(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (+ test-divisor 1)))))

(define (smallest-divisor n)
  (find-divisor n 2))

(define (prime? n)
  (= n (smallest-divisor n)))

;; 素数を求めるために必要な時間を計測する
(define (timed-prime-test n)
  (newline)
  (display n)
  (start-prime-test n (current-milliseconds)))

(define (start-prime-test n start-time)
  (if (prime? n)
      (report-prime (- (current-milliseconds) start-time))))

(define (report-prime elapsed-time)
  (display " *** ")
  (display elapsed-time))

(define (search-for-primes start end)
  (cond ((< start end)
         (timed-prime-test start)
         (search-for-primes (+ start 1) end))
        (else 
         (newline)
         (display "end"))))

問題文では, 素数を求める範囲を1000から始めていますが, 時間を計測できないため, もっと大きな数字から始めます. 測定結果は次のとおりです. 見やすいように, 不要な結果を削除しています.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (search-for-primes 100000000000 100000000100)
100000000003 *** 45
100000000019 *** 29
100000000057 *** 28
100000000063 *** 29
end
> (search-for-primes 1000000000000 1000000000100)
1000000000039 *** 97
1000000000061 *** 90
1000000000063 *** 90
end
> (search-for-primes 10000000000000 10000000000100)
10000000000037 *** 293
10000000000051 *** 302
10000000000099 *** 307
end
> (search-for-primes 100000000000000 100000000000100)
100000000000031 *** 939
100000000000067 *** 1397
100000000000097 *** 949
end
> (search-for-primes 1000000000000000 1000000000000200)
1000000000000037 *** 2891
1000000000000091 *** 2928
1000000000000159 *** 2919
end
> (search-for-primes 10000000000000000 10000000000000100)
10000000000000061 *** 9149
10000000000000069 *** 9166
10000000000000079 *** 9174
end

結果を見ると, nを１桁増やすと時間は約３倍になっており、 nを２桁増やすと時間は約10倍になっています. このことから, プログラムが計算に必要なステップ数に比例した時間で走るという考えに合っていると言えます.

問題文には, search-for-primesは「指定範囲の連続する奇整数について」とありますが, この条件を満たしていません. 実害はないのでこのままとします.

[SICP] 問題 1.21 : 最小除数の探索

smallest-divisor手続きを使い, 次の数の最小除数を見つけよ: 199, 1999, 19999.

テキストの手続きを入力します.

(define (square x) (* x x))

(define (divides? a b)
  (= (remainder b a) 0))

(define (find-divisor n test-divisor)
  (cond ((> (square test-divisor) n) n)
        ((divides? test-divisor n) test-divisor)
        (else (find-divisor n (+ test-divisor 1)))))

(define (smallest-divisor n)
  (find-divisor n 2))

(define (prime? n)
  (= n (smallest-divisor n)))

評価してみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (smallest-divisor 199)
199
> (smallest-divisor 1999)
1999
> (smallest-divisor 19999)
7
> 

[SICP] 問題 1.19 : Fibonacci数を対数的ステップ数で計算する

Fibonacci数を対数的ステップ数で計算するうまいアルゴリズムがある. 1.2.2 節のfib-iterプロセスで使った状態変数aとbの変換: a ← a + bとb ← aに注意しよう. この変換をTと呼ぶ. 1と0から始め, Tを繰り返してn回作用させると, Fib(n + 1)とFib(n)の対が出来る. いいかえれば, Fibonacci数は対(1, 0)にTⁿ, つまり変換Tのn乗を作用させると得られる. さて, T_pqは対(a, b)をa ← bq + aq + apとb ← bp + aqに従って変換するものとし, Tを変換族T_pqのp = 0とq = 1の特例としよう. 変換T_pqを二回使うとその効果は同じ形式の変換T_p'q'を一回使ったのと同じになることを示し, p'とq'をp, qを使って表せ. これで変換を二乗する方法が得られ, fast-exptのように逐次平方を使い, Tⁿを計算することが出来る. これらをすべてまとめ, 対数的ステップ数の以下の手続きを完成せよ.

(define (fib n)
  (fib-iter 1 0 0 1 n))

(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                   ＜??＞      ; p'を計算
                   ＜??＞      ; q'を計算
                   (/ count 2)))
        (else (fib-iter (+ (* b q) (* a q) (* a p))
                        (+ (* b p) (* a q))
                        p
                        q
                        (- count 1)))))

変換T_pqは次のような行列として表現できます.

Tpq=	[	p+q	q	]
		q	p

変換T_p'q'はT_pq²となります. T_pq²を計算し整理すると次のようになります.

Tp'q'=	[	(p2+q2) + (2pq+q2)	2pq+q2	]
		2pq+q2	p2+q2

変換T_p'q'とT_pqを見比べると, p'とq'はpとqを用いて次のように表せます.

p'=	p2+q2
q'=	2pq+q2

上で求めたp'とq'を使って手続きを定義すると次のようになります. 計算結果を検証するために, p.22の手続きもfib0として定義します.

(require (lib "racket/trace.ss"))

(define (fib0 n)
  (define (iter a b count)
    (if (= count 0)
        b
        (iter (+ a b) a (- count 1))))
  (iter 1 0 n))

(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                   (+ (* p p) (* q q))
                   (+ (* 2 p q) (* q q))
                   (/ count 2)))
        (else
         (fib-iter (+ (* b q) (* a q) (* a p))
                   (+ (* b p) (* a q))
                   p
                   q
                   (- count 1)))))

(define (fib n)
  (fib-iter 1 0 0 1 n))

;;(trace fib-iter)

計算結果が正しいことを確認します.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (map fib0 '(0 1 2 3 4 5 6 7))
(0 1 1 2 3 5 8 13)
> (map fib '(0 1 2 3 4 5 6 7))
(0 1 1 2 3 5 8 13)
> 

トレースしてみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fib 30)
>(fib-iter 1 0 0 1 30)
>(fib-iter 1 0 1 1 15)
>(fib-iter 2 1 1 1 14)
>(fib-iter 2 1 2 3 7)
>(fib-iter 13 8 2 3 6)
>(fib-iter 13 8 13 21 3)
>(fib-iter 610 377 13 21 2)
>(fib-iter 610 377 610 987 1)
>(fib-iter 1346269 832040 610 987 0)
<832040
832040
> 

トレースの結果を見ると, 対数的ステップ数で計算していることが分かります.

[SICP] 問題 1.18 : 対数的ステップ数の乗算手続き（反復的プロセス）

問題1.16, 1.17の結果を使い, 加算, 二倍, 二分による, 対数的ステップ数の, 二つの整数の乗算の反復的プロセスを生成する手続きを工夫せよ.

問題 1.16の結果を参考にして手続きを定義します. 0や負の数をかけた場合の扱いは, ループに入る前に判定してます.

(require (lib "racket/trace.ss"))

(define (double x) (* x 2))

(define (halve x) (/ x 2))

(define (fast-mult-iter a b p)
  (cond ((= b 0) p)
        ((even? b) (fast-mult-iter (double a) (halve b) p))
        (else (fast-mult-iter a (- b 1) (+ a p)))))

(define (fast-mult a b)
  (cond ((< b 0) (* -1 (fast-mult-iter a (abs b) 0)))
        ((= b 0) 0)
        (else (fast-mult-iter a b 0))))

(trace fast-mult-iter)

p+a×bを考え, pに計算の途中経過をため込むことにします. pを0から始めます. p + a×bの計算について, 奇数と偶数の場合に分けて, 計算をどのように進めるかを考えます.

1) bが奇数の場合, b = m+1とする. 次の式が成立する.

p + a×(m+1) = (p+a) + (a×m)

(p+a)が次のpとなります.
2) nが偶数の場合, b = 2mとする. 次の式が成立する.

P + a×(2m) = p + double(a)×m

double(a)が次のaとなります.

実行結果は次のとおりです.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fast-mult 3 45)
>(fast-mult-iter 3 45 0)
>(fast-mult-iter 3 44 3)
>(fast-mult-iter 6 22 3)
>(fast-mult-iter 12 11 3)
>(fast-mult-iter 12 10 15)
>(fast-mult-iter 24 5 15)
>(fast-mult-iter 24 4 39)
>(fast-mult-iter 48 2 39)
>(fast-mult-iter 96 1 39)
>(fast-mult-iter 96 0 135)
<135
135
> (fast-mult 4 0)
0
> (fast-mult 5 -8)
>(fast-mult-iter 5 8 0)
>(fast-mult-iter 10 4 0)
>(fast-mult-iter 20 2 0)
>(fast-mult-iter 40 1 0)
>(fast-mult-iter 40 0 40)
<40
-40
> 

トレースの結果から, 反復的プロセスを生成していることが分かります.

[SICP] 問題 1.29 : Simpsonの公式

Simpsonの公式は上に示したのより更に正確な数値積分法である.
Simpsonの公式を使えば, aからbまでの関数fの積分は, 偶整数nに対し, h = (b - a)/nまたy_k = f(a + kh)として

(h/3)[y₀ + 4y₁ + 2y₂ + 4y₃ + 2y₄ + ・・・ + 2y_n-2 + 4y_n-1 + y_n]

で近似出来る. (nの増加は近似の精度を増加する.)
引数としてf, a, bとnをとり, Simpson公式を使って計算した積分値を返す手続きを定義せよ. その手続きを使って(n = 100とn = 1000で)0から1までcubeを積分し, またその結果を上のintegral手続きの結果と比較せよ.

Simpsonの公式の大カッコの中身をじっと見つめると, 次のようなパターンが見えてきます.

(y₀ + 4y₁ + y₂) + (y₂ + 4y₃ + y₄) + ・・・ + (y_n-4 + 4y_n-3 + y_n-2) + (y_n-2 + 4y_n-1 + y_n)

y₂に注目してください. 公式では係数が2となっていますが, それを2つに分けています.

つまり, 大カッコの中は３つの項をひとまとめにするとスッキリと書けます. 手続きを定義します.

(require (lib "racket/trace.ss"))

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

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

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

(define (simpson f a b n)
  (define h (/ (- b a) n))
  (define (y k) (f (+ a (* k h))))
  (define (term k) (+ (y k) (* 4.0 (y (+ k 1))) (y (+ k 2))))
  (define (next k) (+ k 2))
  (* (/ h 3.0)
     (sum term 0 next (- n 2))))

0から1までのcubeを積分してみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (integral cube 0 1 0.01)
0.24998750000000042
> (integral cube 0 1 0.001)
0.249999875000001
> (simpson cube 0 1 100)
0.2500000000000001
> (simpson cube 0 1 1000)
0.2500000000000001
> 

計算結果を見ると精度が全く異なっています.

[SICP] 問題 1.17 : 対数的ステップ数の乗算手続き（再帰的プロセス）

本節のべき乗アルゴリズムは, 乗算の繰返しによるべき乗の実行に基づいていた. 同様に整数の乗算を加算の繰返しで実行出来る. 次の乗算手続きは(この言語には加算はあるが, 乗算はないと仮定する)expt手続きに似たものである:

(define (* a b)
  (if (= b 0)
      0
      (+ a (* a (- b 1)))))

このアルゴリズムはbに線形のステップ数をとる. 加算の他に整数を二倍するdouble演算と(偶数の)整数を2で割るhalve演算もあるとしよう. これらを使い, fast-exptと類似の対数的ステップ数の乗算手続きを設計せよ.

p.25のfast-exptを参考に手続きを定義します.

(require (lib "racket/trace.ss"))

(define (double x) (* x 2))

(define (halve x) (/ x 2))

(define (fast-mult a b)
  (cond ((< b 0) (* -1 (fast-mult a (abs b))))
        ((= b 0) 0)
        ((= b 1) a)
        ((even? b) (double (fast-mult a (halve b))))
        (else (+ a (fast-mult a (- b 1))))))

(trace fast-mult)

fast-exptと同じような構造をしています. 負の数をかけた場合と０をかけた場合の扱いを追加しています.

a×bの計算について, 奇数と偶数の場合に分けて, 計算をどのように進めるかを考えます.

1) bが奇数の場合, b = m+1とする. 次の式が成立する.

a×(m+1) = a + (a×m)

2) nが偶数の場合, b = 2mとする. 次の式が成立する.

a×(2m) = double(a×m)

実行結果は次のとおりです.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fast-mult 3 45)
>(fast-mult 3 45)
> (fast-mult 3 44)
> >(fast-mult 3 22)
> > (fast-mult 3 11)
> > >(fast-mult 3 10)
> > > (fast-mult 3 5)
> > > >(fast-mult 3 4)
> > > > (fast-mult 3 2)
> > > > >(fast-mult 3 1)
< < < < <3
< < < < 6
< < < <12
< < < 15
< < <30
< < 33
< <66
< 132
<135
135
> (fast-mult 4 0)
>(fast-mult 4 0)
<0
0
> (fast-mult 5 -8)
>(fast-mult 5 -8)
> (fast-mult 5 8)
> >(fast-mult 5 4)
> > (fast-mult 5 2)
> > >(fast-mult 5 1)
< < <5
< < 10
< <20
< 40
<-40
-40
> 

トレースの結果から、対数的ステップ数で計算していることが分かります.

[SICP] 問題 1.16 : べき乗を求める反復的アルゴリズム

fast-exptのように, 逐次平方を使い, 対数的ステップ数の反復的べき乗プロセスを生成する手続きを設計せよ.

(ヒント: (b^n/2)² = (b²)^n/2に注意し, 指数nと底bの他にもう一つの状態変数aを用意する.
状態の移変りで積abⁿが不変であるように状態遷移を定義する.
プロセス開始時にaを1とし, プロセス終了時にaが結果になるようにする.
一般に, 状態の移変りに不変のままの不変量(invariant quantity)を定義する技法は, 反復的アルゴリズムの設計に有用である.)

問題文のヒントを参考に手続きを定義します.

(require (lib "racket/trace.ss"))

(define (square x) (* x x))

(define (fast-expt-iter b n p)
  (cond ((= n 0) p)
        ((even? n) (fast-expt-iter (square b) (/ n 2) p))
        (else (fast-expt-iter b (- n 1) (* b p)))))

(define (fast-expt b n)
  (fast-expt-iter b n 1))

(trace fast-expt-iter)

p×bⁿを考え, pに結果をため込むようにします. 奇数と偶数の場合に分けて, pとbがどのように変化するかを考えます.

1) nが奇数の場合, n = m+1とする. 次の式が成立する.

p×b^m+1 = (p×b)×b^m

右辺のカッコ内の(p×b)が次のステップのpとなる.

2) nが偶数の場合, n = 2mとする. 次の式が成立する.

p×b^2m = p×(b²)^m

右辺のカッコ内の(b²)が次のステップのbとなる.

実行結果は次のとおりです.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (fast-expt 2 8)
>(fast-expt-iter 2 8 1)
>(fast-expt-iter 4 4 1)
>(fast-expt-iter 16 2 1)
>(fast-expt-iter 256 1 1)
>(fast-expt-iter 256 0 256)
<256
256
> 

トレースの結果から, 反復的プロセスを生成しており, 対数的ステップ数で計算していることが分かります.

[SICP] 問題 1.15 : 正弦の計算

(ラジアンで表す)角度の正弦は, xが十分小さい時, sin x ≒ xの近似と, 正弦の引数の値を小さくするための三角関係式

sin(x) = 3sin(x/3) - 4sin³(x/3)

を使って計算出来る. (この問題のためには, 角度の大きさが0.1ラジアンより大きくなければ「十分小さい」と考える.) この方法は次の手続きに採用してある:

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

(define (p x) (- (* 3 x) (* 4 (cube x))))

(define (sine angle)
   (if (not (> (abs angle) 0.1))
       angle
       (p (sine (/ angle 3.0)))))

a. (sine 12.15)の評価で, 手続きpは何回作用させられたか.

b. (sine a)の評価で, 手続きsineの生成するプロセスが使うスペースとステップ数の増加の程度は(aの関数として)何か.

問題文の手続きを定義します. 手続きpを作用させた回数を調べるためにトレースを使います.

(require (lib "racket/trace.ss"))

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

(define (p x) (- (* 3 x) (* 4 (cube x))))

(define (sine angle)
  (if (not (> (abs angle) 0.1))
      angle
      (p (sine (/ angle 3.0)))))

(trace p)

計算結果は次のとおりです.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (sin 12.5)
-0.06632189735120068
> (sine 12.5)
>(p 0.051440329218107005)
<0.153776521096694
>(p 0.153776521096694)
<0.4467840153484446
>(p 0.4467840153484446)
<0.9836111719853284
>(p 0.9836111719853284)
<-0.8557060643295382
>(p -0.8557060643295382)
<-0.060813768577286265
-0.060813768577286265
> 

トレースの結果から, pを作用させた回数は５回です.

手続きの定義から, aの値が0.1より大きい時は手続きを再帰的に呼び出し, 次のaの値は1/3になっています. 「十分小さい」と判断する値をe, ステップ数をnとすると次の関係が成り立ちます.

(a/3)n	<	e
a/e	<	3n
log3(a/e)	<	n

このことから, ステップ数の増加の程度はlog aであるといえます. また, 手続きは末尾再帰になっていないため, スペースの増加の程度もステップの増加の程度と同じです.

[SICP] 問題 1.14 : 両替の計算が生成するプロセス

1.2.2節のcount-change手続きが11セントの両替の場合に生成するプロセスを示す木構造を描け. 両替の金額の増加につれ, このプロセスが使うスペースとステップ数の増加の程度は何か.

図を書きます. 木の節点に金額と両替に使えるコインを括弧付きで描いています.

 11 -----------------  -39
(50, 25, 10, 5, 1)     (50, 25, 10, 5, 1)
  |
  |
 11 -----------------  -14
(25, 10, 5, 1)         (25, 10, 5, 1)
  |
  |
 11 --------- 1 ------  -9
(10, 5, 1)  (10, 5, 1) (10, 5, 1)
  |           |
  |           |
  |           1 ------  -4
  |          (5, 1)     (5, 1)
  |           |
  |           |
  |           1  ------  0
  |          (1)        (1)
  |
 11 --------  6 -------  1 ----  -4
 (5, 1)      (5, 1)     (5, 1)   (5, 1)
  |           |          |
  |           |          |
  |           |          1 -----  0
  |           |         (1)      (1)
  |           |          |
  |           |          |
  |           |          1
  |           |         ( )
  |           |
  |           6 -- 5 -- 4 -- 3 -- 2 -- 1 -- 0
  |          (1)  (1)  (1)  (1)  (1)  (1)  (1)
  |           |    |    |    |    |    |
  |           |    |    |    |    |    |
  |           6    5    4    3    2    1
  |          ( )  ( )  ( )  ( )  ( )  ( )
  | 
 11 -- 10 -- 9 -- 8 -- 7 -- 6 -- 5 -- 4 -- 3 -- 2 -- 1 -- 0
 (1)   (1)  (1)  (1)  (1)  (1)  (1)  (1)  (1)  (1)  (1)  (1)
  |     |    |    |    |    |    |    |    |    |    |
  |     |    |    |    |    |    |    |    |    |    |
 11    10    9    8    7    6    5    4    3    2    1
 ( )   ( )  ( )  ( )  ( )  ( )  ( )  ( )  ( )  ( )  ( )

このプロセスが使うスペースとステップ数の増加の程度については後で考えます.

[SICP] 問題 1.13 : フィボナッチ数を求める代数方程式

Φ= (1+√5)/2としてFib(n)がΦⁿ/√5に最も近い整数であることを証明せよ. ヒント: Ψ= (1-√5)/2とする. 帰納法とFibonacci数の定義(1.2.2節参照)を用い, Fib(n)=(Φⁿ-Ψⁿ)/√5を証明せよ.

まず, Fib(n)=(Φⁿ-Ψⁿ)/√5を証明する.

(1) n=1のとき

Fib(1)	=	{Φ1-Ψ1} / √5
	=	{(1+√5)/2 - (1-√5)/2} / √5
	=	2√5 / √5
	=	1

(2) n=2のとき

Φ2	=	{(1+√5)/2}2
	=	{1 + 2√5 + 5} / 4
	=	{2 + 2√5 + 4} / 4
	=	{1 + √5} / 2 + 1
	=	Φ + 1

Ψ2	=	{(1-√5)/2}2
	=	{1 - 2√5 + 5} / 4
	=	{2 - 2√5 + 4} / 4
	=	{1 - √5} / 2 + 1
	=	Ψ + 1

Fib(2)	=	{Φ2-Ψ2} / √5
	=	{(Φ + 1) - (Ψ + 1)} / √5
	=	{Φ - Ψ} / √5
	=	1

(3) n=k, k+1のときFib(n)=(Φⁿ-Ψⁿ)/√5が成立すると仮定する. n=k+2のとき

Fib(k+2)	=	Fib(k+1) + Fib(k)
	=	(Φk+1-Ψk+1)/√5 + (Φk-Ψk)/√5
	=	{(Φk+1+Φk) - (Ψk+1+Ψk)} / √5
	=	{Φk(Φ+1) - Ψk(Ψ+1)} / √5
	=	{ΦkΦ2 - ΨkΨ2} / √5
	=	{Φk+2 - Ψk+2} / √5

(4) 以上より, 1以上の整数nについて, Fib(n)=(Φⁿ-Ψⁿ)/√5が成立する.

次に, Ψ/√5の値を計算する.

Ψ/√5 = (1-√5)/2 ≒ -0.2764

Ψ/√5 < 1/2 であるから, 1以上の整数nについて, Ψⁿ/√5 < 1/2 である.
よって、Fib(n)がΦⁿ/√5に最も近い整数であるといえる.