SICP 問題 1.19 (逐次平方変換を利用したFibonacci数を求める末尾再帰処理)

これ結構楽しかった。


【問題】
Fibonacci数を対数的ステップ数で計算するうまいアルゴリズムがある。
1.2.2節の fib-iter プロセスで使った状態変数 a, b の変換、


a ← a + b
b ← a
に注意しよう。この変換を T と呼ぶ。
1 と 0 から始め、T を繰り返して n 回作用させると、Fib(n+1) と Fib(n) の対ができる。
言い換えればFibonacci数は、対 (1, 0) に T^n、つまり変換 T の n 乗を作用させると得られる。
さて、T_pq は対 (a, b) を、

a ← bq + aq + ap
b ← bp + aq
にしたがって変換するものとし、T を変換族 T_pq の p=0, q=1 の特例としよう。

変換 T_pq を2回使うとその効果は同じ形式の変換 T_p'q' を1回使ったのと同じになることを示し、p', q' を、p, q を使って表せ。
これで変換を二乗する方法が得られ、fast-exptのように逐次平方を使い、T^n を計算することができる。
これらをすべてまとめ、対数的ステップ数の以下の手続きを完成せよ。

(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'を計算
                   <??>  ;q'を計算
                   (/ count 2)))
        (else
         (fib_iter (+ (* b q) (* a q) (* a p))
                   (+ (* q a) (* p b))
                   p
                   q
                   (- count 1)))))


【解答】

問題文には求める計算式のほぼ全容が記載されており、p' と q' の穴埋め問題になっていた。
が、もったいないのでその部分は見ずに、関数を生成するところまで含めて全部自前で解いてみた。
なので、以下の解答の記述は一部問題文の関数と表記が異なる箇所があるがご了承を。

今回の柱は、

1.(T_pq)^2 と T_p'q' の変換結果から、p', q' を求める。
2.逐次平方による末尾再帰のFibonacci数を算出する関数を求める。

の二本立て。



■(T_pq)^2 と T_p'q' の変換結果から、p', q' を求める。

(1)まずは「(T_pq)^2」から。

T_pq を1回作用させたものを a' b'とすると、問題文の変換式より次のようになる。


a' = bq + aq + ap
b' = bp + aq
次に、a', b' にもう1回 T_pq を作用させたものを a'' b'' とすると、

a'' = b'q + a'q + a'p
b'' = b'p + a'q
と表記できる。
これらから a'', b'', a 及び b の関係を求めると、次のようになる。


a'' = b'q + a'q + a'p
= (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p

b'' = b'p + a'q
= (bp + aq)p + (bq + aq + ap)q


(2)次に「T_p'q'」を求める。

T_p'q' を1回作用させたものは、a'' b'' になり、かつ T_pq とは同じ変換方式になるので、


a'' = bq' + aq' + ap'
b'' = bp' + aq'


(3)「(T_pq)^2 = T_p'q'」から、a, b それぞれの係数を比較。
(1)、(2)の結果それぞれを、「X・a + Y・b」の形式に変換することにより、
a, b それぞれの係数を比較することで p', q', p, q の関係式が導き出される。

では(1)から。


a'' = b'q + a'q + a'p
= (bp + aq)q + (bq + aq + ap)q + (bq + aq + ap)p
= bpq + aq^2 + bq^2 + aq^2 + apq + bpq + apq + ap^2
※a, b で項をまとめる。
= (2q^2 + 2pq + p^2)・a + (2pq + q^2)・b      ・・・(i)

b'' = b'p + a'q
= (bp + aq)p + (bq + aq + ap)q
= bp^2 + apq + bq^2 + aq^2 + apq
※a, b で項をまとめる。
= (2pq + q^2)・a + (p^2 + q^2)・b          ・・・(ii)

次に(2)を変換。

a'' = bq' + aq' + ap'
= (p' + q')・a + q'・b               ・・・(iii)

b'' = bp' + aq'
= q'・a + p'・b                   ・・・(iv)

(i)〜(iv)までのそれぞれの係数を比べてみると、つぎのようになる。

p' + q' = 2q^2 + 2pq + p^2 ・・・(i)と(iii)の a の係数の比較
q' = 2pq + q^2 ・・・(i)と(iii)の b の係数の比較
q' = 2pq + q^2 ・・・(ii)と(iv)の a の係数の比較
p' = p^2 + q^2 ・・・(ii)と(iv)の b の係数の比較
関係式が過剰供給気味だが、上記から、次の関係式が得られる。

p' = p^2 + q^2
q' = q^2 + 2pq
これが1本目の柱の答えね。


■逐次平方による末尾再帰のFibonacci数を算出する関数を求める。

さて、再帰処理の登場人物を考えよう。


a : 1つ前の項の数値(と同時に解答でもある)
b : 2つ前の項の数値
p : 変換係数その1
q : 変換係数その2
i : あと何回計算するべきか。
こんなもんか?登場人物が出揃ったところで関数のおおまかな全体像を記述してみよう。

(define (Fib_log a b p q i)
  (cond ((= i 0) b)
        ((even? i) )
        (else )))

こんな感じになる。こいつの導出はもういいよね?メインディッシュは、★と☆の導出の仕方だし。
で、その?と?にはそれぞれこんな処理をさせたい。


★ : T_p'q'に相当する変換係数 p', q' を求める処理を行ってから、
同じ a, b でFib_logを呼び出す。
その際、i を半分にしてしまう。
☆ : T_pqに相当する変換処理を行ってからFib_logを呼び出す。
その際、i を 1 減算する。
具体的に実装してみよう。


★ = (Fib_log a
b
(+ (* p p) (* q q))
(+ (* q q) (* 2 p q))
(/ i 2))

☆ = (Fib_log (+ (* b q) (* a q) (* a p))
(+ (* q a) (* p b))
p
q
(- i 1))

おお!これでいけるっしょ!?
完成形を書いてみるぜ。

(define (Fib_log a b p q i)
  (cond ((= i 0) b)
        ((even? i)
         (Fib_log a
                  b
                  (+ (* p p) (* q q))
                  (+ (* q q) (* 2 p q))
                  (/ i 2)))
        (else
         (Fib_log (+ (* b q) (* a q) (* a p))
                  (+ (* q a) (* p b))
                  p
                  q
                  (- i 1)))))

実行してみる。


gosh> (map (lambda (x) (Fib_log 1 0 0 1 x)) '(0 1 2 3 4 5 6 7 8 9 10))
(0 1 1 2 3 5 8 13 21 34 55)
gosh>
おお〜。今回はかなりすんなりできたなぁ。