整数の公式でフィボナッチ数列を求める

次のコードを用いると、なんとフィボナッチ数列が生成できます。

def fib(n):
    return (4 << n*(3+n)) // ((4 << 2*n) - (2 << n) - 1) & ((2 << n) - 1)

この記事では、その導き方と振る舞いを説明しましょう。

具体的な説明に入る前に、背景としてフィボナッチ数列の概要と計算方法を駆け足で紹介します。すでに数学の専門知識がある方は、導入部分はほとんど飛ばして、「母関数」のセクションをざっと読んでから、「整数の公式」に進んでいただいて構いません。

概要

フィボナッチ数列とは、言わずと知れた以下の数列です。

$$1, 1, 2, 3, 5, 8, 13, 21, 34, 55, \ldots$$

この数列の$n$番目の値は、その直前の2つの値の和です。公式で表現すると、以下のような漸化式になります。

ここでは数列の添え字を、一般的な1ではなく、0から始めています。

フィボナッチ数列の計算

この数列には、わりと有名な計算方法が何種類かあります。以下は明示的な再帰的呼び出しを行う方法で、この実装だと処理速度は遅くなります。

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

次の反復実装の場合、計算量は$O(n)$になります。

    def fib_iter(n):
        a, b = 1, 1
        for _ in xrange(n):
            a, b = a + b, a
        return b

上記よりも少し知名度が下がるのが、以下の行列のべき乗を使用した実装で、計算量は$O(\mathrm{log}\ n)$になります。

    def fib_matpow(n):
        m = numpy.matrix('1 1 ; 1 0') ** n
        return m.item(0)

最後の方法は、fib_iter内のaとbを数列として見なす方法で、以下のように表記できます。

これから以下の式が得られます。

よって、もしなら、$b_n = m_{11}$になります（Pythonと違って、行列の添え字は通常1が基準になることに注意してください）。

NumPy行列のべき乗が繰り返し二乗法のような振る舞いをすると想定すると、計算量は$O(\mathrm{log}\ n)$になります。

さらに、漸化式を解くために、閉じた式を見つける方法もあります。これにより、次の実数値の公式が導かれます：$\phi = (1 + \sqrt{5}) / 2$>、$\psi = (1 - \sqrt{5}) / 2$とすると、$\mathrm{Fib}(n) = (\phi^{n+1} - \psi^{n+1}) / \sqrt{5})$。この手法には、任意精度の実数計算を要するという実用上の欠点がありますが、$n$の値が小さければ問題はありません。

    def fib_phi(n):
        phi = (1 + math.sqrt(5)) / 2.0
        psi = (1 - math.sqrt(5)) / 2.0
        return int((phi ** (n+1) - psi ** (n+1)) / math.sqrt(5))

母関数

任意の数列$a_n$の母関数は、$\Sigma_n a_nx^n$の無限和です。フィボナッチ数列の場合は$\Sigma_n \mathrm{Fib}(n)x^n$になります。つまり、これは無限に続くべき級数であり、$x^n$の係数は$n$番目のフィボナッチ数に相当します。

フィボナッチ数には以下の関係があります。

$$\mathrm{Fib}(n+2) = \mathrm{Fib}(n+1) + \mathrm{Fib}(n)$$

この式に$x^{n+2}$をかけて$n$全体で和をとると、以下の式が得られます。

$$\Sigma_n\mathrm{Fib}(n+2)x^{n+2} = \Sigma_n\mathrm{Fib}(n+1)x^{n+2} + \Sigma_n\mathrm{Fib}(n)x^{n+2}$$

$F(x)$を$Fib$の母関数として、それを$\Sigma_n\mathrm{Fib}(n)x^n$と定義すると、上の式は次のように簡略化できます。
$$F(x) - x - 1 = x(F(x) - 1) + x^2F(x)$$

これは、さらに整理できます。

$$F(x) = xF(x) + x^2F(x) + 1$$

これを$F$について解くと以下の式が得られます。

$$F(x) = \frac{1}{1 - x - x^2}$$

驚くことに、全てのフィボナッチ数列を網羅する簡単でシンプルな公式を、なんとか導くことができました。この式の利用方法については、次のセクションで詳しく説明しましょう。

今はひとまず技術的なことは置いておき、いくつかの$x$の値について$F$を評価してみましょう。そのためには、べき級数が収束する必要があります。フィボナッチ数列は$\phi^n$のように増大して、$|a|$ < 1なら等比級数$\Sigma_n a^n$は収束します。よって、$|x|$ < $1/\phi \simeq 0.618$から、べき級数は収束することが分かります。

整数の公式

ここまでで、Pythonのコードを理解するための準備は整いました。

まずは、この公式を直感的にとらえるため、$10^{-3}$で母関数$F$を評価してみましょう。

興味深いことに、小数展開した部分に$1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89$と、フィボナッチ数列が現れています。魔法のような結果に驚いてしまいますが、その理由は次の式から分かります。
$F(10^{-3}) = \mathrm{Fib}(0) + \mathrm{Fib}(1)/10^3 + \mathrm{Fib}(2)/10^6 + \mathrm{Fib}(3)/10^9 + \ldots$

この例では、フィボナッチ数列が次々と$1/1000$倍されて並んでいきます。つまり、その値が一旦1000を超えると、隣り合う数に影響を及ぼし始めるということです。この現象は、上記の$F(10^{-3})$の計算で988から確認できます。正しいフィボナッチ数は987ですが、数列の次の数から1だけオーバーフローが発生しています。その結果Off-by-oneエラーが発生し、以降はパターンが崩れてしまうのです。

しかし、いかなる$n$の値に対しても、10の負の指数を十分大きく取れば、たとえオーバーフローが発生したとしても、$n$番目のフィボナッチ数に悪影響が出ることはありません。ここでは、ある$k$という値について$10^{-k}$が妥当な値になると仮定しましょう。この値は後ほど選定します。

さらに整数で計算をしたいので（その方がコーディングしやすいので）、全体を$10^{kn}$倍して$n$番目のフィボナッチ数が整数の範囲にくるようにして、式を整理します。

$$10^{kn} F(10^{-k}) = \frac{10^{kn}}{1 - 10^{-k} - 10^{-2k}} = \frac{10^{kn+2k}}{10^{2k} - 10^{k} - 1}$$

この結果を$10^k$を法として見ると、$n$番目のフィボナッチ数が得られます（先ほども書きましたが、$k$には十分大きな値を選んだものと想定しています）。

では先に進む前に、底を10から2に変換しましょう。理由は単にプログラミングを楽にするためです。

$$2^{kn} F(2^{-k}) = \frac{2^{k(n+2)}}{2^{2k} - 2^{k} - 1}$$

あとは$\mathrm{Fib}(n+1)$ < $2^k$になるように、$k$を十分大きく取るだけです。フィボナッチ数列は$\phi^n$のように増大して、$\phi$ < $2$なので、$k = n+1$とすれば安全です。

これらを統合すると以下の式が得られます。

Pythonでは左シフトの表記法が使えるので、$a$ << $k = a \cdot 2^k$の場合は次のように書けます。

$\mathrm{mod}\ (2$ << $n)$の部分は、$(2$ << $n) - 1$ とのビット単位の論理和（&）として表せるので、元のPythonプログラムは次のように書き直せます。

  def fib(n):
        return (4 << n*(3+n)) // ((4 << 2*n) - (2 << n) - 1) & ((2 << n) - 1)

非反復型で閉じた解が得られたのは興味深いですが、これは全く実用的な手法ではありません。ここではサイズが$O(n^2)$ビットの整数を用いて、整数演算を実行しています。でも実際には、最終的にビット単位で論理和を取る前に、最初の$n$個のフィボナッチ数が全て連結した整数値を取得しているのです。