対数ステップでフィボナッチ数を計算する

30分プログラム、その581。対数ステップでフィボナッチ数を計算できるらしいので、やってみる。
SICPの演習問題で担当になったところに、フィボナッチ数を対数ステップで計算するアルゴリズムが載っていたので試してみた。

とりあえず、配布資料を貼って置きますね。

対数ステップでフィボナッチ数を計算するアルゴリズムが存在する。

fib-iterではaa := a + b、b := aという変形を行なっていた。
この変形をTと呼ぶと、Tはn回適用されている。
つまり、(1,0)にT^nを適用すると、Fib(n+1)とFib(n)が得られる。

ここで、Tを

  • a := bq + aq + ap
  • b := bp + aq

というTpqの特殊な場合(p=0,q=1)と考える。

Tpqを二回適用したのと同じになるTp'q'を考え、プロシージャを完成させなさい。

で、ごにょごきょ計算すると、p'=p^2+q^2、q'=q^2+2pqということが分かる。あとは、これを使って関数を作ってやればいい。

使い方

*Main> fib 10
89

-- 確認のため1〜10までを計算してみる
*Main> map fib [1..10]
[1,2,3,5,8,13,21,34,55,89]

ソースコード

isEven :: Int -> Bool
isEven n = n `mod` 2 == 0

t p q (a,b) = (b*q + a*q + a*p,b*p + a*q)

fib' p q 0 (a,b) = (a,b)
fib' p q n (a,b) =
    if isEven n then
        fib' ((p*p) + (q*q)) ((q*q) + 2 * p * q) (n `div` 2) (a,b)
    else
        t p q $ fib' p q (n-1) (a,b)

fib n = fst $ fib' 0 1 n (1,0)