OCaml/SICP-1b

練習問題の概要と解答案 (1章)

正解かどうかの保証はしない。問題の解釈を間違えている可能性もある。

1.11

再帰関数。 値を再利用するから memoize (メモ化) したほうがよいのだろうが そんなことはここではまだ問われていないはず。

let rec fr n =
  if n<3 then n else
    fr (n-1) + 2 * fr (n-2) + 3 * fr (n-3)

iteration的 (場合分けは後ろでやっても良いか)。

let fi n =
  let rec fii n1 n2 n3 n =
    if n=0 then n3
    else if n=1 then n2
    else if n=2 then n1
    else fii (n1 + 2*n2 + 3*n3) n1 n2 (n-1) in
    fii 2 1 0 n

で、本当の繰り返し :-)

let fl n =
  let n1 = ref 2 in
  let n2 = ref 1 in
  let n3 = ref 0 in
  let f  = ref 2 in
  for i = 3 to n do
    f :=  (!n1 + 2* !n2 + 3* !n3) ;
    n3 := !n2;
    n2 := !n1;
    n1 := !f;
  done;
  if n<3 then n else ! f;;

iteration 的に再帰するのってわかりにくいけど for loop なら直感的。

1.12

パスカルの三角形。

let rec pasctr n r =
  if r=0 then 1 
  else if n=1 then 1 
  else if n=r then 1
  else (pasctr (n-1) (r-1)) + (pasctr (n-1) r)

1.13

この問題はプログラムを書くのではなく、数学。

三項間漸化式の解き方なんか忘れてしまったなぁ。がんばってやってみよう。

f(n) = f(n-1) + f(n-2) の特性方程式 x^2 = x + 1 を解くと、x = (1+-√5)/2 (=Φ,Ψ)。

したがって、 f(n)-Ψf(n-1) = Φ(f(n-1) - Ψf(n-2)) であるから、 f(n)-Ψf(n-1) = Φ^(n-2) (f(2)-Ψf(1)) = Φ^(n-2)(1-Ψ)=Φ^(n-1)。

(追記: 双対から f(n)-Φf(n-1)=Ψ^(n-1) で連立方程式を解く方が早い。以下は遠回り)

g(n) = f(n)/Φ^n とすると、Φg(n)-Ψg(n-1)=1。特性方程式は x = 1/(Φ-Ψ) = 1/√5

ゆえに Φ(g(n) - 1/√5) = Ψ(g(n-1) - 1/√5) 。 g(n)-1/√5 = (Ψ/Φ)^n (g(0)-1/√5) = (Ψ/Φ)^n (-1/√5)

f(n) = Φ^n/√5 + Ψ^n*(-1/√5) = (Φ^n - Ψ^n)/√5。

これでフィボナッチ数列の一般式が証明できた。

ここで -0.6 > Ψ > -0.7 なので |Ψ^n /√5| < 0.5 が成り立つので Q.E.D.

1.14

OCaml では前方参照(後方参照？)出来ないので記述の順番に注意。

let first_denomi kind =
  if kind = 1 then 1
  else if kind = 2 then 5
  else if kind = 3 then 10
  else if kind = 4 then 25
  else if kind = 5 then 50
  else -1 (* 本来は例外を上げるべき *)

let rec cc amount kind =
  if amount = 0 then 1
  else if (amount < 0) || (kind = 0) then 0
  else (cc amount (kind - 1)) + (cc (amount - (first_denomi kind)) kind)

let count_change amount = cc amount 5

これでツリーを書いてみると、

count_change <-- 11
cc 11 5
cc 11 4 + cc -39 5 
 |        0
cc 11 3 + cc -14 4
 |        0
cc 11 2 + cc 1 3
 |        cc 1 2 + cc -9 3
 |         |        0
 |        cc 1 1 + cc -4 2
 |         |        0
 |        cc 1 0 + cc 0 1
 |         0        1
cc 11 1 + cc 6 2
 |        cc 6 1 + cc 1 2
 |         |       cc 1 1 + cc -4 2
 |         |        |        0
 |         |       cc 1 0 + cc 0 1
 |         |        0        1
 |        cc 6 0 + cc 5 1
 |         0       cc 5 0 + cc 4 1
 |                  0       cc 4 0 + cc 3 1
 |                           0       cc 3 0 + cc 2 1
 |                                    0       cc 2 0 + cc 1 1
 |                                             0       cc 1 0 + cc 0 1
 |                                                      0        1
cc 11 0 + cc 10 1
 0        cc 10 0 + cc 9 1
               ...
                    cc 1 0 + cc 0 1
                     0        1

count_change --> 4
- : int = 4

枝わかれ多いな。オーダーどれくらいなんだ？最低でも線形オーダだが、べきオーダに見えるな。

cc a b = cc a b-1 + cc a-k5 b 
= ... + cc a-k1 b-4 + cc a-k2 b-3 + cc a-k3 b-2 + cc a-k4 b-1 + cc a-k5 b

O(cc a b) = ΣO(cc a-δ i) i=1..5 

b=5 で決め打つと a->a+1 で a の定数倍だけ増えるように見える (b が小さい奴は b が大きい奴よりも極限的にはオーダーは同程度以下)。つまり、これを積分することで、計算量は自乗オーダ。であってるかなぁ？

メモ化すれば速くなりそうではある。

1.15

とりあえず OCaml 版...って手を加えてるけど。

let rec sine angle =
  let abs x = if x>0.0 then x else -. x in
  let sine3 x = sine (x /. 3.0) in
  let reduc f = f *. 3.0 -. 4.0 *. f *. f *. f in 
    if (abs angle) < 0.0000001 then angle 
    else reduc (sine3 angle) ;;

a. 5回かな。

b. 計算量は log オーダ。同じだけスタックを食うので領域も log。

1.16

イテレータ型で冪乗計算しろ、と言っている。 まずは再帰型。

let rec bekir a x =
  let modulo = x mod 2 in
  let hbekir a x = bekir a (x/2) in 
  let hbekir2 a x = (hbekir a x) * (hbekir a x) in
    if x=1 then a 
    else if modulo = 1 then a * (hbekir2 a x)
    else (hbekir2 a x)

再帰型が頭にあるとイテレータは難しい。問題のヒントと、下位ビットから処理していくことを考慮する。 (再帰型でもスタックの奥底で下位ビットから戻ってくるわけだからこの部分を念頭に置くとよさげ)

let bekii a x =
  let rec b_iter a x s =
    if x=0 then s
    else if (x mod 2) = 1 then b_iter a (x-1) (s*a)
    else b_iter (a*a) (x/2) (s)
  in
    b_iter a x 1

あるいはこっちの方が普通か？

let bekii2 a x =
  let rec b_iter a x s =
    if x=0 then s
    else if (x mod 2) = 1 then b_iter (a*a) (x/2) (s*a)
    else b_iter (a*a) (x/2) s
  in b_iter a x 1

てか命令型で書いた方が早いよ。

let bekil a x =
  let n = ref x in
  let res = ref 1 in
  let bs = ref a in
  while !n > 0 do
    if (!n mod 2) = 1
    then (n := !n - 1; res := !res * !bs)
    else (n := !n / 2; bs := !bs * !bs)
  done;
    ! res;;

let bekil2 a x =
  let n = ref x in
  let res = ref 1 in
  let pow = ref a in
  while !n > 0 do
    if (!n mod 2) = 1
    then (res := !res * !pow;)
    else ();
    n := !n/2;
    pow := !pow * !pow 
  done;
    ! res;;

教訓：イテレータ型は命令型で書いてからコンバートする方が早い。

1.17

冪乗のやりかたと同じように掛け算も足し算の繰り返しで定義できる。 1ずつ足すバージョンは例のとおりだから、log オーダの高速版を作れ。

書き直すだけだからチョロい。

let multi a b =
  let rec m_iter a b p =
    if b=0 then p
    else if (b mod 2) = 1 then m_iter (a*2) (b/2) (p+a)
    else m_iter (a*2) (b/2) p in
    m_iter a b 0

1.18

1.17 をイテレータで...って、 1.17 は recursive で書くことを要請してたのか？ じゃ逆にこっちを recursive 版で。

let rec multr a x =
  let modulo = x mod 2 in
  let hmultr a x = bekir a (x/2) in 
  let hmultr2 a x = (hmultr a x) * 2 in
    if x=1 then a 
    else if modulo = 1 then a + (hmult2 a x)
    else (hmult2 a x)

動作確認してない。これはもういいや。

1.19

フィボナッチ数列をlogオーダで計算する方法。フィボナッチの作用素Tとすると T^n(1,0) でf(n)を得る。T = Tpq in p=0,q=1 とした一般作用素 Tpq を考える。定義を以下と置くとき、

Tpq(a,b): a <- bq+aq+ap, b <- bp+aq

(Tpq)^2 = Tp'q'としたときの p' q' の関係式を求めて、それを元にフィボナッチ関数を作れ。

問題文が長いが、実質数学の問題っていうか数式計算。

T^2pq(a,b) はどうなるかというと、

a <- a((p+q)^2 + q^2) + b(2pq+q^2)
b <- b(p^2 + q^2) + a(2pq+q^2)

なので、T^2pq = Tp'q' となる p',q' は

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

従って空き部分はこの式を突っ込めばおしまい。

1.20

OCaml は applicative-order 動作をする、といういい方でいいのか？

let rem a b = a mod b
let rec gcd a b =
  if b=0 then a else gcd b (rem a b);;

# #trace gcd;;
gcd is now traced.
# #trace rem;;
rem is now traced.
# gcd 206 40;;
gcd <-- 206
gcd --> <fun>
gcd* <-- 40
rem <-- 206
rem --> <fun>
rem* <-- 40
rem* --> 6
gcd <-- 40
gcd --> <fun>
gcd* <-- 6
rem <-- 40
rem --> <fun>
rem* <-- 6
rem* --> 4
gcd <-- 6
gcd --> <fun>
gcd* <-- 4
rem <-- 6
rem --> <fun>
rem* <-- 4
rem* --> 2
gcd <-- 4
gcd --> <fun>
gcd* <-- 2
rem <-- 4
rem --> <fun>
rem* <-- 2
rem* --> 0
gcd <-- 2
gcd --> <fun>
gcd* <-- 0
gcd* --> 2
gcd* --> 2
gcd* --> 2
gcd* --> 2
gcd* --> 2
- : int = 2