2015-05-17 の日記
このところ学校やボランティアで作業に参加しているクリエイターズネットワークで、これまで触れてこなかった技術に触れる機会が多くて、その度にいろいろなことができるようになっている。
学校では ET ロボコン 2015に向けての勉強が始まった。研究室の先生によるモデリングの特別講義を受けたり、その後で Git 勉強会を開いてみんなでコミットできるようになったり。今週末にはメンバーのひとりが名古屋に行って研修を受けるのだけれど、そのときの使用言語が学校では習っていない C++ で、慌てて一緒に勉強したりとか。とにかく授業だけでは不十分なのが実感できるのだけれど、メンバー全員が結構興味を持って取り組めているので、非常に良い雰囲気。できるだけ続けていきたい。
クリエイターズネットワークの方では最近 HTTP サーバーの調子が悪いのでいろいろ調整しているのだけれど、負荷削減やアクセスログ解析等のために Nginx を使って Web アプリを動かすための作業を行う機会が増えた。WordPress やアクセスログ解析の AWStats といったアプリごとにつまづく点があって苦労したけれど、変更するファイルの書き方を少しずつ理解できてきたので、今後同様の作業をするときには多少楽に進められるようになりそう。なんとか安定動作に持っていきたい。
日々が充実していて時間が足りないのが一番の悩みだ。
SICP: Exercise 1.32
問題 1.32
sumと（問題 1.31 の）productは、一般的なアキュムレーションの関数：(accumulate combiner null-value term a next b)を使い、項の集りを組み合せる
accumulateという更に一般的なものの特殊な場合であることを示せ。accumulateは引数としてsumやproductと同様、項と範囲指定と、先行する項のアキュムレーションと現在の項をどう組み合せるかを指定する（二引数の）combiner手続き、項がなくなった時に使う値を指定するnull-valueをとる。accumulateを書き、sumやproductがaccumulateの単なる呼出しで定義出来ることを示せ。上の
accumulateが再帰的プロセスを生成するなら、反復的プロセスを生成するものを書け。反復的プロセスを生成するなら、再帰的プロセスを生成するものを書け。
脚注 51 で触れられているように、後に並び（リスト）が出たとき、アキュミュレーションは fold とか reduce としてよく使われる。
まずは再帰的プロセスで書いてみる。
(define (accumulate combiner null-value term a next b) (if (> a b) null-value (combiner (term a) (accumulate combiner null-value term (next a) next b))))
次に反復的プロセスで書く。問題 1.30 と同様。
(define (accumulate combiner null-value term a next b) (define (iter a result) (if (> a b) result (iter (next a) (combiner result (term a))))) (iter a null-value))
sum や product を accumulate を使って書くとこうなる。
(define (sum term a next b) (accumulate + 0 term a next b)) (define (product term a next b) (accumulate * 1 term a next b))
とても簡潔。
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 \) が十分に大きければ予想の通りになると考えられる。
