昨日追い切れなかった問題 1.41 からの続き。

問題解答

問1.41
結局、自力で完全に展開するのは無理だったので答えを写経しながら展開。それでも必ずどこかで間違える。

まあ、大体の流れは把握できたので OK とする。

休日やら差し込み仕事やらで気づいたら間が結構空いてしまいました。気を取り直して再開。

本文

1.3.4 値として返される手続き
不動点の探索、平均緩和法、y→x/y を使って平方根や立方根を求めるコードを導き出している。

Newton 法
1.1.7 の例と比べると格段にシンプルな記述になっている。

抽象と第一級手続き

問題解答

問1.40

(define (cubic a b c)
  (newtons-method (lambda (x)
                    (+ (* x x x) (* a x x) (* b x) c))
                   1.0))

こんな感じ?

問1.41

(define (double f)
  (display f)(newline)
  (lambda (x) (f (f x))))

(define (inc x)
  (display x)(newline)
  (+ x 1))

(((double (double double)) inc) 5)

実行すると 16 回 inc が呼ばれるが、8 回しか呼ばれないはずなのになぜ?
→解答を見ながら手で展開してみるも、挫折。明日再チャレンジしてみる。

問題解答

問1.37

(define (cont-frac n d k)
  (define (cont-sum i)
    (if (= i k)
        (/ (n i) (d i))
        (+ (/ (n i) (+ (d i) (cont-sum (+ i 1)))))))

  (cont-sum 1))

再帰版はうまく作れなかったのでカンニング。

(define (cont-frac-iter n d k)
  (define (cf-iter v i)
    (if (= i 0)
        v
        (cf-iter (/ (n i) (+ (d i) v)) (- i 1))))

  (cf-iter (/ (n k) (d k)) k))

4 桁の精度の近似を得る。

(cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 9)
→0.6181818181818182
(cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 10)
→0.6179775280898876
(cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 11)
→0.6180555555555556
(cont-frac (lambda (i) 1.0) (lambda (i) 1.0) 12)
→0.6180257510729613

11 くらいあれば OK?

問1.38

(+ 2 (cont-frac (lambda (i) 1.0)
           (lambda (i)
             (if (= (mod (+ i 1) 3) 0)
                 (* (/ (+ i 1) 3) 2.0)
                 1.0))
           12))

問1.39

(define (tan-cf x k)
  (define (cont-tan i)
    (let ((squareX (* x x))
          (d (- (* i 2) 1)))
      (cond ((= i k) (/ squareX d))
            ((= i 1) (/ x (- d (cont-tan (+ i 1)))))
            (else (/ squareX (- d (cont-tan (+ i 1))))))))

  (cont-tan 1))

(tan-cf 1.0 10)

何だか、他の SICP 読書ブログを見てると 4,5 日でこのあたりまで来てる人達ばっかりで、しかもちゃんと問題も解いてるっぽいんですが自分の進みが恐ろしく遅いんですかね。。。。自分は一ヵ月でここなんですが(白目

本文

1.3.3 一般的方法としての手続き
区間二分法をテーマにした内容。一気にプログラムっぽい感じのコードが出てきた。個人的にはこちらの方が慣れてる感じがするのでとっつきやすい。

関数の不動点の探索

不動点 : x が f(x) = x を満たす時、x を関数 f の不動点という。ある初期値に対して f を繰り返し作用させることで求めることができる。

この不動点を求める方法を使って収束計算をおこなうことができる。ただし、これだけでは値が振動してしまう場合もある。それを緩和するために予測値との平均をとる方法があり、これを平均緩和法という。

問題解答

問1.35
φ^2=φ+1 より、両辺をφで割って φ=1+1/φ。

(define (goldenRatio x)
  (fixed-point (lambda (y) (+ 1.0 (/ 1.0 y)))
               1.0))

(goldenRatio 1.0)

問1.36

(define (fixed-point f first-guess)
  (define (close-enough? v1 v2)
    (< (abs (- v1 v2)) tolerance))
  (define (try guess)
    (display guess)
    (newline)
    (let ((next (f guess)))
      (if (close-enough? guess next)
          next
          (try next))))
  (try first-guess))

平均緩和未使用

(fixed-point (lambda (x) (/ (log 1000) (log x))) 10.0)

gosh> 10.0
2.9999999999999996
6.2877098228681545
3.7570797902002955
5.218748919675316
4.1807977460633134
4.828902657081293
4.386936895811029
4.671722808746095
4.481109436117821
4.605567315585735
4.522955348093164
4.577201597629606
4.541325786357399
4.564940905198754
4.549347961475409
4.5596228442307565
4.552843114094703
4.55731263660315
4.554364381825887
4.556308401465587
4.555026226620339
4.55587174038325
4.555314115211184
4.555681847896976
4.555439330395129
4.555599264136406
4.555493789937456
4.555563347820309
4.555517475527901
4.555547727376273
4.555527776815261
4.555540933824255
4.555532257016376

平均緩和使用

(fixed-point (lambda (x) (/ (+ x (/ (log 1000) (log x))) 2)) 10.0)
gosh> 10.0
6.5
5.095215099176933
4.668760681281611
4.57585730576714
4.559030116711325
4.55613168520593
4.555637206157649
4.55555298754564
4.555538647701617
4.555536206185039

本文

1.3.2 lambda を使う手続きの構築
ここでついに lambda が登場。

非 lambda 版だと

(define (pi-sum a b)
  (define (pi-term a)
    (/ 1.0 (* a (+ a 2))))

  (define (pi-next n)
    (+ n 4))

  (sum pi-term a pi-next b))

こうなるのが、 lambda 版だと

(define (pi-sum a b)
  (sum (lambda (x) (/ 1.0 (* x (+ x 2))))
       a
       (lambda (x) (+ x 4))
       b))

こうなるよと。

局所変数を作り出す let の使い方

ここでやっと局所変数が使えるように。この説明が無いためにどれだけ大変だったか。。。。

解説の、lambda 版と let 版のコードの構造がわかり辛かった(カッコ多すぎやねん!!)ので自分がわかりやすいように整形してみた。

lambda
(define (f x y)
  (
    (lambda (a b)
      (+ (* x (square a))
         (* y b)
         (* a b))
      )
      (+ 1 (* x y))
      (- 1 y)
    )
  )
)

let
(define (f x y)
  (let
    (
      (a (+ 1 (* x y)))
      (b (- 1 y))
    )

    (+ (* x (square a))
       (* y b)
       (* a b)
    )
  )
)

問題解答

問1.34
(f f)
(f 2)
(2 2)

となってエラーになる。
何故ここでこの問題をするのか?という意図がわからない。

どうにも、SICP 内だけで説明されている Scheme の構文だけだといろいろ辛い気がするんですが、それでもまずはこの範囲内でやりくりした方がいいんだろうか?というのがそこはかとない疑問。

問題解答

問1.31
再帰版

(define (product term a next b)
  (define (iter a result)
    (if (> a b)
        result
        (iter (next a) (* result (term a)))))

  (iter a １))

反復版

(define (product term a next b)
  (if (> a b)
      1
      (*  (term a)
          (product term (next a) next b))))

(define (pi-next i) (+ i 2))
(define (pi-factor n) (/ (* (- n 1) (+ n 1)) (* n n)))

(* (product pi-factor 3.0 pi-next 1000) 4.0)

数列の作り方が思い浮かばなかった。。。。これはプログラミング力以前の算数力の低下が。。。

問1.32

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

(define (product term a next b)
  (accumulate * 1 term a next b))

(* (product pi-factor 3.0 pi-next 1000) 4.0)

問1.33

(define (filtered-accumulate filter combiner null-value term a next b)
  (define (iter a result)
    (cond ((> a b) result)
          ((filter a) (iter (next a) (combiner result (term a))))
          (else  (iter (next a) result))))

  (iter a null-value))

a.

(define (next i) (+ i 1)
(define (square x) (* x x))

(define (prime-squaresum a b)
  (define (prime? n)
    ...)

  (filtered-accumulate prime? + 0 square a next b))

問題回答

問1.29
先日から引き続きコードを見直し。おかしいところは見つけたもののそれでもなおらず。ギブアップ。
回答を見ながらコードを修正。

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

関数だけではなく、変数もこのように定義できるのか。

そして、コードを直していくうちにそもそも出題されている内容から離れてしまっていた模様。
回答を見ながら自分のコードを修正して、正しい結果が得られるようになった。

(define (sum term a next b)
  (if (> a b)
      0
      (+  (term a)
          (sum term (next a) next b))))


(define (cube x)
  (* x x x))


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

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

  (define (term k)
    (* (cond ((odd? k) 4)
          ((or (= k 0) (= k n)) 1)
          ((even? k) 2))
       (y k)))

  (* (/ h 3.0) (sum term 0 next n)))

問1.30

(define (sum term a next b)
  (define (iter a result)
    (if (> a b)
        result
        (iter (next a) (+ result (term a)))))

  (iter a 0))

この問題はサラッと。

問題回答

問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。
たしかにだまされている。