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