SICP: Exercise 1.36
問題 1.36
問題 1.22 で示した基本の
newlineとdisplayを使い、生成する近似値を順に印字するようfixed-pointを修正せよ。次に \( x \mapsto \log(1000)/\log(x) \) の不動点を探索することで、\( x^x = 1000 \) の解を見つけよ。（自然対数を計算する Scheme の基本log手続きを使う。）平均緩和を使った時と使わない時のステップ数を比べよ。（fixed-pointの予測値を 1 にして始めてはいけない。\( \log(1) = 0 \) による除算を惹き起すからだ。）
fixed-point の方は、printf デバッグのように display を追加するだけ。
(define tolerance 0.00001) (define (fixed-point f guess) (define (close-enough? v1 v2) (< (abs (- v1 v2)) tolerance)) (define (try guess) (let ((next (f guess))) (display "guess: ") (display guess) (display ", next: ") (display next) (newline) (if (close-enough? guess next) next (try next)))) (try guess))
平均緩和を行わない場合は
(define (f x) (/ (log 1000) (log x))) (define (solve-x^x=1000) (fixed-point f 2.0)) (solve-x^x=1000) ; guess: 2.0, next: 9.965784284662087 ; guess: 9.965784284662087, next: 3.004472209841214 ; guess: 3.004472209841214, next: 6.279195757507157 ; guess: 6.279195757507157, next: 3.759850702401539 ; guess: 3.759850702401539, next: 5.215843784925895 ; guess: 5.215843784925895, next: 4.182207192401397 ; guess: 4.182207192401397, next: 4.8277650983445906 ; guess: 4.8277650983445906, next: 4.387593384662677 ; guess: 4.387593384662677, next: 4.671250085763899 ; guess: 4.671250085763899, next: 4.481403616895052 ; guess: 4.481403616895052, next: 4.6053657460929 ; guess: 4.6053657460929, next: 4.5230849678718865 ; guess: 4.5230849678718865, next: 4.577114682047341 ; guess: 4.577114682047341, next: 4.541382480151454 ; guess: 4.541382480151454, next: 4.564903245230833 ; guess: 4.564903245230833, next: 4.549372679303342 ; guess: 4.549372679303342, next: 4.559606491913287 ; guess: 4.559606491913287, next: 4.552853875788271 ; guess: 4.552853875788271, next: 4.557305529748263 ; guess: 4.557305529748263, next: 4.554369064436181 ; guess: 4.554369064436181, next: 4.556305311532999 ; guess: 4.556305311532999, next: 4.555028263573554 ; guess: 4.555028263573554, next: 4.555870396702851 ; guess: 4.555870396702851, next: 4.555315001192079 ; guess: 4.555315001192079, next: 4.5556812635433275 ; guess: 4.5556812635433275, next: 4.555439715736846 ; guess: 4.555439715736846, next: 4.555599009998291 ; guess: 4.555599009998291, next: 4.555493957531389 ; guess: 4.555493957531389, next: 4.555563237292884 ; guess: 4.555563237292884, next: 4.555517548417651 ; guess: 4.555517548417651, next: 4.555547679306398 ; guess: 4.555547679306398, next: 4.555527808516254 ; guess: 4.555527808516254, next: 4.555540912917957 ; guess: 4.555540912917957, next: 4.555532270803653 4.555532270803653
となった。34 ステップ。
\( x \mapsto \log(1000)/\log(x) \) の不動点探索で平均緩和を行う場合は、本文にあるとおり両辺に \( x \) を足して 2 で割る。 $$ \frac{x + x}{2} = \frac{x + \log(1000)/\log(x)}{2} \quad \therefore\ x = \frac{x + \log(1000)/\log(x)}{2} $$
コードを書いて実行すると
(define (solve-x^x=1000-with-average-damping) (fixed-point (lambda (x) (average x (f x))) 2.0)) (solve-x^x=1000-with-average-damping) ; guess: 2.0, next: 5.9828921423310435 ; guess: 5.9828921423310435, next: 4.922168721308343 ; guess: 4.922168721308343, next: 4.628224318195455 ; guess: 4.628224318195455, next: 4.568346513136242 ; guess: 4.568346513136242, next: 4.5577305909237005 ; guess: 4.5577305909237005, next: 4.555909809045131 ; guess: 4.555909809045131, next: 4.555599411610624 ; guess: 4.555599411610624, next: 4.5555465521473675 ; guess: 4.5555465521473675, next: 4.555537551999825 4.555537551999825
となった。9 ステップと大幅にステップ数が縮まった。
SICP: Exercise 1.35
問題 1.35
（1.2.2 節の）黄金比 \( \phi \) が変換 \( x \mapsto 1 + \frac{1}{x} \) の不動点であることを示し、この事実を使い
fixed-point手続きにより \( \phi \) を計算せよ。
\( f(x) = x \) となる \( x \) を関数 \( f \) の不動点という。グラフを考えると、不動点において $$ \begin{cases} y = f(x) \\ y = x \end{cases} $$ と書けるから、これはグラフ \( y = f(x) \) と直線 \( y = x \) の交点である。また、 $$ f(x),\,f(f(x)),\,f(f(f(x))),\,\cdots $$ と \( f \) を繰り返し作用させていけば不動点 \( x_0 \) に収束していくとき、\( x_0 \) を吸引的不動点というらしい。
方程式の解の探索を不動点探索によって行える場合がある。方程式を変形して \( x = f(x) \) と書くことができ、\( f(x) \) が吸引的不動点を持つならば、不動点探索によって解 \( x \) を求めることができる。
さて、これを黄金比 \( \phi \) を求めるのに利用してみる。1.2.2 節より \( \phi^2 = \phi + 1 \) だったから、両辺を \( \phi \) で割って $$ \phi = 1 + \frac{1}{\phi} $$ したがって、\( f(x) = 1 + \frac{1}{x} \)、つまり \( x \mapsto 1 + \frac{1}{x} \) とおけば \( \phi \) はその不動点である。
実際に計算してみると
;; 不動点探索で求めた黄金比 (define golden-ratio-as-fixed-point (fixed-point (lambda (x) (+ 1 (/ 1 x))) 1.0)) 1.6180327868852458 ;; 正確な黄金比 (define accurate-golden-ratio (/ (+ 1 (sqrt 5)) 2.0)) 1.618033988749895 ;; 誤差 (- golden-ratio-as-fixed-point accurate-golden-ratio) -1.2018646491362972e-6
そこそこの精度で求められている。
黄金比に収束することの証明
詳しい話が Amaca, Edgar Gilbuena, "On rational functions with Golden Ratio as fixed point" (2008). Thesis. Rochester Institute of Technology. に載っていた。問題 1.13 とこの論文を参考にして、上記の関数を繰り返し作用させれば黄金比に収束することが理解できた。
本文の平方根の計算と同様に、上記の \( f \) を \( x \) に \( n \) 回作用させた結果を \( y_n \) と書くことにする。まず、\( n \) が小さいときの \( y_n \) を簡単な分数で表してみる。 \begin{align*} y_1 &= 1 + \frac{1}{x} = \frac{x + 1}{x} \\ y_2 &= 1 + \frac{1}{\displaystyle 1 + \frac{1}{x}} = 1 + \frac{1}{\displaystyle \frac{x + 1}{x}} = 1 + \frac{x}{x + 1} = \frac{2x + 1}{x + 1} \\ y_3 &= 1 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{x}}} = 1 + \frac{1}{\displaystyle \frac{2x + 1}{x + 1}} = 1 + \frac{x + 1}{2x + 1} = \frac{3x + 2}{2x + 1} \\ y_4 &= 1 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{x}}}} = 1 + \frac{1}{\displaystyle \frac{3x + 2}{2x + 1}} = 1 + \frac{2x + 1}{3x + 2} = \frac{5x + 3}{3x + 2} \\ y_5 &= 1 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{\displaystyle 1 + \frac{1}{x}}}}} = 1 + \frac{1}{\displaystyle \frac{5x + 3}{3x + 2}} = 1 + \frac{3x + 2}{5x + 3} = \frac{8x + 5}{5x + 3} \end{align*}
この辺りまで出してみれば十分だろう。フィボナッチ数を \( F_n \) として、すべての自然数 \( n \) に対して以下が成り立つことが予想できる。 $$ y_n = \frac{F_{n + 1} x + F_n}{F_n x + F_{n - 1}} $$ ただし、 $$ \begin{cases} F_0 = 0 \\ F_1 = 1 \\ F_n = F_{n - 1} + F_{n - 2} \quad (n \geq 2) \end{cases} $$ である。
予想が正しいことを数学的帰納法で示す。
\( n = 1 \) のときは成立している。\( n = k \) のときに成立する、すなわち $$ y_k = \frac{F_{k + 1} x + F_k}{F_k x + F_{k - 1}} $$ であると仮定する。このとき $$ \begin{align*} y_{k + 1} &= 1 + \frac{1}{y_k} = 1 + \frac{F_k x + F_{k - 1}}{F_{k + 1} x + F_k} = \frac{(F_{k + 1} + F_k)x + (F_k + F_{k - 1})}{F_{k + 1} x + F_k} \\ &= \frac{F_{k + 2} x + F_{k + 1}}{F_{k + 1} x + F_k} \end{align*} $$ であるから、仮定の下では \( n = k + 1 \) でも成立する。以上から、すべての自然数 \( n \) に対して $$ y_n = \frac{F_{n + 1} x + F_n}{F_n x + F_{n - 1}} $$ が成立することが示された。
次に $$ \lim_{n \to \infty} y_n = \phi $$ を示す。
問題 1.13 より $$ F_n = \frac{\phi^n - \psi^n}{\sqrt{5}} \quad \left( \phi = \frac{1 + \sqrt{5}}{2}, \, \psi = \frac{1 - \sqrt{5}}{2} \right) $$ であった。
まず、以下の極限を求める。 $$ \lim_{n \to \infty} \frac{F_{n + 1}}{F_n} = \lim_{n \to \infty} \frac{\displaystyle \frac{\phi^{n + 1} - \psi^{n + 1}}{\sqrt{5}}}{\displaystyle \frac{\phi^n - \psi^n}{\sqrt{5}}} = \lim_{n \to \infty} \frac{\phi^{n + 1} - \psi^{n + 1}}{\phi^n - \psi^n} $$ \( \psi \) について $$ -1 = \frac{1 - \sqrt{9}}{2} < \frac{1 - \sqrt{5}}{2} = \psi < 0 \quad \therefore\ |\psi| < 1 $$ したがって \( \lim_{n \to \infty} \psi^n = 0 \) であるから $$ \begin{align*} \lim_{n \to \infty} \frac{F_{n + 1}}{F_n} &= \frac{\displaystyle \lim_{n \to \infty} \phi^{n + 1} - \lim_{n \to \infty} \psi^{n + 1}}{\displaystyle \lim_{n \to \infty} \phi^n - \lim_{n \to \infty} \psi^n} = \frac{\displaystyle \lim_{n \to \infty} \phi^{n + 1} - 0}{\displaystyle \lim_{n \to \infty} \phi^n - 0} = \frac{\displaystyle \lim_{n \to \infty} \phi^{n + 1}}{\displaystyle \lim_{n \to \infty} \phi^n} \\ &= \lim_{n \to \infty} \frac{\phi^{n + 1}}{\phi^n} = \lim_{n \to \infty} \phi = \phi \end{align*} $$
これを利用すると $$ \begin{align*} \lim_{n \to \infty} y_n &= \lim_{n \to \infty} \frac{F_{n + 1} x + F_n}{F_n x + F_{n - 1}} = \lim_{n \to \infty} \frac{F_{n + 1} x + F_n}{F_n x + F_{n - 1}} \cdot \frac{\displaystyle \frac{1}{F_n}}{\displaystyle \frac{1}{F_n}} = \lim_{n \to \infty} \frac{\displaystyle \frac{F_{n + 1}}{F_n} x + 1}{\displaystyle x + \frac{F_{n - 1}}{F_n}} \\ &= \frac{\phi x + 1}{\displaystyle x + \frac{1}{\phi}} = \frac{\phi (\phi x + 1)}{\phi x + 1} = \phi \end{align*} $$
よって、\( \phi \) は関数 \( f(x) = 1 + \frac{1}{x} \) の吸引的不動点である。
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 テストを改良することができた。
