« OpenGLは右手座標系? | メイン | 10キーで棋譜入力するiモードアプリ »

 フィボナッチ数列の任意の項を高速に計算する方法

そんなプログラミングの問題を見かけて、反射的に思ったのは、一般項を用いることだった。
フィボナッチ数列の漸化式はan=an-1+an-2なので、一般項anを求めることは可能である。

3項を含む漸化式は、特性方程式を用いて解くのが基本であるが、やり方を忘れてしまったので、自己流で解いてみた。
a_{n+1} = a_{n} + a_{n-1}

a_{n+1} + Aa_n = (1+A)(a_n+ Aa_{n-1})
という形にできると、
(1+A)(a_n+ Aa_{n-1}) = (1+A)^2 (a_{n-1}+ Aa_{n-2}) = \cdots = (1+A)^n(a_1+Aa_0)
であることと、a1=1, a0=0であることを用いて
a_{n+1} + Aa_n = (1+A)^n
とできる。このAは、an-1の係数を見ると、(1+A)A=1となるAなので、A=(-1±√5)/2である。従って、

\left\{
\begin{array}{c}
a_{n+1} + \frac{-1+\sqrt{5}}{2} a_n = \left(1 + \frac{-1+\sqrt{5}}{2}\right)^n \\
a_{n+1} + \frac{-1-\sqrt{5}}{2} a_n = \left(1 + \frac{-1-\sqrt{5}}{2}\right)^n
\end{array}
\right.
であり、上の式から下の式を引いて整理すると、
a_n = \frac{1}{\sqrt{5}} \left( \left(\frac{1+\sqrt{5}}{2}\right)^n - \left(\frac{1-\sqrt{5}}{2}\right)^n \right) と求まる。
念の為申し上げるが、このanはnが何であろうとも整数である。実にミラクルである。

√5のn乗なんてのがあるのでは、誤差が出まくるだろうとは思ったが、CやJavaでは浮動小数点の精度が低くて無理でも、PythonのmathパッケージにMathematicaやMaxima並の驚異の性能が無いとは限らないので、試しに

#!/usr/bin/python
from math import pow
from math import sqrt

def fib(n):
    return int(1/sqrt(5)*(pow((1+sqrt(5))/2,n)-pow((1-sqrt(5))/2,n)))

an_2 = 0  #a_{n-2}
an_1 = 0  #a_{n-1}
an = 1    #a_{n}
for i in range(1, 100):
    print "n =", i, "a_n =", an, "fib(n) =", fib(i)
    if fib(i) != an:
        print "error!"
    an_2 = an_1
    an_1 = an
    an = an_1 + an_2
とやってみた所、72項目で誤差が出た。
Maximaはn>200000でも正確な整数を弾き出す(下記入出力参照、下3桁が正しいことを確認)が、さすがにPythonでは無理だった。おそらく、Maximaは内部で√5を√5として保存しながら数式処理をしてるからであろう。Pythonの言語仕様上、pow()の結果は一度float値にしないといけないので、正確な整数値を出すには無限の精度が必要になってしまう。

Maximaへの入力とその出力
(%i1) 1/sqrt(5)*(((1+sqrt(5))/2)^n-((1-sqrt(5))/2)^n),n=234567,ratsimp;
(%o1) 179607932890115544888825400217[48962 digits]014673023583001136817746498978

MathematicaやMaximaを使って良いなら、一般項を用いる方法もあるが、今回の問題は、「複数言語選択可能」としながら、MathematicaもMaximaも選択肢に無かった。 フィボナッチ数列のn項目を高速に計算する方法は、他には思い付かなかったのだが、少し調べると、計算量がO(log N)になる簡単な方法があった。

M=\pmatrix{1 & 1 \cr 1 & 0} とすると、Mn-1の{1,1}要素がフィボナッチ数列のn項目だというのが使える。
さらに、M=M*Mという計算を繰り返すことにより、M^(2k)は計算量少なく求められるので、Mn-1をM^(2k)の積とすることによって、高速に計算できる。
例えば、
M5 = 0*M8 · 1*M4 · 0*M2 · 1*M
M10 = 1*M8 · 0*M4 · 1*M2 · 0*M
であるように、2進数にした時の1のビットに対応するM^(2k)を掛け合わせれば良いのである。

従って、例えば次のようにすれば、n=10000でも高速に計算される。(Pythonで扱える整数の上限はメモリが許す限り無限である)

#!/usr/bin/python

class multipliable_2x2matrix(object):
    def __init__(self, a, b, c, d):
        self.a = a;
        self.b = b;
        self.c = c;
        self.d = d;
    def mul(self, b):
        return multipliable_2x2matrix(
            self.a * b.a + self.b * b.c,
            self.a * b.b + self.b * b.d,
            self.c * b.a + self.d * b.c,
            self.c * b.b + self.d * b.d)
    
def fib(n):
    if n == 0: return 0
    m = multipliable_2x2matrix(1, 1, 1, 0)
    f = multipliable_2x2matrix(1, 0, 0, 1)
    n = n - 1
    while n > 0:
        if n & 1 == 1:
            f = f.mul(m)
        n = n >> 1
        m = m.mul(m)
    return f.a

トラックバック

このエントリーのトラックバックURL:
http://ynomura.dip.jp/cgi-bin/mt/mt-tb.cgi/227

コメント投稿フォーム

※投稿されたコメントはオーナーが承認するまで表示されません。


Powered by Movable Type 3.35