フィボナッチのロジック最適化

フィボナッチに慣れてきました。
ここまで来ると最速を知りたくなります。
最近の最速はどうやら固定行列を使った繰り返し自乗法の極まったやつなようです。
バイブルはこちら。

今日はこれを読み解いていきます。

Textbook recursive(教科書的な再帰)

よく書かれているコードですね。再帰で計算するやつで絶望的に遅く、スタックオーバーフローしちゃうやつです。黄金比を使った計算も出来ますが、誤差が発生するため実用には向きません。

import timeit

def fibonacci(n):
    if n < 2:
        return n
    else:
        return fibonacci(n-2) + fibonacci(n-1)

def measure(n):
    t = timeit.timeit(f'fibonacci({n})', globals=globals(), number=1)
    return f'fibonacci({n}): {t}'

print(fibonacci(10))
# 55
print(measure(35))
# fibonacci(35): 3.0515678379997553

Dynamic programming(動的プログラミング)

再帰計算をforなどを使ってループさせるやつです。私が書いてたのもコレ。ロジック的には何の工夫もないので遅いです。

import timeit

def fibonacci(n):
    if n < 2:
        return n
    v0, v1 = 0, 1
    for i in range(n-1):
        v0, v1 = v1, v0 + v1
    return v1

def measure(n):
    t = timeit.timeit(f'fibonacci({n})', globals=globals(), number=1)
    return f'fibonacci({n}): {t}'

print(fibonacci(10))
# 55
print(measure(35))
# fibonacci(35): 8.58999919728376e-06
print(measure(500000))
# fibonacci(500000): 2.5509368060011184

とはいえ再帰に比べると6桁くらいは高速ですが、遅いです。

Matrix exponentiation(行列のべき乗)

ここからが本題になります。

フィナボッチは以下のFを満たします。

\begin{pmatrix}
1 & 1 \\
1 & 0 
\end{pmatrix}
^n
=
\begin{pmatrix}
F(n+1) & F(n) \\
F(n) & F(n-1)
\end{pmatrix}

証明については上のバイブルに載っているので省略。
この式の左辺を使って計算するのですが、ここで繰り返し自乗法を使います。原理はコチラ。

自乗の自乗の…を計算することで計算を端折れます。

import timeit

def fibonacci(n):
    if n < 2:
        return n
    x = [1, 1, 1, 0]
    y = [1, 0, 0, 1]
    n -= 1
    while n > 0:
        if n % 2:
            mul(x, y, y)
        mul(x, x, x)
        n //= 2
    return y[0]

def mul(x, y, z):
    z[0], z[1], z[2], z[3] =\
        x[0]*y[0]+x[1]*y[2],\
        x[0]*y[1]+x[1]*y[3],\
        x[2]*y[0]+x[3]*y[2],\
        x[2]*y[1]+x[3]*y[3]

def measure(n):
    t = timeit.timeit(f'fibonacci({n})',  globals=globals(), number=1)
    return f'fibonacci({n}): {t}'

print(fibonacci(10))
# 55
print(measure(500000))
# fibonacci(500000): 0.14098625300175627
print(measure(3000000))
# fibonacci(3000000): 2.76616032200036

行列計算は重いわけですが、掛け算を端折れる分Nが大きくなるとすこぶる速いです。N=500000でもう1つ前のやつの10倍以上高速になってます。

Fast doubling(高速ダブル)

ついに本稿最速です。
行列のべき乗で使用した式を書き下して以下を抽出したようです。証明はバイブルを参照してください。

\begin{eqnarray}
F(2k+1) = F(k+1)^2 + F(k)^2 \\
F(2k) = F(k)(2F(k+1) - F(k)^2)
\end{eqnarray}

これを使うと、$F(2k+1)$か$F(2k)$を$F(k)$と$F(k+1)$で計算可能になります。つまり、元の定義式を使うと$N$=$2k+1$か$2k$を引数にして、$F(N)$,$F(N+1)$を返す再帰関数が実装できます。すると$F(1)=1$,$F(0)=0$に当たるまで再帰させれば計算できることになります。事実上行列のべき乗から極限まで肉を削いだ計算になってるように思います(計算量のオーダーは同じ)。

これを実装した綺麗なコードがバイブルにpublic domainで載っているのでそのまま利用させてもらうと…

import timeit

def fibonacci(n):
    if n < 0:
        raise ValueError("Negative arguments not implemented")
    return _fib(n)[0]

def _fib(n):
    if n == 0:
        return (0, 1)
    else:
        a, b = _fib(n // 2)
        c = a * (b * 2 - a)
        d = a * a + b * b
        if n % 2 == 0:
            return (c, d)
        else:
            return (d, c + d)

def measure(n):
    t = timeit.timeit(f'fibonacci({n})',  globals=globals(), number=1)
    return f'fibonacci({n}): {t}'

print(fibonacci(10))
# 55
print(measure(3000000))
# fibonacci(3000000): 0.36685713999759173
print(measure(10000000))
# fibonacci(10000000): 2.3166784589993767

さらに桁が違う速度になりました。天晴。

最後にまとめというか感想

最初たったN=35個で数秒かかってた処理が、同程度の時間でN=10000000まで出来るようになるのは感無量。行列のべき乗までは検索するとチラホラ出てきましたが、その高速化ロジック(最速)は日本語だと見かけなかったので、そこは良かったかなと思います。

ロジックが簡単で重い処理というのはベンチマーク向きなので、最適化されたロジックを使うと少々味気ない感じもします。ただ言語比較までするようなベンチマークでは、あまりに重いと何を測ってるのか分からない状態になるので、最低限のパフォーマンス考慮はした方がいいのかもしれません。