フィボナッチ数列の計算法（1）

http://blog.livedoor.jp/dankogai/archives/50958771.html

O(1)の計算方法が紹介されているが、roundをつかうのはちょっと気持ち悪い（64ビット整数を求めるのに、80ビットの浮動小数点数の計算では追いつかない？）という人のために、O(logn)で計算する方法を紹介する。
たぶん、上では省略されただけだろう。

まず、次を考える。

$\frac{a + b\sqrt{5}}{2}$

~~ここで、a, bは奇数。~~
この形で得られる集合は乗法について閉じていることが簡単に確認できる。
間違えた。a, bは共に奇数、または偶数。すなわち、差が偶数。

$a_{m+n} - b_{m+n} = \frac{(a_m+b_m)(a_n-b_n)+4b_nb_m}{2}$

右辺の分子の第1項はどちらも偶数で、差は偶数となる。

$c \equiv \frac{1 + \sqrt{5}}{2}$

とする。

$c^n \equiv \frac{a_n + b_n\sqrt{5}}{2}$

これが計算できればよい。

$c^{m + n} = c^mc^n$

より、

$a_{m + n} = \frac{a_ma_n + 5b_mb_n}{2}$
$b_{m + n} = \frac{a_mb_n + a_nb_m}{2}$

となるので、cⁿはO(logn)で計算できる（下のコード参照）。
これで、a_n,b_nを求めて、

$F_n = \frac{1}{\sqrt{5}}(\frac{a_n + b_n\sqrt{5}}{2} - \frac{a_n - b_n\sqrt{5}}{2}) = b_n$



function fib(n) {

    var c = [ 1, 1 ];       // (1 + √5) / 2

    var cn = pow(c, n);     // cのべき乗

    

    return cn[1];

}
// べき乗

function pow(c, n) {

    var x = [ c[0], c[1] ]; // コピー

    var m = calcFigure(n);  // nの2進数での桁数

    for(var i = m - 2; i >= 0; i--) {

        x = mult(x, x);

        if((n >> i) & 1)

            x = mult(x, c);

    }

    

    return x;

}
// 乗算

function mult(x, y) {

    return [ (x[0] * y[0] + 5 * x[1] * y[1]) / 2,

             (x[0] * y[1] +     x[1] * y[0]) / 2 ];

}
// 桁を数える

function calcFigure(n) {

    var m;

    for(m = 1; n >> m != 0; m++)

        ;

    return m;

}