1.2 手続きとその生成するプロセス

視覚化大事ね。

1.2.1 線形再帰と反復

遅延演算が山なりに増えて減るのが再帰的プロセス。特に遅延演算の列の長さが入力に対して線形に増えていくと線形再帰的プロセス。

プロセスの長さが伸び縮みしないのが反復的プロセス。とくに必要なステップ数が線形に増えるのが線形反復的プロセス。反復的プロセスではプログラム変数がプロセスの完全な状態を保持している。

ふぅ。

問題1.9

前者の方は、再帰的プロセスでを生成する。

(+ 4 5)
==> (if (= 4 0) 5 (inc (+ (dec 4) 5)))
==> (inc (+ (dec 4) 5))
==> (inc (+ 3 5))
==> (inc (inc (+ (dec 3) 5)))
==> (inc (if (= 3 0) 5 (inc (+ (dec 3) 5))))
==> (inc (inc (+ (dec 3) 5)))
==> (inc (inc (+ 2 5)))
==> (inc (inc (if (= 2 0) 5 (inc (+ (dec 2) 5)))))
==> (inc (inc (inc (+ (dec 2) 5))))
==> (inc (inc (inc (+ 1 5))))
==> (inc (inc (inc (if (= 1 0) (inc (+ (dec 1) 5))))))
==> (inc (inc (inc (inc (+ (dec 1) 5)))))
==> (inc (inc (inc (inc (+ 0 5)))))
==> (inc (inc (inc (inc (if (= 0 0) 5 (inc (+ (dec 0) 5)))))))
==> (inc (inc (inc (inc 5))))
==> (inc (inc (inc 6)))
==> (inc (inc 7))
==> (inc 8)
==> 9

後者の方は、反復的プロセスを生成する。

(+ 4 5)
==> (if (= 4 0) 5 (+ (decr 4) (incr 5)))
==> (+ (decr 4) (incr 5))
==> (+ 3 6)
==> (if (= 3 0) 6 (+ (decr 3) (incr 6)))
==> (+ (decr 3) (incr 6))
==> (+ 2 7)
==> (if (= 2 0) 7 (+ (decr 2) (incr 7)))
==> (+ (decr 2) (incr 7))
==> (+ 1 8)
==> (if (= 1 0) 8 (+ (decr 1) (incr 8)))
==> (+ (decr 1) (incr 8))
==> (+ 0 9)
==> (if (= 0 0) 9 (+ (dec 0) (incr 9)))
==> 9

問題1.10

;; Ackermann関数
(define (A x y)
  (cond ((= y 0) 0)
        ((= x 0) (* 2 y))
        ((= y 1) 2)
        (else (A (- x 1)
                 (A x (- y 1))))))
(A 1 10) ==> 1024
(A 2 4) ==> 65536
(A 3 3) ==> 65536

式を置き換えていく。condで、条件が必ず#fもしくは到達しない項は消していくことにする。

(define (f n) (A 0 n))
==> (define (f n)
      (cond ((= n 0) 0)
            ((= 0 0) (* 2 n))
            ((= n 1) 2)
            (else (A (- 0 1)
                     (A 0 (- n 1))))))
==> (define (f n)
      (cond ((= n 0) 0)
            ((= 0 0) (* 2 n))))

fの値はは n=0のとき0, それ以外は、(* 2 n)になる。
よって、f(n) = 2n

(define (g n) (A 1 n))
==> (define (g n)
      (cond ((= n 0) 0)
            ((= 1 0) (* 2 n))
            ((= n 1) 2)
            (else (A (- 1 1)
                     (A 1 (- n 1))))))
==> (define (g n)
      (cond ((= n 0) 0)
            ((= n 1) 2)
            (else (A 0
                     (A 1 (- n 1))))))
(A 1 (- n 1)) = (g (- n 1))なので、
==> (define (g n)
      (cond ((= n 0) 0)
            ((= n 1) 2)
            (else (A 0 (g (- n 1))))))
==> (define (g n)
      (cond ((= n 0) 0)
            ((= n 1) 2)
            (else (f (g (- n 1))))))

gは2のn乗を返す.ただし、n=0の時は0を返す.
$g(0) = 0$
$g(n) = 2^n (n > 0)$

(define (h n) (A 2 n))
==> (define (h n)
      (cond ((= n 0) 0)
            ((= 2 0) (* 2 n))
            ((= n 1) 2)
            (else (A (- 2 1)
                     (A 2 (- n 1))))))
==> (define (h n)
      (cond ((= n 0) 0)
            ((= n 1) 2)
            (else (A 1
                     (A 2 (- n 1))))))
(A 2 (- n 1)) = (h (- n 1))なので、
==> (define (h n)
      (cond ((= n 0) 0)
            ((= n 1) 2)
            (else (A 1
                     (h (- n 1))))))
==> (define (h n)
      (cond ((= n 0) 0)
            ((= n 1) 2)
            (else (g (h (- n 1))))))

$h(0) = 0$
$h(1) =2$
$h(n) = 2^{h(n-1)} (n > 1)$
なので
$h(0) = 0$
$h(n) = 2^{2^{2^{...}}} (n > 0)$ (2はn個)

1.2.2 木構造 再帰

木構造再帰は呼び出し関係のグラフを書くと木構造になるのか。

両替の計算

ためしに、実行してみる。硬貨は日本のものを使う.

(define (count-change amount)
  (cc amount 6))

(define (cc amount kinds-of-coins)
  (cond ((= amount 0) 1)
        ((= amount 0) 1)
        ((or (< amount 0) (<= kinds-of-coins 0)) 0)
        (else (+ (cc amount (- kinds-of-coins 1))
                 (cc (- amount (first-denomination kinds-of-coins))
                     kinds-of-coins)))))

(define (first-denomination kinds-of-coins)
  (cond ((= kinds-of-coins 1) 1)
        ((= kinds-of-coins 2) 5)
        ((= kinds-of-coins 3) 10)
        ((= kinds-of-coins 4) 50)
        ((= kinds-of-coins 5) 100)
        ((= kinds-of-coins 6) 500)))
gosh> (count-change 10)
4
gosh> (count-change 20)
9
gosh> (count-change 100)
159

これだと、(count-change 1000) は計算に長い時間が掛かるので配列を使ってメモ化してみる。これで、計算量はO(n)にまで抑えられるはず。

(define (count-change amount)
  (define table (make-vector (+ amount 1)))
  (do ((i 0 (+ i 1)))
      ((>= i (vector-length table)))
    (vector-set! table i (make-vector 7 #f)))
  (cc amount 6 table))

(define (cc amount kinds-of-coins table)
  (cond ((= amount 0) 1)
        ((or (< amount 0) (<= kinds-of-coins 0)) 0)
        (else
         (if (vector-ref (vector-ref table amount) kinds-of-coins)
             (vector-ref (vector-ref table amount) kinds-of-coins)
             (let ((num
                    (+ (cc amount (- kinds-of-coins 1) table)
                       (cc (- amount (first-denomination kinds-of-coins))
                           kinds-of-coins
                           table))))
               (vector-set! (vector-ref table amount) kinds-of-coins  num)
               num)))))

問題1.11

再帰的プロセス

(define (f n)
  (if (< n 3)
      n
      (+ (f (- n 1))
         (* 2 (f (- n 2)))
         (* 3 (f (- n 3))))))

反復的プロセス

(define (f n)
  (define (f-iter i a b c)
    (if (= i 2)
        c
        (f-iter (- i 1)
                b
                c
                (+ c (* 2 b) (* 3 a)))))
  (if (< n 3)
      n
      (f-iter n 0 1 2)))

問題1.12

左端からx番目、上からy番目(0から数える)。

(define (pascal-triangle x y)
  (cond ((or (< x 0) (< y 0) (> x y)) 0)
        ((or (= x 0) (= x y)) 1)
        (else (+ (pascal-triangle (- x 1) (- y 1))
                 (pascal-triangle x (- y 1))))))

(define (pascal-triangle x y)
  (cond ((or (< x 0) (< y 0) (> x y)) 0)
        ((or (= x 0) (= x y)) 1)
        (else (+ (pascal-triangle (- x 1) (- y 1))
                 (pascal-triangle x (- y 1))))))
(do ((i 0 (+ i 1)))
    ((> i 10))
  (do ((j 0 (+ j 1)))
      ((> j 10))
    (let ((n (pascal-triangle i j)))
      (format #t "~3@a " (if (= n 0) #\space n))))
  (newline))
  1   1   1   1   1   1   1   1   1   1   1 
      1   2   3   4   5   6   7   8   9  10 
          1   3   6  10  15  21  28  36  45 
              1   4  10  20  35  56  84 120 
                  1   5  15  35  70 126 210 
                      1   6  21  56 126 252 
                          1   7  28  84 210 
                              1   8  36 120 
                                  1   9  45 
                                      1  10 
                                          1

なんとか三角形に見えるかな。

問題1.13

(define phi (/ (+ 1 (sqrt 5)) 2))
(define (fib n)
  (inexact->exact (round (/ (expt phi n) (sqrt 5)))))

おお、確かにこれで求められる.だけど、桁が大きくなると誤差が出て誤った数になるね.

証明：
$\phi = \frac{1 + \sqrt[]{5}}{2}$
$\psi = \frac{1 - \sqrt[]{5}}{2}$
として、まず
$Fib(n) = \frac{\phi^n - \psi^n}{\sqrt[]{5}}$ 　… A
を帰納法を用いて証明する。
n=0のとき、
$Fib(0) = \frac{\phi^0 - \psi^0}{\sqrt[]{5}} = 0$
となり、成り立つ.
n=1のとき、
$Fib(1) = \frac{\phi^1 - \psi^1}{\sqrt[]{5}} = 1$
となり、成り立つ.
n=kのとき、n=k-1, n=k-2で成り立つとすると、
$Fib(k) = Fib(k-1) + Fib(k-2) = \frac{\phi^{k-1} - \psi^{k-1}}{\sqrt[]{5}} + \frac{\phi^{k-2} - \psi^{k-2}}{\sqrt[]{5}}$
でこれを計算していくと、
$Fib(k) = \frac{\phi^k - \psi^k}{\sqrt[]{5}}$
となる。よって、Aは証明された.

さて、
$f(n) = \frac{\phi^n}{\sqrt[]{5}}$
とすると、
$\left| \frac{\phi^n - \psi^n}{\sqrt[]{5}} - \frac{\phi^n}{\sqrt[]{5}} \right| = \left| \frac{\psi^n}{\sqrt[]{5}} \right|$ 　… B
となる。 $|\psi| < 0.5$ , $\sqrt[]{5} > 2$ より、B < 0.5となる。よって、Fib(n)はf(n)に最も近い整数である.

1.2.3 増加の程度

関数の増え型のの記法は、Wikipediaによると、
f(x) = O(g(x))のとき、f(x)はg(x)の定数倍以下。f(x) = Ω(g(x))のとき、f(x)はg(x)の定数倍以上。f(x) = Θ(g(x))のとき、f(x)はg(x)の定数倍の範囲に収まる.これは、f(x) = O(g(x))かつ f(x) = Ω(g(x)) であること等価。
という感じ.

問題1.14

http://d.hatena.ne.jp/taisuke_h/20070303 を使った。

(cc 11 4) = 4
|
+-(cc 11 3) = 4
| |
| +-(cc 11 2) = 3
| | |
| | +-(cc 11 1) = 1
| | | |
| | | +-(cc 11 0) = 0
| | | |
| | | `-(cc 10 1) = 1
| | |   |
| | |   +-(cc 10 0) = 0
| | |   |
| | |   `-(cc 9 1) = 1
| | |     |
| | |     +-(cc 9 0) = 0
| | |     |
| | |     `-(cc 8 1) = 1
| | |       |
| | |       +-(cc 8 0) = 0
| | |       |
| | |       `-(cc 7 1) = 1
| | |         |
| | |         +-(cc 7 0) = 0
| | |         |
| | |         `-(cc 6 1) = 1
| | |           |
| | |           +-(cc 6 0) = 0
| | |           |
| | |           `-(cc 5 1) = 1
| | |             |
| | |             +-(cc 5 0) = 0
| | |             |
| | |             `-(cc 4 1) = 1
| | |               |
| | |               +-(cc 4 0) = 0
| | |               |
| | |               `-(cc 3 1) = 1
| | |                 |
| | |                 +-(cc 3 0) = 0
| | |                 |
| | |                 `-(cc 2 1) = 1
| | |                   |
| | |                   +-(cc 2 0) = 0
| | |                   |
| | |                   `-(cc 1 1) = 1
| | |                     |
| | |                     +-(cc 1 0) = 0
| | |                     |
| | |                     `-(cc 0 1) = 1
| | |
| | `-(cc 6 2) = 2
| |   |
| |   +-(cc 6 1) = 1
| |   | |
| |   | +-(cc 6 0) = 0
| |   | |
| |   | `-(cc 5 1) = 1
| |   |   |
| |   |   +-(cc 5 0) = 0
| |   |   |
| |   |   `-(cc 4 1) = 1
| |   |     |
| |   |     +-(cc 4 0) = 0
| |   |     |
| |   |     `-(cc 3 1) = 1
| |   |       |
| |   |       +-(cc 3 0) = 0
| |   |       |
| |   |       `-(cc 2 1) = 1
| |   |         |
| |   |         +-(cc 2 0) = 0
| |   |         |
| |   |         `-(cc 1 1) = 1
| |   |           |
| |   |           +-(cc 1 0) = 0
| |   |           |
| |   |           `-(cc 0 1) = 1
| |   |
| |   `-(cc 1 2) = 1
| |     |
| |     +-(cc 1 1) = 1
| |     | |
| |     | +-(cc 1 0) = 0
| |     | |
| |     | `-(cc 0 1) = 1
| |     |
| |     `-(cc -4 2) = 0
| |
| `-(cc 1 3) = 1
|   |
|   +-(cc 1 2) = 1
|   | |
|   | +-(cc 1 1) = 1
|   | | |
|   | | +-(cc 1 0) = 0
|   | | |
|   | | `-(cc 0 1) = 1
|   | |
|   | `-(cc -4 2) = 0
|   |
|   `-(cc -9 3) = 0
|
`-(cc -14 4) = 0

金額をn,硬貨の種類の数をdとする。
スペースの増加の程度は、上記の木構造の下方向の深さがΘ(d), 横方向の深さがΘ(n)になるので、Θ(n + d)。この場合は硬貨の種類は固定なのでΘ(n)になる。
次に、ステップ数。f(x, d)は金額x硬貨の種類がdの時のステップ数.以下、非常にいい加減な計算の元にステップ数の増加の程度を計算してみる.

f(x, 0) = 1
f(x, 1) = 1 + f(x-1, 1) + f0(x, 0)
        = x + 1
f(x, 2) = 1 + f(x-5, 2) +  f(x, 1)
        = x/5 + f(x, 1) + f(x-5, 1) + f(x-10, 1) + ...
        = x/5 + (x/5 + x + x-5 + x-10 + ...)
        = x/5 + (x/5 + Θ(x^2))
        = Θ(x^2)

f(x, 3) = 1 + f(x-10, 3) + f(x, 2)
        = x/10 + f(x-10, 2) + f(x-20, 2) + ...
        = f/10 + Θ(x^3)
f(x, 4) = 1 + f(x-25, 4) + f(x, 3)
        = Θ(x^4)
f(x, 5) = 1 + f(x-50, 5) + f(x, 4)
        = Θ(x^5)

というわけで、ステップの増加の程度はΘ(n^d), この場合はΘ(n^5)。（あんまり自信が無い）

問題1.15

(define (cube x) (* x x x))
(define (p x) (- (* 3 x) (* 4 (cube x))))
(define (sine angle)
  (if (not (> (abs angle) 0.1))
      angle
      (p (sine (/ angle 3.0)))))

a. 5回
b.
sineの処理は再帰的にsineを呼び出している以外はすべてΘ(1)の処理である.sineの呼出は末尾再帰では無い.よって、再帰の深さの増加の程度がスペース、ステップ数の増加の程度となる。sineの再帰が1深くなるごとにangleは1/3になる。
以下の式のlogの底は3。
再帰の深さをnとすると、nは
angle / 3^n <= 0.1
となる最小の正の数である。
angle <= 0.1 * 3^n
log(angle) <= log(0.1) + log(3^n)
log(angle) <= log(0.1) + n
log(angle) - log(0.1) <= n
よって、スペース、ステップ数の増加の程度はΘ(log(angle))となる。

1.2.4　べき乗

fast-exptの乗算の回数は原書の補足によると
$\lfloor \log[2](n) \rfloor + number-of-1-in-binary$
となっている。
最下位ビットが１の時(奇数の時）それを0にして、fast-exptを呼ぶ.
最下位ビットが０のとき（偶数の時）右に１ビットシフトしてfast-exptを呼ぶ.
という処理を0になるまで行うので、奇数になる回数は2進数での１の数.０より大きい偶数になる回数はnに必要なビット数-1（式のlogの部分）となる。

問題1.16

fast-expt-iter-subの呼出においてab^nが常に一定になっている。これが、不変量か。

(define (fast-expt-iter b n)
  (fast-expt-iter-sub b n 1))
(define (fast-expt-iter-sub b n a)
  (cond ((= n 0) a)
        ((even? n) (fast-expt-iter-sub (* b b) (/ n 2) a))
        (else (fast-expt-iter-sub b (- n 1) (* b a)))))

問題1.17

(define (halve n) (/ n 2))
(define (double n) (* n 2))
;; 再帰版
(define (fast-mult-rec a b)
  (cond ((= b 0) 0)
        ((even? b)
         (double (fast-mult-rec a (halve b))))
        (else
         (+ a (fast-mult-rec a (- b 1))))))
;; 反復版
; ab+cが不変になるようにする
; 偶数の時は a * b + c = 2a * b/2 + c
; 奇数の時は a * b + c = a * (b-1) + a+c
(define (fast-mult-iter a b)
  (define (fast-mult-iter-sub a b c)
    (cond ((= b 0) c)
          ((even? b) (fast-mult-iter-sub (double a) (halve b) c))
          (else (fast-mult-iter-sub a (- b 1) (+ a c)))))
  (fast-mult-iter-sub a b 0))

問題1.18

あれ、1.17でやっちゃった。コピーしておく。

; ab+cが不変になるようにする
; 偶数の時は a * b + c = 2a * b/2 + c
; 奇数の時は a * b + c = a * (b-1) + a+c
(define (fast-mult-iter a b)
  (define (iter a b c)
    (cond ((= b 0) c)
          ((even? b) (iter (double a) (halve b) c))
          (else (iter a (- b 1) (+ a c)))))
  (iter a b 0))

問題1.19

p', q'をp, qを用いて表す.
本文中のaの変換TpqをTa(a, b, p, q).bの変換をTb(a, b, p, q)と書く.
Ta(a, b, p, q) = bq + aq + ap
Tb(a, b, p, q) = bp + aq
変換を2回適用する.

Ta(Ta(a, b, p, q), Tb(a,b, p, q), p, q)
  = Ta(bq + aq + ap, bp + aq, p, q)
  = (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
  = bpq + aqq + bqq + aqq + apq + bpq + apq + app
  = b(2pq + qq) + a(2qq + 2pq + pp)
  = (a + b)(2pq + qq) + a(qq + pp)
  = Ta(a, b, (qq + pp), (2pq + qq))

Tb(Ta(a, b, p, q), Tb(a,b, p, q), p, q)
  = Ta(bq + aq + ap, bp + aq, p, q)
  = (bp + aq)p + (bq + aq + ap)q
  = bpp + apq + bpq + aqq + apq
  = b(pp + pq) + a(2pq + pq)
  = Tb(a, b, (qq + pp), (2pq + qq))

よって、

p' = p^2 + q^2
q' = 2pq + q^2

次に、プログラム

(define (fib n)
  (fib-iter 1 0 0 1 n))
(define (fib-iter a b p q count)
  (cond ((= count 0) b)
        ((even? count)
         (fib-iter a
                   b
                   (+ (* p p) (* q q))
                   (+ (* 2 p q) (* q q))
                   (/ count 2)))
        (else
         (fib-iter (+ (* b q) (* a q) (* a p))
                   (+ (* b p) (* a q))
                   p
                   q
                   (- count 1)))))

上記の式は行列で考えるとわかりやすい.

Fib(n) = (a, b) ((p+q), q)^n
                 (q,   , p)

プログラムはこの式を計算しているに他ならない.
2x2行列の2乗は
(pp+2pq+2qq, 2pq+qq) = (p' + q', q')
(2pq+qq, qq+pp,) (q' , p')
となるので、pとqを変換するだけだ.

1.2.5 最大公約数

問題1.20

ifでbが評価される.

(gcd a b)
(if (= 40 0) 206 (gcd 40 (rem 206 40)))
(gcd 206, 40)
= (if (= 40 0) 206 (gcd 40 (rem 206 40)))
= (gcd 40 (rem 206 40))
= (if (= (rem 206 40) 0) 40 (gcd (rem 206 40) (rem 40 (rem 206 40))))
1
= (if (= 6 0) 40 (gcd (rem 206 40) (rem 40 (rem 206 40))))
= (gcd (rem 206 40) (rem 40 (rem 206 40)))
= (if (= (rem 40 (rem 206 40)) 0)
      (rem 206 40)
      (gcd (rem 40 (rem 206 40)) (rem (rem 206 40) (rem 40 (rem 206 40)))))
= (if (= 4 0)
      (rem 206 40)
      (gcd (rem 40 (rem 206 40)) (rem (rem 206 40) (rem 40 (rem 206 40)))))
2
= (gcd (rem 40 (rem 206 40)) (rem (rem 206 40) (rem 40 (rem 206 40))))
= (if (= (rem (rem 206 40) (rem 40 (rem 206 40))) 0)
      (rem 40 (rem 206 40))
      (gcd (rem (rem 206 40) (rem 40 (rem 206 40)))
           (rem (rem 40 (rem 206 40)) (rem (rem 206 40) (rem 40 (rem 206 40))))))
4
= (if (= 2 0)
      (rem 40 (rem 206 40))
      (gcd (rem (rem 206 40) (rem 40 (rem 206 40)))
           (rem (rem 40 (rem 206 40)) (rem (rem 206 40) (rem 40 (rem 206 40))))))

= (gcd (rem (rem 206 40) (rem 40 (rem 206 40)))
       (rem (rem 40 (rem 206 40)) (rem (rem 206 40) (rem 40 (rem 206 40)))))
= (if (= (rem (rem 40 (rem 206 40)) (rem (rem 206 40) (rem 40 (rem 206 40)))) 0)
      (rem (rem 206 40) (rem 40 (rem 206 40)))
      (gcd (rem (rem 40 (rem 206 40)) (rem (rem 206 40) (rem 40 (rem 206 40))))
           (rem (rem (rem 206 40) (rem 40 (rem 206 40)))
                (rem (rem 40 (rem 206 40))
                     (rem (rem 206 40) (rem 40 (rem 206 40)))))))
7
= (if (= 0 0)
      (rem (rem 206 40) (rem 40 (rem 206 40)))
      (gcd (rem (rem 40 (rem 206 40)) (rem (rem 206 40) (rem 40 (rem 206 40))))
           (rem (rem (rem 206 40) (rem 40 (rem 206 40)))
                (rem (rem 40 (rem 206 40))
                     (rem (rem 206 40) (rem 40 (rem 206 40)))))))
= (rem (rem 206 40) (rem 40 (rem 206 40)))
4
= 2

正規順序評価では、1+2+4+7+4=18回.同じ式が一回しか評価されないとすると4回。
作用的順序では4回。

1.2.6 例: 素数制のテスト

√nまでに除数を持つというのは、除数をaとbとするとa>=√nの時、b<=√nになるからだな.

(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))))
(use srfi-27) ; 乱数のために必要
(define (fermat-test n)
  (define (try-it a)
    (= (expmod a n n) a))
  (try-it (+ 1 (random-integer (- n 1)))))

(define (fast-prime? n times)
  (cond ((= times 0) #t)
        ((fermat-test n) (fast-prime? n (- times 1)))
        (else #f)))

問題1.21

(define (smallest-divisor n)
  (find-divisor n 2))
(define (divide? a b)
  (= 0 (remainder a b)))

(define (find-divisor n test-divisor)
  (cond ((< n (square test-divisor)) n)
        ((divide? n test-divisor) test-divisor)
        (else (find-divisor n (+ 1 test-divisor)))))
(define (prime? n)
  (= n (smallest-divisor n)))

(smallest-divisor 199) ==> 199
(smallest-divisor 1999) ==> 1999
(smallest-divisor 19999) ==> 7

問題1.22

問題に書いてあるプログラムは使いにくいので中身をいじりました.時間計測にはGaucheの機能を使っています.

(define (runtime)
  (time->seconds (current-time)))

(define-syntax %time
  (syntax-rules ()
    ((_ count f)
     (let ((time 0))
       (do ((i 0 (+ i 1)))
           ((>= i count))
         (let* ((start (runtime))
                (v f)
                (end (runtime)))
           (set! time (+ time (- end start)))))
       (/ time count)))))

(define (start-prime-test n)
  (if (prime? n)
      (or (report-prime n (%time 10 (prime? n)))
          #t)
      #f))

(define (report-prime n elapsed-time)
  (display n)
  (display " *** ")
  (display elapsed-time)
  (newline))

(define (search-for-primes n)
  (define (iter n cnt)
    (cond ((= cnt 0))
          (else
           (let ((t (start-prime-test n)))
             (iter (+ n 2)
                   (- cnt (if t 1 0)))))))
  (iter (if (even? n) (+ n 1) n)  3))

(do ((i 1000 (* i 10)))
    ((> i 100000000))
  (search-for-primes i)
  (newline))

実行結果.時間は秒単位。時間は10回実行した平均を取っている.

1009 *** 2.3894309997558593e-5
1013 *** 2.427816390991211e-5
1019 *** 2.254009246826172e-5

10007 *** 7.630348205566407e-5
10009 *** 7.204771041870117e-5
10037 *** 6.950616836547851e-5

100003 *** 2.2489070892333984e-4
100019 *** 2.3761510848999023e-4
100043 *** 2.2112369537353516e-4

1000003 *** 7.245206832885742e-4
1000033 *** 5.836963653564453e-4
1000037 *** 5.877399444580078e-4

10000019 *** 0.0015128350257873535
10000079 *** 0.0019036698341369628
10000103 *** 0.0015006232261657714

100000007 *** 0.004454286098480225
100000037 *** 0.004819636344909668
100000039 *** 0.004433591365814209

問題1.23

(define (next n)
  (if (= n 2)
      3
      (+ n 2)))
(define (find-divisor n test-divisor)
  (cond ((< n (square test-divisor)) n)
        ((divide? n test-divisor) test-divisor)
        (else (find-divisor n (next test-divisor)))))

前回の実行時間と今回改良したものの実行時間の比較.数が大きくなると実行時間は半分に近づいている.しかし、全体として0.5よりも数値が高いのはCUPのキャッシュの関係なのだろうか?

N	Naive [μs]	Improved [μs]	Naive / Improved
1009	23.89	15.01	0.63
1013	24.28	12.89	0.53
1019	22.54	17.43	0.77
10007	76.30	40.17	0.53
10009	72.05	38.90	0.54
10037	69.51	39.49	0.57
100003	224.89	126.66	0.56
100019	237.62	129.26	0.54
100043	221.12	129.59	0.59
1000003	724.52	430.63	0.59
1000033	583.70	388.76	0.67
1000037	587.74	382.82	0.65
10000019	1512.84	766.93	0.51
10000079	1903.67	1001.47	0.53
10000103	1500.62	833.84	0.56
100000007	4454.29	2300.75	0.52
100000037	4819.64	2305.68	0.48
100000039	4433.59	2302.61	0.52

問題1.24

$\lfloor \log[2](n) \rfloor + number-of-1-in-binary$
fast-prime?の計算量はこの式なのでだいたい、底が2の2log(n)になっている。nが1000から1000000になると 2log(1000000) / 2log(1000) = 2 倍になると予想する.では、やってみよう.
start-prime-testを以下のように変える。

(define (start-prime-test n)
  (if (fast-prime? n 1)
      (or (report-prime n (%time 10 (fast-prime? n 1)))
          #t)
      #f))

実行結果は

N	time [μs]
1009	19.93
1013	19.34
1019	20.05
10007	22.08
10009	22.89
10037	26.56
100003	46.73
100019	52.86
100043	53.00
1000003	63.20
1000033	63.75
1000037	65.80
10000019	83.33
10000079	84.21
10000103	92.41
100000007	101.09
100000037	91.29
100000039	97.92

になった。1000の時を、20、1000000の時を、64とすると、3.2倍になっている.予想と全然違う.。でも、100000と100000000を比べると、大体2倍になっているからいいのかな?というか数値がバラバラ過ぎて、よくわからない.

問題1.25

返す値は同じになる.ただし、非常に大きな数の計算をしないとならなくなるので効率が悪い。二つの実行時間を計測してみた。

(use gauche.time)
(define (fast-expt base exp)
  (cond ((= exp 0) 1)
        ((even? exp) (fast-expt (* base base) (/ exp 2)))
        (else (* base (fast-expt base (- exp 1))))))

(define (expmod base exp m)
  (remainder (fast-expt base exp) m))

(define (square n) (* n n))
(define (fast-expmod base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (remainder (square (fast-expmod base (/ exp 2) m))
                    m))
        (else
         (remainder (* base (fast-expmod base (- exp 1) m))
                    m))))
(define-syntax %time
  (syntax-rules ()
    ((_ . body)
     (let1 cnt (make <user-time-counter>)
       (with-time-counter cnt . body)
       (exact->inexact (time-counter-value cnt))))))

(format #t "~10a~10a~10a\n" "N" "expmod" "fast-expmod")
(do ((n 10 (* n 2)))
    ((> n 1e7))
  (format #t "~10a~10a~10a\n"
          n
          (%time (expmod 2 n 5))
          (%time (fast-expmod 2 n 5))))

実行結果

N	expmod	fast-expmod
10	0.0	0.0
20	0.0	0.0
40	0.0	0.0
80	0.0	0.0
160	0.0	0.0
320	0.0	0.0
640	0.0	0.0
1280	0.0	0.0
2560	0.0	0.0
5120	0.0	0.0
10240	0.0	0.0
20480	0.0	0.0
40960	0.02	0.0
81920	0.06	0.0
163840	0.25	0.0
327680	1.0	0.0
655360	3.97	0.0
1310720	15.84	0.0
2621440	63.66	0.0
5242880	254.37	0.0

というように、乗数が大きくなるとexpmodの方は著しく計算に時間が掛かるようになる.

問題1.26

本問の案だとexpが偶数の時にexpmodを2回呼んでしまっている。
以前のexmpdの偶数での乗算の回数をf(n)とすると、(logの底は2)

f(n) = 1 + f(n/2)
f(8) = 1 + f(4)
     = 1 + 1 + f(2)
     = 1 + 1 + 1 + f(1)
     = 1 + 1 + 1 + 1
     = 4
f(n) = log(n) + 1

本問の偶数での乗算の回数をg(n)とすると

g(n) = 1 + 2*g(n/2)
g(8) = 1 + 2g(4)
     = 1 + 2(1 + 2g(2))
     = 1 + 2(1 + 2(1 + 2g(1)))
     = 1 + 2(1 + 2(1 + 2))
     = 15
g(n) = 1 + 2 + 4 + ... + 2^log(n)
     = 1 + 2 + 4 + ... + 2^log(n)
     = 2^(log(n)+1) - 1
     = 2^log(n) * 2 - 1
     = n * 2 - 1

よって、g(n) = Θ(n)となる。上の計算は大体の計算です。

問題1.27

(define (square n) (* n n))
(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 (check n)
  (define (iter a)
    (if (< a n)
        (and (= a (expmod a n n))
             (iter (+ a 1)))
        #t))
  (iter 1))

(and
 (check 561)
 (check 1105)
 (check 1729)
 (check 2465)
 (check 2821)
 (check 6601)) ==> #t

すべてのチェックを通過した.

問題1.28

問題文の通りに0シグナルを出す。

(use srfi-27)
(define (square n) (* n n))

(define (expmod2 base exp m)
  (cond ((= exp 0) 1)
        ((even? exp)
         (let* ((x (expmod2 base (/ exp 2) m))
                (r (remainder (square x) m)))
           (if (and (< 1 x)
                    (< x (- m 1))
                    (= r 1))
               0
               r)))
        (else
         (remainder (* base (expmod2 base (- exp 1) m))
                    m))))
(define (miller-rabin-test n)
  (define (try-it a)
    (= (expmod2 a (- n 1) n) 1))
  (try-it (+ 1 (random-integer (- n 1)))))
(define (miller-prime? n times)
  (cond ((= times 0) #t)
        ((miller-rabin-test n)
         (miller-prime? n (- times 1)))
        (else #f)))

どれくらい、判定漏れがあるかチェックしてみる

;; チェック
(define (try-miller start end times)
  (when (< start end)
    (let ((p (prime? start))
          (m (miller-prime? start times)))
      (cond ((and p (not m))
             (format #t "~8@a ->\n" start)) ; 逸脱
            ((and (not p) m)
             (format #t "~8@a <-\n" start)))) ; 混入
    (try-miller (+ start 1) end times)))

結果は

(try-miller 2 1000 1)
        38 <-
        70 <-
        76 <-
        91 <-
       262 <-
       339 <-
       427 <-
       532 <-
       889 <-
       976 <-

(try-miller 2 100000 2)
       2047 <-


(try-miller 2 100000 3)
nothing