SICP: Exercise 1.34
「1.3.2 lambda
を使う手続きの構築」でついに無名関数 lambda
が出てきた。先に JavaScript で見慣れていたから驚きは少なかったけれど、他の言語でプログラミングを始めていたら大きく驚いたかも(とはいえ、例えば C 言語でも関数ポインタとかあるしなぁ…)。だんだん関数型プログラミングらしくなっていく。
あと、局所変数束縛の特殊形式 let
も出てきた。
(let ((var1 exp1) (var2 exp2) ... (varn expn)) body)
は
((lambda (var1 var2 ... varn) body) exp1 exp2 ... expn)
のシンタックスシュガー。「変数の値は let
の外側で計算される」ということは覚えておいた方がよさそう。
というわけで問題。
問題 1.34
手続き
(define (f g) (g 2))
を定義したとする。その時
(f square) 4 (f (lambda (z) (* z (+ z 1)))) 6
解釈系に組合せ
(f f)
を(意地悪く)評価させるとどうなるか。説明せよ。
置換モデルで少しずつ還元していく。ややこしいので「1.1.5 手続き作用の置換えモデル」の手順で一歩ずつ。
(f f) (g 2) ; g:仮引数 (f 2) ; 仮引数 g -> f (g 2) ; g:仮引数 (2 2) ; 仮引数 g -> 2 ; エラー ; 2 は手続きでないので作用させることができない
SICP: Exercise 1.33
問題 1.33
組み合せる項にフィルタ(filter)の考えを入れることで、
accumulate
(問題 1.32)の更に一般的なものが得られる。つまり範囲から得られて、指定した条件を満した項だけを組み合せる。出来上ったfiltered-accumulate
抽象は、accumulate
と同じ引数の他、フィルタを指定する一引数の述語をとる。filtered-accumulate
を手続きとして書け。filtered-accumulate
を使い、次をどう表現するかを示せ。
- 区間 \( a \), \( b \) の素数の二乗の和(
prime?
述語は持っていると仮定する。)- \( n \) と互いに素で、\( n \) より小さい正の整数(つまり \( i < n \) で \( \mathrm{GCD}(i, n) = 1 \) なる全整数 \( i \))の積
これもリストの扱いでよくありそうな処理、filter
してから reduce
。問題文にある通り、注目している番号が指定した条件を満たしているときのみ組合せを行うようにすればよい。
まずは再帰的プロセス版。
(define (filtered-accumulate combiner null-value pred term a next b) (if (> a b) null-value (if (pred a) (combiner (term a) (filtered-accumulate combiner null-value pred term (next a) next b)) (filtered-accumulate combiner null-value pred term (next a) next b))))
続いて反復的プロセス版。こちらの方が簡潔に見える。
(define (filtered-accumulate combiner null-value pred term a next b) (define (iter a result) (if (> a b) result (iter (next a) (if (pred a) (combiner (term a) result) result)))) (iter a null-value))
細かいことだけれど、引数の pred
を入れる場所としてどこが最適かで悩んだ。この順番だと自然に読めるように感じたが、どうだろうか。
a
素数の 2 乗の和。filtered-accumulate
のおかげでとても短く書ける。
(define (sum-of-squared-primes a b) (filtered-accumulate + 0 prime? square a inc b))
正しく答えが出れば prime?
はどれでも良いので、1.2.6 節からどれかを持ってくる。ここでは本文と問題 1.23 のもので試してみる。
;; 反復的プロセスの filtered-accumulate (define (filtered-accumulate combiner null-value pred term a next b) (define (iter a result) (if (> a b) result (iter (next a) (if (pred a) (combiner (term a) result) result)))) (iter a null-value)) ;; 素数判定 (define (prime? n) (= n (smallest-divisor n))) (define (smallest-divisor n) (find-divisor 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 (divides? a b) (= (remainder b a) 0)) (define (next n) (if (= n 2) 3 (+ n 2))) ;; 素数の 2 乗の和 (define (sum-of-squared-primes a b) (filtered-accumulate prime? + 0 square a inc b)) (define (square x) (* x x)) (define (inc x) (+ x 1)) ;; 2 以上 10 以下の素数(2, 3, 5, 7)の 2 乗の和を求めてみる (+ 4 9 25 49) ; => 87 (sum-of-squared-primes 2 10) ; => 87
合っていた。
b
\( n \) と互いに素で、\( n \) より小さい正の整数の積。一見ややこしい条件だけれど、互いに素(最大公約数が 1)かどうかを返す手続きを作れば分かりやすくなる。
(define (product-of-relatively-prime-numbers n) (define (relatively-prime? i) (= (gcd n i) 1)) (filtered-accumulate * 1 relatively-prime? identity 1 inc (- n 1)))
こちらもテストしてみる。
;; 10 と互いに素で、それより小さい正の整数の積を求めてみる (* 1 3 7 9) ; => 189 (product-of-relatively-prime-numbers 10) ; => 189
大丈夫だった。
抽象化のおかげで、いずれの場合も必要な部分だけ簡潔に書くことができた。無駄がなくてすっきり。
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) $$
となる。