[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)のプロセスにしてしまい, 計算に時間がかかるようになります.