十と三: SICP

ラベル SICP の投稿を表示しています。すべての投稿を表示

2014年12月6日土曜日

[SICP] 問題 1.40 : Newton法

newtons-methodの手続きと一緒に

(newtons-method (cubic a b c) 1)

の形の式で使い, 三次式 x³ + ax² + bx + c の零点を近似する手続き cubicを定義せよ.

Newton法のための手続きと問題文で指定されたcubic手続きを定義します.

(define dx 0.00001)

(define tolerance 0.00001)

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

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

(define (deriv g)
  (lambda (x)
    (/ (- (g (+ x dx)) (g x))
       dx)))

(define (newton-transform g)
  (lambda (x)
    (- x (/ (g x) ((deriv g) x)))))

(define (newtons-method g guess)
  (fixed-point (newton-transform g) guess))

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

(newtons-method (cubic 0 0 -10) 1)

実行例として10の三乗根を求めてみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (newtons-method (cubic 0 0 -10) 1)
2.154434690031893
> (expt 2.154434690031893 3)
10.000000000000131
>

2014年12月4日木曜日

[SICP] 問題 1.39 : 正弦関数の連分数展開

正接関数の連分数展開は1770年にドイツの数学者 J. H. Lambertが発表した. xをラジアンで表し,
Lambertの式に基づいて正接関数の近似値を計算する手続き(tan-cf x k) を定義せよ. kは問題1.37と同様, 計算する項数を指定する.

問題文の式から手続きを定義します. 分子がx²となる項からcont-fracを使います.

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

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

(define (tan-cf x k)
  (/ x
     (+ 1 (cont-frac (lambda (i) (- (square x)))
                     (lambda (i) (+ (* i 2) 1))
                     k))))

tanは定義済みですので, その結果と比較します.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (tan 1.0)
1.5574077246549023
> (tan-cf 1.0 10)
1.557407724654902
> (tan-cf 1.0 100)
1.557407724654902
>

2014年12月2日火曜日

[SICP] 問題 1.38 : 連分数による自然対数の底eの近似

1737年, スイスの数学者 Leonhard Eulerは De Fractionibus Continuisというメモを発表した. その中にeを自然対数の底としてe - 2 の連分数展開がある. この分数ではNiはすべて1, Diは順に1, 2, 1, 1, 4, 1, 1, 6, 1, 1, 8, ... . 問題1.37のcont-fracを使い, Eulerの展開によりeを近似するプログラムを書け.

D_iの値は, iを3で割った余りが2のときに2, 4, 6, 8と増えています. 手続きを定義すると次のようになります.

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

(define (euler k)
  (+ 2
     (cont-frac (lambda (i) 1.0)
                (lambda (i) 
                  (if (= (remainder i 3) 2)
                      (* 2.0 (/ (+ i 1) 3))
                      1.0))
                k)))

実行してみます.自然対数の底は定義済みでしたので,その結果と比較します.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> e
2.718281828459045
> (euler 10)
2.7182817182817183
> (euler 100)
2.7182818284590455
>

それらしい値が得られています.

2014年11月29日土曜日

[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にすれば良さそうです.

2014年11月28日金曜日

[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以下になっています.

2014年11月27日木曜日

[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ですから, 求める結果が得られています.

2014年11月26日水曜日

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

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

2014年11月25日火曜日

[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で割った余りを求められることが分かります.

2014年11月22日土曜日

[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が返ってくることを確認しています.

2014年11月21日金曜日

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

2014年11月20日木曜日

[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の定義が必要となります.

2014年11月19日水曜日

[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
>

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

2014年11月17日月曜日

[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であるかどうかを判定しているからでしょう.

2014年11月15日土曜日

[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は「指定範囲の連続する奇整数について」とありますが, この条件を満たしていません. 実害はないのでこのままとします.

2014年11月14日金曜日

[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
>

2014年11月13日木曜日

[SICP] 問題 1.33 : filtered-accumulate

組み合せる項にフィルタ(filter)の考えを入れることで, accumulate(問題1.32)の更に一般的なものが得られる. つまり範囲から得られて, 指定した条件を満した項だけを組み合せる. 出来上ったfiltered-accumulate抽象は, accumulate と同じ引数の他, フィルタを指定する一引数の述語をとる. filtered-accumulateを手続きとして書け. filtered-accumulateを使い, 次をどう表現するかを示せ.

a. 区間a, bの素数の二乗の和(prime?述語は持っていると仮定する.)
b. nと互いに素で, nより小さい正の整数(つまりi < nでGCD(i, n)=1なる全整数i)の積

フィルターに対応したaccumulateを定義します. 考え方としては, filterが真を返すときはaccumulateと同様の処理をし, 偽を返すときはcombinerを使った計算をしないようにします. 手続きは次のようになります.

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

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

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

(define (cube x) (* 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 (filtered-accumulate-r filter combiner null-value term a next b)
  (cond ((> a b) 
         null-value)
        ((filter a)
         (combiner (term a)
                   (filtered-accumulate-r filter combiner null-value term (next a) next b)))
        (else
         (filtered-accumulate-r filter combiner null-value term (next a) next b))))

;; 反復的プロセスを生成する.
(define (filtered-accumulate-i filter combiner null-value term a next b)
  (define (iter a result)
    (cond ((> a b) 
           result)
          ((filter a) 
           (iter (next a) (combiner (term a) result)))
          (else
           (iter (next a) result))))
  (iter a null-value))

;; 区間a, bの素数の二乗の和
(define (sum-prime-square a b)
  (filtered-accumulate-r prime?
                         +
                         0
                         square
                         a
                         inc 
                         b))

;; nと互いに素で, nより小さい正の整数の積
(define (product-disjoint n)
  (define (disjoint? x) (= 1 (gcd n x)))
  (filtered-accumulate-r disjoint?
                         *
                         1
                         identity
                         1
                         inc
                         (- n 1)))

区間a, bの素数の二乗の和を求めてみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (sum-prime-square 2 10)   
87
> (+ (square 2) (square 3) (square 5) (square 7))
87
>

2から10の間の素数は, 2, 3, 5, 7ですから, 計算は合っています.

nと互いに素で, nより小さい正の整数の積を求めてみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (product-disjoint 10)
189
> (* 1 3 7 9)
189
>

10よりも小さい正の整数で10と互いに素になる数は, 1, 3, 7, 9ですから, 計算は合っています.

2014年11月12日水曜日

[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は次のような行列として表現できます.

T_pq=	[	p+q	q	]
		q	p

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

T_p'q'=	[	(p²+q²) + (2pq+q²)	2pq+q²	]
		2pq+q²	p²+q²

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

p'=	p²+q²
q'=	2pq+q²

上で求めた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
>

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

2014年11月8日土曜日

[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
>

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

2014年11月7日金曜日

[SICP] 問題 1.32 : accumulate

a. sumと(問題1.31の)productは, 一般的なアキュムレーションの関数:

(accumulate combiner null-value term a next b)

を使い, 項の集りを組み合せるaccumulateという更に一般的なものの特殊な場合であることを示せ.
accumulateは引数としてsumや productと同様, 項と範囲指定と, 先行する項のアキュムレーションと現在の項をどう組み合せるかを指定する(二引数の)combiner手続き, 項がなくなった時に使う値を指定するnull-valueをとる. accumulateを書き, sumやproductがaccumulateの単なる呼出しで定義出来ることを示せ.

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

手続きsumとproductを比較すると, 「＋」と「＊」がcombiner, 「０」と「１」がnull-valueに対応します. 手続きを定義すると次のようになります.

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

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

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

;; 再帰的プロセスを生成する.
(define (accumulate-r combiner null-value term a next b)
  (if (> a b)
      null-value
      (combiner (term a)
                (accumulate-r combiner null-value term (next a) next b))))

;; 反復的プロセスを生成する.
(define (accumulate-i 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))

;; 総和の定義
(define (sum term a next b)
  (accumulate-r + 0 term 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 (product term a next b)
  (accumulate-i * 1 term a next b))

;; 円周率を求める. 
(define (my-pi n)
  (define (term x) 
    (/ (* x (+ x 2)) (square (+ x 1))))
  (define (next x)
    (+ x 2))
  (* 4.0 (product term 2 next n)))

実行してみます.

ようこそ 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
> (my-pi 10000)
3.141749705738071
> (my-pi 100000)
3.141608361278168
>

これまでの問題と同様に計算できています.

2014年11月5日水曜日

[SICP] 問題 1.31 : productの定義

a. sum手続きは高階手続きとして書ける, 同様な数多い抽象の最も単純なものに過ぎない.51 与えられた範囲の点での関数値の積を返すproductという似た手続きを書け. productを使って, factorialを定義せよ. また式

π　　２・４・４・６・６・８・・・
−＝−−−−−−−−−−−−−−− 
４　　３・３・５・５・７・７・・・

によってπの近似値を計算するのにproductを使え.
b. 上のproductが再帰的プロセスを生成するなら, 反復的プロセスを生成するものを書け. 反復的プロセスを生成するなら, 再帰的プロセスを生成するものを書け.

与えられた範囲の関数値の積を求める手続きproductを定義します. 和を積に変えているだけなので, sumの定義を利用できます. 手続きは次のようになります.

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

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

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

;; p.18のfactorialの定義
(define (factorial-0 n)
  (if (= n 1)
      1
      (* n (factorial-0 (- n 1)))))

;; productの定義. 再帰的プロセス
(define (product-r term a next b)
  (if (> a b)
      1
      (* (term a)
         (product-r term (next a) next b))))

;; productを用いたfactorialの定義
(define (factorial n)
  (product-r identity 1 inc n))

;; 反復的プロセスを生成
(define (product-i term a next b)
  (define (iter a result)
    (if (> a b)
        result
        (iter (next a) (* (term a) result))))
  (iter a 1))

;; 円周率を求める. nが大きくなれば精度が良くなる.
(define (my-pi n)
  (define (term x) 
    (/ (* x (+ x 2)) (square (+ x 1))))
  (define (next x)
    (+ x 2))
  (* 4.0 (product-i term 2 next n)))

変更内容は, 「＋」を「＊」に, 初期値を0から1に変えているぐらいです. factorialの定義には再帰的プロセスを生成するproductを使用し, 円周率の計算には反復的プロセスを生成productを使用しています.

また, 円周率を計算する式は, 次のようにカッコをつけると見やすくなります.

π　　２・４　　　４・６　　　６・８
−＝（−−−）・（−−−）・（−−−）・・・
４　　３・３　　　５・５　　　７・７

計算してみます.

ようこそ DrRacket, バージョン 6.1 [3m].
言語: Pretty Big; memory limit: 2048 MB.
> (factorial-0 6)
720
> (factorial 6)
720
> (my-pi 10)
3.2751010413348074
> (my-pi 100)
3.1570301764551676
> (my-pi 1000)
3.1431607055322663
> (my-pi 10000)
3.1417497057380523
> (my-pi 100000)
3.1416083612781764
>

実行結果から, 円周率を計算できていそうです.

登録: 投稿 (Atom)