SICP: Exercise 1.31
問題 1.31
sum手続きは高階手続きとして書ける、同様な数多い抽象の最も単純なものに過ぎない。与えられた範囲の点での関数値の積を返すproductという似た手続きを書け。productを使って、factorialを定義せよ。また式 $$ \frac{\pi}{4} = \frac{2 \cdot 4 \cdot 4 \cdot 6 \cdot 6 \cdot 8 \cdots}{3 \cdot 3 \cdot 5 \cdot 5 \cdot 7 \cdot 7 \cdots} $$ によって \( \pi \) の近似値を計算するのにproductを使え。上の
productが再帰的プロセスを生成するなら、反復的プロセスを生成するものを書け。反復的プロセスを生成するなら、再帰的プロセスを生成するものを書け。
a
和と同様に積 $$ \prod_{n = a}^b\,f(n) = f(a) \times \cdots \times f(b) $$ を表す手続きを書けばよい。
(define (product term a next b) (if (> a b) 1 (* (term a) (product term (next a) next b))))
階乗
階乗は sum-integers と同様になる。
(define (factorial n) (product identity 1 inc n)) (factorial 0) ; 1 (factorial 1) ; 1 (factorial 5) ; 120 (factorial 10) ; 3628800
π の近似値
\( \pi \) の近似値は、初見時は一般項の表現に困ったが、 $$ \frac{2 \cdot 4}{3 \cdot 3} \cdot \frac{4 \cdot 6}{5 \cdot 5} \cdot \frac{6 \cdot 8}{7 \cdot 7} \cdots $$ と見たら分かりやすくなった。左上の数を \( k \) とおくと $$ \frac{k(k + 2)}{(k + 1)^2} $$ と表せる。このように表現したときは、\( k \) に 2 を加えて次の項を作る。
(define (pi-prod n) (define (pi-term x) (/ (* x (+ x 2.0)) (square (+ x 1)))) (define (pi-next x) (+ x 2)) (product pi-term 2 pi-next n))
実際に近似値を計算すると
(* 4 (pi-prod 10000000)) ; 少し時間がかかる ; 3.1415928106708866
3.141592 まで合った。
b
a の product は再帰的プロセスなので反復的プロセスで書く。問題 1.30 と同様に書けばよい。
(define (product term a next b) (define (iter a result) (if (> a b) result (iter (next a) (* (term a) result)))) (iter a 1))
SICP: Exercise 1.30
問題 1.30
上の
sumの手続きは線形再帰を生成する。総和が反復的に実行出来るように手続きを書き直せる。次の定義の欠けているところを補ってこれを示せ：(define (sum term a next b) (define (iter a result) (if ⟨??⟩ ⟨??⟩ (iter ⟨??⟩ ⟨??⟩))) (iter ⟨??⟩ ⟨??⟩))
元の sum は再帰的プロセス。
(define (sum term a next b) (if (> a b) 0 (+ (term a) (sum term (next a) next b))))
例えば 1〜3 の和を求めるときは、置換モデルだと以下のようになる。
(sum identity 1 inc 3) (+ 1 (sum identity 2 inc 3)) (+ 1 (+ 2 (sum identity 3 inc 3))) (+ 1 (+ 2 (+ 3 (sum identity 4 inc 3)))) (+ 1 (+ 2 (+ 3 0))) (+ 1 (+ 2 3)) (+ 1 5) 6
反復的プロセスで書くと
(define (sum term a next b) (define (iter a result) (if (> a b) result (iter (next a) (+ (term a) result)))) (iter a 0))
上の例は以下の結果になる。
(sum identity 1 inc 3) (iter 1 0) (iter 2 1) (iter 3 3) (iter 4 6) 6
SICP: Exercise 1.29
1.3 節「高階手続きによる抽象」から、ついに lambda が登場する。SICP を始める前には主に JavaScript を使っていたのだけれど、『JavaScript: The Good Parts』で匿名関数や高階関数の強力さが力説されているのを見てきたので、この辺りはかなり楽しく進めていた記憶がある。
というわけで、問題 1.29。
問題 1.29
Simpson の公式は上に示したのより更に正確な数値積分法である。Simpson の公式を使えば、\( a \) から \( b \) までの関数 \( f \) の積分は、偶整数 \( n \) に対し、\( h = \frac{b - a}{n} \) また \( y_k = f(a + kh) \) として $$ \frac{h}{3} [ y_0 + 4y_1 + 2y_2 + 4y_3 + 2y_4 + \cdots + 2y_{n - 2} + 4y_{n - 1} + y_n ] $$ で近似出来る。（\( n \) の増加は近似の精度を増加する。）引数として \( f \), \( a \), \( b \) と \( n \) をとり、Simpson 公式を使って計算した積分値を返す手続きを定義せよ。その手続きを使って（\( n = 100 \) と \( n = 1000 \) で）0 から 1 まで
cubeを積分し、またその結果を上のintegral手続きの結果と比較せよ。
高階手続きの導入に伴って、これまでのものより抽象的な手続きが出てきた。
sum は総和
$$ \sum_{n = a}^b\,f(n) = f(a) + \cdots + f(b) $$
を表す手続き。
;; term および next は手続き (define (sum term a next b) (if (> a b) 0 (+ (term a) (sum term (next a) next b))))
inc は引数を 1 増やす手続き。
(define (inc n) (+ n 1))
恒等手続き。
(define (identity x) x)
定積分は細かく区間を分けて長方形で近似する。\( \Delta x \) を小さな実数として $$ \int_a^b f(x)\,dx \approx \left\{ f \left( a + \frac{\Delta x}{2} \right) + f \left( a + \Delta x + \frac{\Delta x}{2} \right) + f \left( a + 2\Delta x + \frac{\Delta x}{2} \right) + \cdots \right\} \Delta x $$
手続きとして書くと
(define (integral f a b dx) (define (add-dx x) (+ x dx)) (* (sum f (+ a (/ dx 2.0)) add-dx b) dx))
解答
\( n \) は偶数のみを想定している。奇数のときにエラーとする手もありそうだけれど、係数が 4, 2, 4, … となるだけなので、チェックは省いた。
(define (integral-simpson f a b n) (let ((h (/ (- b a) n))) (define (y k) (f (+ a (* k h)))) (define (term k) (* (cond ((or (= k 0) (= k n)) 1.0) ((even? k) 2.0) (else 4.0)) (y k))) (* (/ h 3) (sum term 0 inc n))))
integral との比較
0 から 1 まで cube を積分した結果で、integral と integral-simpson を比較した。正確な値は 0.25。
;; 積分値、誤差を計算する ;; n = 100 (integral cube 0 1 0.01) ; 0.24998750000000042 (- (integral cube 0 1 0.01) 0.25) ; -1.249999999958229e-5 (integral-simpson cube 0 1 100) ; 0.24999999999999992 (- (integral-simpson cube 0 1 100) 0.25) ; -8.326672684688674e-17 ;; n = 1000 (integral cube 0 1 0.001) ; 0.249999875000001 (- (integral cube 0 1 0.001) 0.25) ; -1.2499999899051595e-7 (integral-simpson cube 0 1 1000) ; 0.2500000000000002 (- (integral-simpson cube 0 1 1000) 0.25) ; 2.220446049250313e-16
長方形近似と比べると、Simpson の公式を使えばかなり精度が上がることが分かる。面白いのが、\( n = 100 \) のときの方が \( n = 1000 \) のときより真の値に近いこと。丸め誤差が蓄積したのだろうか。
SICP: Exercise 1.28
問題 1.28
だまされない Fermat テストの変形の一つは Miller-Rabin テスト（Miller-Rabin test）（Miller 1976; Rabin 1980）という。これは、\( n \) が素数であり、\( a \) を \( n \) より小さい任意の正の整数とすると、\( a \) の \( (n - 1) \) 乗は \( n \) を法として 1 と合同であるという、Fermat の小定理のもう一つの方から始める。Miller-Rabin テストで数 \( n \) の素数性をテストするには、\( a < n \) をランダムにとり、
expmod手続きを使って \( n \) を法とした \( a \) の \( (n - 1) \) 乗を作る。しかし、expmodで二乗を計算しながら「\( n \) を法とした 1 の自明でない平方根」、つまり 1 でも \( n - 1 \) でもないがその二乗が \( n \) を法として 1 になる数がなかったか調べる。そういう、\( n \) を法とした 1 の自明でない平方根があれば \( n \) が素数でないことは証明出来る。また \( n \) が素数でない奇数なら、\( a < n \) の少なくても半分は \( a^n-1 \) を計算する時、\( n \) を法とした 1 の自明でない平方根が出てくることも証明出来る。（これが Miller-Rabin テストがだまされない理由である。）\( n \) を法とした 1 の自明でない平方根を見つけたらシグナルを出すようにexpmod手続きを修正せよ。これを使いfermat-testに似た Miller-Rabin テストの手続きを実装せよ。分っている素数、非素数をテストし、手続きを検査せよ。ヒント：expmodにシグナルを出させる方法の一つは 0 を返させることである。
「素数性のテスト」の最後の問題は、改良された素数性テストを実装するもの。問題文は長いが、丁寧に説明されているので、それに従えば素直に実装することができる。
「1 の自明でない平方根」の話は、Wikipedia に載っていた証明が分かりやすかった。
まず、有限体 \( \mathbb{Z}_p \) の単位元の平方根についての補題を考える。ここで、\( p \) は奇素数である。\( p \) を法とした剰余\( \pmod{p} \) において、1 と -1 の二乗は常に 1 となる。これらを 1 の「自明な」平方根と呼ぶ。\( p \) は素数なので 1 の「自明でない」平方根は存在しない。これを示すため、1 mod \( p \) の非自明な平方根を \( x \) とする。このとき、以下が成り立つ： $$ x^2 \equiv 1\pmod{p} $$ $$ x^{2} - 1 \equiv 0\pmod{p} $$ $$ (x - 1)(x + 1) \equiv 0\pmod{p} $$ \( p \) は素数なので、これは \( x - 1 \) または \( x + 1 \) が \( p \) で割り切れることを示す。一方、\( x \) は 1 でも -1 でもないので\( \pmod{p} \)、\( x - 1 \) も \( x + 1 \) も \( p \) で割り切れない。すなわち、二乗して \( 1 \pmod{p}\) となるのは、1 および -1 だけということになる。
解答
(define (miller-rabin-test n) (define (expmod base exp m) ;; 1 の非自明な平方根かどうかを調べる (define (check-nontrivial-sqrt-1 x) (cond ((= x 1) 1) ; 1 の自明な平方根 ((= x (- m 1)) 1) ; 1 の自明な平方根 (else (let ((r (remainder (square x) m))) (if (= r 1) 0 ; 1 の非自明な平方根であることを示すシグナル r))))) (cond ((= exp 0) 1) ((even? exp) (check-nontrivial-sqrt-1 (expmod base (/ exp 2) m))) (else (remainder (* base (expmod base (- exp 1) m)) m)))) (define (try-it a) (= (expmod a (- n 1) n) 1)) (try-it (+ 1 (random (- n 1))))) (define (fast-prime? n times) (cond ((= times 0) #t) ((miller-rabin-test n) (fast-prime? n (- times 1))) (else #f))) (define (prime? n) ;; Miller-Rabin テストを 20 回行うので誤りの確率は ;; (1/4)^20 ≈ 10^(-12) (fast-prime? n 20))
テスト
素数、非素数でテストせよとのことなので、とりあえず 100 未満の数から素数を選んでみる。
(use srfi-1) (define numbers (iota 99 2)) ; 2〜100 までの整数のリスト (filter prime? numbers) ; (2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97)
何回か試したが、結果は 100 未満の素数のリストになり、大丈夫そう。
問題 1.27 で扱った Carmichael 数はどうだろうか。
(define carmichael-numbers '(561 1105 1729 2465 2821 6601)) (map prime? carmichael-numbers) ; (#f #f #f #f #f #f)
Carmichael 数にもだまされていない！ したがって、Fermat テストを改良することができた。
SICP: Exercise 1.27
問題 1.27
脚注 47 の Carmichael 数が本当に Fermat テストをだますことを示せ。それには整数 \( n \) をとり、\( a < n \) なるすべての \( a \) で、\( a^n \) が \( n \) を法として \( a \) の合同になるかどうか見る手続きを書き、Carmichael 数でその手続きを使ってみる。
脚注 47 で挙げられている Carmichael 数は 561, 1105, 1729, 2465, 2821 および 6601。問題文に従ってテストする。
Fermat テストをだます数（Fermat の小定理より素数もそのような数になる）かどうかを調べる手続きは
(define (fool-fermat-test? n) (define (iter a) (cond ((= a 0) #t) ((= (expmod a n n) a) (iter (- a 1))) (else #f))) (iter (- n 1)))
それではまとめて調べてしまおう。
(define carmichael-numbers '(561 1105 1729 2465 2821 6601)) (map fool-fermat-test? carmichael-numbers) ; (#t #t #t #t #t #t)
確かに脚注 47 の Carmichael 数はすべて Fermat テストをだますようだ。
SICP: Exercise 1.26
問題 1.26
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 がいった。「手続きをこのように書き替えたので、\( \Theta(\log n) \) のプロセスを \( \Theta(n) \) のプロセスにしてしまった。」説明せよ。
expmod の exp に注目し、どのように手続きが呼び出されるか考える。
square を使う
(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))))
の場合、以下の図のように呼び出される。段階数は \( \Theta(\log n) \)。
一方、Louis の
(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))))
の場合、(*) の部分に到達するたびに expmod が 2 回呼び出されるので、以下の図のようになる。
図で比較すると、Louis の書き方のまずさがよく分かる。(*) に到達する回数は \( \Theta(\log_2 n) \) であり、そこで経路が 2 つに分かれるので、全体の段階数は
$$ \Theta\left( 2^{\log_2 n} \right) = \Theta(n) $$
となる。
SICP: Exercise 1.25
問題 1.25
Alyssa P. Hacker は
expmodを書くのに多くの余分なことをやったと不満で、結局べき乗の計算法は知っているから(define (expmod base exp m) (remainder (fast-expt base exp) m))と単純に書ける筈だといった。これは正しいか。これも高速素数テストと同じに使えるか、説明せよ。
正しくない。
これが正しくないことは、expmod を大きな exp に対して行うことで分かる。試しに \( 5^{10000000} \mod 10000000 \) を計算してみよう。まず問題 1.16 から繰り返しプロセスの fast-expt を持ってきて、Alyssa が示した expmod を実行する。
(define (fast-expt b n) (define (fast-expt-iter a b n) (cond ((= n 0) a) ((even? n) (fast-expt-iter a (square b) (/ n 2))) (else (fast-expt-iter (* a b) b (- n 1))))) (fast-expt-iter 1 b n)) (define (expmod base exp m) (remainder (fast-expt base exp) m)) (expmod 5 10000000 10000000) ; なかなか終わらない…
一方、本文で示された 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)))) (expmod 5 10000000 10000000) ; 2890625（一瞬で出る）
となり、計算時間に大きな差がある。
fast-expt を使うと最終的にとても大きい数を \( m \) で割った余りを求めることになるので、計算が遅くなったり桁あふれが起こる可能性がある。各段階で remainder を適用すれば、脚注 46 で
指数 \( e \) が 1 より大きい場合の簡約化では、任意の整数 \( x \)、\( y \) と \( m \) について、\( x \times y \mod m \) の剰余は、\( x \mod m \) と \( y \mod m \) の剰余を別々に計算し、その結果を掛け合せ、更に積の \( \mod m \) の剰余をとればよいという事実に基づいている。\( e \) が偶数の場合は、\( b^{e/2} \mod m \) の剰余を計算し、二乗し、その \( \mod m \) の剰余をとる。これは \( m \) より遥かに大きい数値を扱わずに計算が実行出来、有用である。
と指摘されているように、\( m^2 \) 以上の数を扱わずに計算できるので、上記の問題が生じない。
SICP: Exercise 1.24
問題 1.24
問題 1.22 の
timed-prime-test手続きをfast-prime?（Fermat 法参照）を使うように変え、その問題で見つけた 12 個の素数をテストせよ。Fermat テストは \( \Theta(\log n) \) の増加なので、1,000,000 近くの素数をテストする時間を 1,000 近くの素数をテストする時間と比べ、どの位と予想するか。データはその予想を支持するか。違いがあれば説明出来るか。
下準備
Gauche で random を使えるようにするために、srfi-27 モジュールを利用する。以下、本文の手続きも含めてまとめて定義する。
(use srfi-27) (define (random n) (random-integer n)) (random-source-randomize! default-random-source) (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 (fermat-test n) (define (try-it a) (= (expmod a n n) a)) (try-it (+ 1 (random (- n 1))))) (define (fast-prime? n times) (cond ((= times 0) #t) ((fermat-test n) (fast-prime? n (- times 1))) (else #t))) (load "./ex-1-22.scm")
解答
とりあえず 10 回 Fermat テストを行うように書いた。
(define (start-prime-test n start-time) (if (fast-prime? n 10) (report-prime n (- (runtime) start-time))))
1,000,000 近くの素数をテストする時間と 1,000 近くの素数をテストする時間の比は、段階数が \( \Theta(\log n) \) なので
$$ \frac{\log 1000000}{\log 1000} = \frac{\log (10^6)}{\log (10^3)} = \frac{6 \log 10}{3 \log 10} = 2 $$
と予想した。
結果
問題 1.23 と同様に 12 個の素数について、素数テストに要する時間を測った。
結果を以下に示す。
| 素数：\( p \) | 問題 1.23：\( t_\mathrm{A}/\mathrm{\mu s} \) | 問題 1.24：\( t_\mathrm{B}/\mathrm{\mu s} \) | \( t_\mathrm{A} / t_\mathrm{B} \) |
|---|---|---|---|
| 1,009 | 7 | 47 | 0.1 |
| 1,013 | 4 | 39 | 0.1 |
| 1,019 | 5 | 42 | 0.1 |
| 10,007 | 13 | 48 | 0.3 |
| 10,009 | 13 | 46 | 0.3 |
| 10,037 | 13 | 47 | 0.3 |
| 100,003 | 38 | 54 | 0.7 |
| 100,019 | 39 | 55 | 0.7 |
| 100,043 | 39 | 55 | 0.7 |
| 1,000,003 | 122 | 60 | 2.0 |
| 1,000,033 | 122 | 62 | 2.0 |
| 1,000,037 | 122 | 62 | 2.0 |
毎回 10 回ずつ Fermat テストを行っているためか \( p \) が小さいうちは遅いけれど、\( p \approx 1,000,000 \) ではかなり速くなっている。
追加でもう一桁大きいときの時間を計測。
| 素数：\( p \) | 問題 1.23：\( t_\mathrm{A}/\mathrm{\mu s} \) | 問題 1.24：\( t_\mathrm{B}/\mathrm{\mu s} \) | \( t_\mathrm{A} / t_\mathrm{B} \) |
|---|---|---|---|
| 10,000,019 | 395 | 74 | 5.3 |
| 10,000,079 | 395 | 75 | 5.3 |
| 10,000,103 | 394 | 74 | 5.3 |
さらに差が広がった。桁数が大きくなるにつれて \( \Theta(\sqrt{n}) \) と \( \Theta(\log n) \) の差は大きくなっていくのだろう。
1,000,000 近くの素数をテストする時間と 1,000 近くの素数をテストする時間の比は、平均をとって計算してみると
$$ \frac{61.3}{42.7} \approx 1.44 $$
おや、思ったより大きくない。改めて fermat-test を見てみる。
(define (fermat-test n) (define (try-it a) (= (expmod a n n) a)) (try-it (+ 1 (random (- n 1)))))
expmod は \( \Theta(\log n) \) だがその寄与が大きくないということは、random で時間がかかっているということだろう。
そこで、p を大きくして expmod に要する時間を長くしてみる。
(search-for-primes 100000000 100000100) ; 100000007 *** 92 ; 100000037 *** 86 ; 100000039 *** 88 ; 平均 88.7 (search-for-primes 100000000000 100000000100) ; 100000000003 *** 180 ; 100000000019 *** 176 ; 100000000057 *** 173 ; 平均 176.3
$$ \frac{176.3}{88.7} \approx 1.99 $$
このようにすると 1,000 倍したときとの時間の比が 2 に近づいた。したがって、\( p \) が十分に大きければ予想の通りになると考えられる。
SICP: Exercise 1.23
問題 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 と違う事情を説明せよ。
解答
問題文の通りに next を定義し、(next test-divisor) を使うように修正する。
(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)))))
結果
問題 1.22 で調べた 12 個の素数について、改良した素数テストに要する時間を測った。
(for-each timed-prime-test '(1009 1013 1019 10007 10009 10037 100003 100019 100043 1000003 1000033 1000037))
結果を以下に示す。
| 素数：\( p \) | 問題 1.22：\( t_\mathrm{A}/\mathrm{\mu s} \) | 問題 1.23：\( t_\mathrm{B}/\mathrm{\mu s} \) | \( t_\mathrm{A} / t_\mathrm{B} \) |
|---|---|---|---|
| 1,009 | 7 | 7 | 1.0 |
| 1,013 | 7 | 4 | 1.8 |
| 1,019 | 7 | 5 | 1.4 |
| 10,007 | 20 | 13 | 1.5 |
| 10,009 | 20 | 13 | 1.5 |
| 10,037 | 21 | 13 | 1.6 |
| 100,003 | 64 | 38 | 1.7 |
| 100,019 | 64 | 39 | 1.6 |
| 100,043 | 65 | 39 | 1.7 |
| 1,000,003 | 199 | 122 | 1.6 |
| 1,000,033 | 209 | 122 | 1.7 |
| 1,000,037 | 198 | 122 | 1.6 |
\( p \) が大きくなると \( t_\mathrm{A} / t_\mathrm{B} \) は 1.6〜1.7 で安定した。2 には届いていないが、ほぼ期待通りといえるだろう。遅くなった原因は、next 内の if や、next の呼び出し自体に時間を要するためと考えられる。
SICP: Exercise 1.22
2 年ぶりに SICP の記事を書くことになった。1 年半前の鬱々としていた日々がなければもっと早く進んだんだけどなぁ…
問題 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 より大きい、小さい方から三つの素数を見つけよ。素数をテストする時間に注意せよ。テストのアルゴリズムは \( \Theta(\sqrt{n}) \) の増加の程度だから、10,000 前後の素数のテストは 1,000 前後の素数のテストの \( \sqrt{10} \) 倍かかると考えよ。時間のデータはこれを支持するか。100,000、1,000,000 のデータは \( \sqrt{n} \) 予想のとおりだろうか。結果はプログラムが計算に必要なステップ数に比例した時間で走るという考えに合っているか。
Gauche で runtime を使えるようにする
「殆んどの Lisp の実装には基本手続き runtime があり」とあるが、Gauche にはなかった。「実行時間の計測 - n-oohiraのSICP日記」より以下を拝借。
(define (runtime) (use srfi-11) (let-values (((a b) (sys-gettimeofday))) (+ (* a 1000000) b)))
解答
素数は 2 以上なので、一応その範囲だけ調べるようにした。また、3 以上の素数は必ず奇数なので、奇数を調べた後は一つ置きに調べるようにした。
(define (search-for-primes start end) (cond ((> start end) 'finished) ((< start 2) (search-for-primes 2 end)) (else (timed-prime-test start) (search-for-primes (+ start (if (even? start) 1 2)) end))))
結果
search-for-primes を評価すると以下のように出力される。
(search-for-primes 20 30) ; 20 ; 21 ; 23 *** 1 ; 25 ; 27 ; 29 *** 2
長いので、1,000 以上のものは最初の素数 3 つのみ書くことにする。
(search-for-primes 1000 1050) ; 1009 *** 7 ; 1013 *** 7 ; 1019 *** 7 (search-for-primes 10000 10050) ; 10007 *** 20 ; 10009 *** 20 ; 10037 *** 21 (search-for-primes 100000 100050) ; 100003 *** 64 ; 100019 *** 64 ; 100043 *** 65 (search-for-primes 1000000 1000050) ; 1000003 *** 199 ; 1000033 *** 209 ; 1000037 *** 198
素数テストに要した時間 [μsec] は、オーダーが同じならばほぼ等しいので、上記の結果の平均をとってみると
- \( n \approx 1000 \)：7.0
- \( n \approx 10000 \)：20.7
- \( n \approx 100000 \)：64.7
- \( n \approx 1000000 \)：202.0
となった。一つ上のオーダーでの値との比を調べると
- \( \frac{20.7}{7.0} \approx 3.0 \)
- \( \frac{64.7}{20.7} \approx 3.13 \)
- \( \frac{202.7}{64.7} \approx 3.13 \)
となり、すべて \( \sqrt{10} \approx 3.16 \) に近い値だった。したがって \( \sqrt{n} \) 予想は成立しており、プログラムが計算に必要なステップ数に比例した時間で走るという考えに合っているといえる。
今回の範囲では \( \sqrt{n} \) 予想はずっと成立していたけれど、今より速度が出ない昔のコンピュータを使っていたらどうなるのだろうか。問題の最後の文を見ると、ちょっと引っかかる。
