問題回答

問1.28
内容的に特に新しいことがなさそうなのでパス。

本文

1.3 高階手続きによる抽象
手続きを行う手続き : 高階手続き
事前に定義された手続きを引数にとって処理をおこなう手続きのこと。

高階関数を使って総和の計算などを抽象的に扱うことができる。

問題回答

問1.29
一応それっぽく書くことができたものの、結果がおかしい。
悩んでいたところでタイムアウトなので、明日に持ち越し。

(define (simpson f a b n)
  (define (h)
    (/ (- b a) n))

  (define (y f k)
    (f (+ a (* k (h)))))
    
  (define (next i)
    (+ i 1))

  (define (sum k)
    (cond ((> k n)
           0)
          ((= k 0)
           (+ (y f k) (sum (next k))))
          ((= (mod n 2) 0)
           (+ (* (y f k) 2) (sum (next k))))
          (else
           (+ (* (y f k) 4) (sum (next k))))))

  (* (/ (h) 3) (sum 0)))

Meadow から Emacs for Windows への移行はちょっと大変そうなので、とりあえずは置いておくことにしました(問題の先送り….)。
それにあわせて DrRacket 導入も見送りで、結局当分は Meadow+Gauche という元鞘に戻りそうな予感。時が満ちたら再チャレンジします。

SICP 本体は、どうも Fermat テストに関する理解ができていないせいでこのあたりの問題に着手すると手が止まってしまいます。
飛ばすのは簡単なんですが、Fermat テストに関係ない本質的なところまで飛ばしてしまうことになるので悩みどころです。
完璧に理解するのは無理として、一度 Fermat テスト周りのコードは全て写経して完全に動くソースコードを作ってみます。

本文

Fermat テスト再考

Fermat テストの理解もざっくりしすぎなので、一度整理。

法として合同 :
二つの数を n で割った時、同じ徐余を持てば、 n を法として合同という。

例 : 8 と 11
それぞれ 3 で割ると余りは 2 なので、8 と 11 は 3 を法として合同である。
また、この時の余り 2 は、3 を法とした 8(11) の徐余もしくは 3 を法とした 8(11) である。

Fermat の小定理 :
n を素数、a を n より小さい正の任意の整数とすると、a の n 乗は n を法として a と合同である。

例 : n=11, a=3
3^11 = 177147
177147 mod 11 = 3
3 mod 11 = 3 → 11 を法として合同

n = 12, a=3
3^12 = 531441
531441 mod 12 = 9 → 12 を法として合同でない

Fermat テストコード

(define (random num)
  (use srfi-27)
  (* (random-integer num) 1.0))


(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 #f)))

問題解答

問1.27
Carmichael 数の確認

 (checkCarmichael n)
  (define (check-iter a n)
    (cond ((= a (- n 1)) #t)
          ((= (expmod a n n) a)
           (check-iter (+ a 1) n))
          (else #f)))

  (check-iter 1 n)

gosh> (checkCarmichael 561)
#t
gosh> (checkCarmichael 1105)
#t
gosh> (checkCarmichael 1729)
#t
gosh> (checkCarmichael 2465)
#t
gosh> (checkCarmichael 2821)
#t
gosh> (checkCarmichael 6601)
#t

例えば 2821 = 7 * 13 * 31。
たしかにだまされている。

問題解答

問1.26

(* (expmod base (/ exp 2) m)
     (expmod base (/ exp 2) m))

この部分で、expmod が 2 度呼ばれてしまうために (/ exp 2) と、 exp の減少率を半分にしている意味がなくなってしまうため。

Racket 環境構築

問1.27 をやるにあたってコードを書かないといけなかったので、ここで Racket を使ってみるために環境構築。
主に「計算機プログラムの構造と解釈」のためのプログラミング環境を参考に環境構築をおこないます。

ただ、SICP Support for DrRacketの環境構築が私のところ(Racket 6.1.1 x86_64)ではうまくいっていなさそうです。最初のインストール時にエラーメッセージが出たり、Language の選択で SICP が対象にならないです。 #lang planet neil/sicp と入力すれば実行できているっぽいもののちょっと気持ち悪いので、 元記事と同じ 5.3.3 を使用することにしました。

5.3.3+SICP Support でもインストール時にエラーがボロボロ出たんですが、メニューからは選ぶことができるようになっているしとりあえずこれで様子を見てみて、やっぱりおかしそうなら 32bit 版を試したりしてみます。

ちょっと歩みが遅すぎて埒が明かなさそうなので問題に関する方針を変えてみることにします。

これまでは一問一問結構粘って、結果一日一問やったら終わりみたいな感じになってしまっていたので、ちょっとサクサク行くようにしてみます。Lisp のエラーとかそういうところでハマって時間を取られてしまったら答えあわせをしてしまって、写経しつつ自分のコードと比較をして理解をするようにします。

また、問題も一問で最大30分を目安にしてそれでもダメだったら諦めて答えを見ながら理解を進めるようにします。

Lisp 関連のエラーで悩まされて時間を潰すことがほとんどだったので、これで随分進みは良くなるんじゃないかな～?あと、ちょっと Gauche+Meadow の環境が辛いので別の環境を試したほうがいいのかも。 SICP コースだと Racket がお勧めっぽいので、そちらも試してみます。

問題回答

問1.24
gauche には random が無いようなので先人の知恵を拝借。
また、Gauche では true/false は #t/#f のようなのでこれは読み替えて使用。

回答からコードを持ってきて実行しても、素数のはずなのに素数ではないと判定されてしまう。何故～?

gosh> (fast-prime? 1000000000000037 10)
#f

仕方ないので解説を読む。

The timing numbers from the Fermat test start out looking pretty poor compared to what we've seen previously. This has mostly to do with the completely arbitrary number of random values I chose to test each prime with. As we increase the size of the numbers we're testing, you can see that the time using the Fermat test barely increases at all. We can verify that the time increase is logarithmic by observing that as the numbers under test increase by a factor of 10, the ratio column increases by a factor of roughly three. This logarithmic characteristic is due to the fact that the expmod procedure used in fermat-test has a logarithmic time complexity.

値が小さいときには非常に芳しくない結果となり、値が大きくなっても log10 ではなく log3 程度の効率化ししかなっていない。
これは、expmod が対数時間複雑性を持っているから、とのこと(P.29 の頭で説明されていた)。

問1.25
先人の知恵を拝借。

注釈 46 に解説があるのであわせて読む。

巨大な整数を扱うのはコストがかかるため、その影響が出てくるとのこと。なるほどー。

今日は方針を決めたりそれをメモしたりするのに時間をとられたのでここまで。よくわからないままウンウン悩むよりはサクサクと
答えを見ながら理解したほうがいい気はするので当分この流れでいってみる。

問題解答

問1.23
next の定義。

(define (next n)
  (if (= n 2)
      3
      (+ n 2)))

next を使って find-divisor をかきかえる

(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)))))

処理時間の比較。最近の計算機だと処理速度が速すぎて小さい値だと正確な処理時間が測れないようなので、大きめの値を使用。

改良前
1000000021 *** 20519#
100000000003 *** 267726#

改良後
1000000021 *** 18565#
100000000003 *** 164153#

処理時間はそれぞれ 0.9 倍、 0.61 倍となっており、 半分にはならない。

明確な理由がわからなかったので先人の知恵を借りる。
計算以外の、処理系のオーバーヘッドが大きな影響を及ぼしているっぽい。言われてみれば確かにそうか。

本文

1.2.6 例:素数性のテスト

Fermat テスト : ある数が素数か否かを確率的に求める。

テストに失敗すれば素数ではなく、テストに成功すれば高い確率で素数となる。更にテストを繰り返すと正確性が上がっていく。

このようなテストの存在が確率的アルゴリズムとして知られるようになった。

問題解答

問1.21
本のプログラムを写経して実行するだけ。

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

問1.22
gauche には runtime が無いようなので、まずはそれを解決。

prime? などは本をそのまま流用して、search-for-primes はこのような感じ。

(define (search-for-primes n)
  (define (search-for-primes-iter i k)
    (cond ((= k 0))
          ((prime? i)
           (display i)
           (newline)
           (search-for-primes-iter (+ i 2) (- k 1)))
          (else (search-for-primes-iter (+ i 2) k))))
                                
  (search-for-primes-iter n 3))

入力された値の偶奇判断とかモロモロ盛り込もうとするともっと複雑になってしまいそう。
あと、Meadow 上で gosh を走らせてコードを実行すると Stack Trace がきちんと表示されなかったり、ステップ実行できなかったりしてエラーが起きたときの対応が辛い。もっとちゃんとした環境を作らないと効率が悪すぎる。

問題解答

問1.17, 1.18
1.17 をやっていたらついでに 1.18 もやってしまったのでまとめて回答。

(define (even? n)
  (= (remainder n 2) 0))

(define (double a)
  (+ a a))

(define (halve a)
  (/ a 2))

(define (mult a n)
  (define (mult-iter a n t)
    (cond ((= n 0) 0)
          ((= n 1) a)
          ((even? n) (mult-iter (double a) (halve n) t))
          (else (+ (mult-iter a (- n 1) t) a))))

  (mult-iter a n 0))

正直、ビット演算とシフトを使いたい。。。

問1.19
プログラム的にやることはこれまでと変わらないのと、Fibonacchi 数の計算についての内容が主なのでパス。

本文

1.2.5 最大公約数
Euclid のアルゴリズムを使用して最大公約数(GCD)を求める

問題解答

問1.20

(gcd 206 40)
((gcd 40 (remainder 206 40)))
(((gcd (remainder 206 40) (remainder 40 (remainder 206 40))))))
:
:

正規順序評価を用いると gcd が延々展開されてしまって正常に評価がおこなわれない
→答え合わせをしたらそんなことなかった。めんどくさがらずにきちんと展開すればできた。

作用的順序評価は以下の通り

(gcd 206 40)
(gcd 40 (remainder 206 40))
(gcd 40 6)
(gcd 6 (remainder 40 6))
(gcd 6 4)
(gcd 4 (remainder 6 4))
(gcd 4 2)
(gcd 2 (remainder (4 2))
(gcd 2 0)

よって、remainder は 4 回呼ばれる。

問題解答

問1.15 a
(sine 12.15)
(sine 4.05)
:

と、sine が呼ばれる度に angle が 1/3 になり、それが 0.1 未満になるまで続くのだから

12.15*(1/3^n) < 0.1 が解ければよい。 121.5 < 3^n log(121.5) < n(log(3)) log(121.5)/log(3) < n 4.369 < n より、n は 5。 問1.15 b 反復的操作のため、スペースは O(1) ステップ数は n が三倍になると一つ増えるので O(log3(n))

本文

1.2.4 べき乗
再帰的プロセス、反復的プロセスを使用した場合のオーダーと、更に逐次平方を利用した場合のオーダーの比較。
O(log n) にできた場合の効果。

問題解答

問1.16
自分の書いたコードは以下。

(define (fast-expt b n)
  (define (even? v)
    (= (remainder v 2) 0))

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

  (expt-iter 1 b n))

答え合わせをして気づいたのだが、 even? が真だった場合に評価されるところが末尾最適化されないので解答としては不正解。

生活習慣が(今のところ)ガラッと変わって夜に寝るのがものすごく早くなったので、当分は朝に勉強時間を設けるようにしてみます。
どういうペースでやるのがいいのか、当分は試行錯誤が続きそうです。

問題解答

問題1.14
木の展開はひたすら手を動かす問題。
オーダーについてはよくわからなかったので先人の記録を頼る(ただ、ちょいちょい間違いがあるような気がするので後述のサイトを参考にした方がよさそう)。また、そこから先に更に詳しい解説が。

cc が呼ばれる度に、 1,0 もしくは二つのサブツリーの和を返す。そして、ツリー中の各ノードでは一処理をおこなう。
amount n と kinds-of-coints が 5 の場合、処理の数は (1 + (cc n 4) の処理の数 + (cc (-n 50 5)) の処理の数)。以後、
N(n,m) を subtree (cc n m) によっておこなわれる処理の数とする。

これを踏まえて N について考える。

N(n,　5) = 1 + N(n,　4) + N((n-50), 5)
N(n,　4) = 1 + N(n,　3) + N((n-25), 4)
N(n,　3) = 1 + N(n,　2) + N((n-10), 3)
N(n,　2) = 1 + N(n,　1) + N((n-5), 2)
N(n,　1) = 1 + N(n,　0) + N((n-1), 1)

N(n,　0) = 1
N((n-1),　1) = N(n, 0) + N((n-2), 1)

(n　-　n) まで展開される。

これを後ろからのぼっていく。

N((n-n, 1) = N(0, 1) = 1
N((n-(n-1), 1) = N(1, 1) = 1 + N(1, 0) + N(0, 1) = 1 + 1 + 1 = 3
N((n-(n-2), 1) = N(2, 1) = 1 + N(2, 0) + N(1, 1) = 1 + 1 + 3 = 5
N((n-(n-3), 1) = N(3, 1) = 1 + N(3, 0) + N(2, 1) = 1 + 1 + 5 = 7

つまり、 N(n, 1) = 2n+1

よって、オーダーは Θ(n)。

同様にして N(n, 2) について考える。

N(n, 2) = 1 + N(n, 1) + N((n-5) 2)
= 1 + N(n, 1) + N(n, 1) + N((n-10) 2)

n-5x <= 0 になるまで展開されるので、ざっくり n/5 回繰り返される。 そして、N(n, 1) は Θ(n) なので Θ(n*n/5) = Θ(n^2)。 同様に N(n, 3), N(n, 4), N(n, 5) を考えると N(n, 5) のオーダーは Θ(n^5)。 うはー。一問、2行に二時間かけてもうた。。。