フィボナッチ数を計算する方法をまとめてみた

Project Eulerで巨大な数のフィボナッチ数(10^14番目ぐらい)の剰余を求める必要があったので、関連するアルゴリズムについてまとめました。行列を使った最適化については@nico_shindanninさんにヒントをいただきました。

定義通り求める

すぐに限界のくる方法。一度の関数呼び出しで2度再帰することと、同じフィボナッチ数を何度も計算するので時間がかかる。

def fib_naive(n):
    if n in [0,1]:
        return 1
    else:
        return fib_naive(n-1) + fib_naive(n-2)

メモ化して求める

ナイーブな方法より計算が減るので速くなるけど大きなフィボナッチ数を求めようとするとメモリをバカ喰いする。

def fib_memo(n, m={0:0,1:1}):
    if not n in m:
        m[n] = fib_memo(n-1, m) + fib_memo(n-2, m)
    return m[n]

一般項から求める

一番速そうだけれど、浮動小数点を使うので大きなnで求めるにはDecimalモジュールでもの凄く大きな桁を使うか、整数化して計算するとか必要で面倒。

def fib_general(n):
    sqrt5 = 5**.5
    g = (1 + sqrt5) / 2
    fn = (g**n - (-g)**-n) / sqrt5
    return int(fn)

足し算で求める

割と速くて、再帰も無いので10万程度なら簡単に計算できる。剰余を計算するときも計算途中の桁数を抑えたまま進められるのもいい。

def fibonacci_sequence():
    a,b = 0,1
    while True:
        yield a
        a,b = b,a+b #b,(a+b)%m
def fib_iter(n):
    def take(g, m):
        for i in xrange(m):
            yield next(g)
    for fnum in take(fibonacci_sequence(), n):
        fn = fnum
    return fn

行列のべき乗で求める

行列M={{1,0},{0,1}}に対してM'={{1,1},{1,0}}をn回かけ算した結果のM[0,1]とM[1,0]がFnになるアルゴリズム。M^2n = M^n * M^nのように、一度計算した結果を使い回して指数をどんどん上げられるので速い。x^yを再帰的に計算していく方法に良く似ている感じ。

x^yの計算

def power(x, y):
    return [lambda:power(x, y/2)**2, lambda: 1][y<2]() * [1,x][y%2]
M=[]
def _fib_matrix(n,m):
    global M
    if n > 1:
        _fib_matrix(n/2,m)
        M = [
            [M[0][0]*M[0][0] + M[0][1]*M[1][0],
             M[0][0]*M[0][1] + M[0][1]*M[1][1],
            ],
            [M[1][0]*M[0][0] + M[1][1]*M[1][0],
             M[1][0]*M[0][1] + M[1][1]*M[1][1],
            ],
        ]
        if m:
            M = [map(lambda x:x%m, v) for v in M]
    if n % 2 == 1:
        N= [[1,1], [1,0]]
        M = [
            [M[0][0]*N[0][0] + M[0][1]*N[1][0],
             M[0][0]*N[0][1] + M[0][1]*N[1][1],
            ],
            [M[1][0]*N[0][0] + M[1][1]*N[1][0],
             M[1][0]*N[0][1] + M[1][1]*N[1][1],
            ],
        ]
        if m:
            M = [map(lambda x:x%m, v) for v in M]

def fib_matrix(n, modulo=None):
    global M
    M = [[1,0], [0,1]]
    _fib_matrix(n, modulo)
    return M[0][1]

てきとうな比較

それぞれの関数でどこまで計算できるか試してみました。予めsysモジュールで再帰の最大深さを上げてあります。

import sys; sys.setrecursionlimit(1000000)

環境はCPython2.6.5です。

関数 n 計算時間 備考
fib_naive 32 2.14
fib_memo 20000 0.059 このあたりから上を実行するとSERVでPythonが落ちる
fib_general 1420 0.000 このあたりから上を実行するとOverflowError
fib_iter 20000 2.368
fib_matrix 1000000 2.150

参考
http://www.ics.uci.edu/~eppstein/161/960109.html