フィボナッチ数 - 名古屋で数学するプログラマ(仮)

今日は、みんな大好きフィボナッチ数列。
柏餅おいしかったとかちまき上品な甘さだったとか実家のタケノコごはんがとかそういうGW後半ネタはすぱっと省略。

フィボナッチ数列。
「1, 1,」から始まって、それ以降は「前の2つの数の和」が続く数列、のこと。数式（漸化式）による定義は以下の通り（0番目を0とした定義）。

$F_0 = 0, F_1 = 1, F_n = F_{n-1} + F_{n-2}(n\geq2)$

最初の方を列挙すると以下の通り。

(0, )1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, …

できるだけ短いコードで列挙する、という所謂CodeGolfのネタにもなっています。
私の最高記録はRubyで27バイト…このあたりまでは誰でも行き着くんですよね、あと4バイトどうやって短くするんだろう(+_+)

でも今日は「短いコード」の話ではなく、一般項の話。つまり、
「整数nを与えたときに、フィボナッチ数列のn番目の数（＝n番目のフィボナッチ数）を返す関数」
を実装する話。
それをアルゴリズムの話だけでなく、もうちょっと数学のスパイスを加えた話をしてみます。

手続き指向の解答例

プログラムを少しでも組んだことのある人なら、きっと誰でも思いつくのは、「定義に従って順番に計算して答えを出す方法」。
Rubyで書くと、こんな感じ。

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0  # for negative n
  a, b = 0, 1
  n.times { a = b + b = a }
  a
end

関数*1内1行目は、n<0だった場合の対応。今回はあまり気にしないでください。本質はここじゃないので。
さて、もちろんこれは正しいフィボナッチ数を返しますが、効率の方はどうでしょう？
ベンチマークを取ってみました。
n=50,100,500,1000,5000,10000 の場合の各計算を1000回ずつ実施してその実行時間を計測（Rubyの benchmark ライブラリを使用）

                      user     system      total        real
fib(50)x1000      0.010000   0.000000   0.010000 (  0.006343)
fib(100)x1000     0.010000   0.000000   0.010000 (  0.011449)
fib(500)x1000     0.200000   0.000000   0.200000 (  0.198729)
fib(1000)x1000    0.480000   0.000000   0.480000 (  0.475711)
fib(5000)x1000    3.460000   0.160000   3.620000 (  3.633362)
fib(10000)x1000   8.080000   0.860000   8.940000 (  8.946042)

うん、ほぼ当然のごとく、nにだいたい比例した数値が出ていますね。つまり計算量はO(n)です。
十分小さなnに対しては問題ありませんが、大きくなるに従ってどんどん時間がかかります。

関数指向の解答例

続いて、少しプログラムやアルゴリズムを勉強して「再帰関数」を覚えた人、もしくは「名古屋人なら関数型言語」な人なら、こんなコードを思いつくでしょう。

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0  # for negative n
  return n if n < 2
  fib(n - 2) + fib(n - 1)
end

つまり「定義をそのまま実装」した例です。
これを実行してみると、もちろん正しい答えは出るのですが…。

              user     system      total        real
fib(10)   0.000000   0.000000   0.000000 (  0.000037)
fib(20)   0.000000   0.000000   0.000000 (  0.004549)
fib(30)   0.400000   0.000000   0.400000 (  0.399396)
fib(40)  48.120000   0.040000  48.160000 ( 48.225261)
# fib(50) takes about 1:40'
# fib(60) takes about 8.6days…

はい。非常に時間がかかります。上のベンチマーク、先ほどと違って1000回繰り返しているのではなく、1回だけしか実行していない処理時間であることに注意。
n=50で約1時間40分、n=60で約8日と14時間かかります。n=100,500,1000,… なんて、恐ろしくて実行できません。計算量はO(eⁿ)。
これは関数が、同じ引数を与えると毎回同じ処理をするから。

関数指向＋メモ化

同じ引数が与えられたら同じ結果を返すようにするために、前に計算した値を覚えておくという方法があります。これを「メモ化」と言います。
先ほどの例に、メモ化の要素を加えた例を示します。

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0  # for negative n
  fib_hash = Hash.new {|h, n|
    h[n] = case n
      when 0, 1
        n
      else
        h[n - 2] + h[n - 1]
    end
  }
  fib_hash[n]
end

Ruby特有の書き方になっているので関数型言語のメモ化とはちょっと違うし分かりにくい実装ですが、やってることは同じ。
「一度計算した（他の）フィボナッチ数は、ハッシュに記憶しておいて、後で再利用できるようにする」コードになっています。
これだけで、実行時間は劇的に短くなります。、が…。

                     user     system      total        real
fib(50)x1000     0.050000   0.000000   0.050000 (  0.050571)
fib(100)x1000    0.100000   0.000000   0.100000 (  0.100457)
fib(500)x1000    0.740000   0.000000   0.740000 (  0.742134)
fib(1000)x1000   1.760000   0.010000   1.770000 (  1.774478)
#fib(5000)x1000  #@> stack level too deep (SystemStackError)

n=1000 を超えたあたりになると、再帰呼び出しの深さが深くなりすぎて、コンピュータが処理できずにエラーになったりします。
それに、計算量そのものは同じO(n)なのですが、ベンチマークの結果を見てみると、最初の「手続き指向」よりもわずかながら時間がかかっています。
つまり関数型指向の実装は、必ずしも効率が良い方法ではない、ということなのでしょうか？
そうとも言えますし、実はそうとばかりは言えない部分もあったりします。

メモ化＋倍数公式

さぁ、ここからがこの記事のメインです。

フィボナッチ数は、その定義となる漸化式から、色々な公式が導き出せます。
そして余り知られていないのですが、2n、2n+1番目のフィボナッチ数について、こんな面白い公式があります。

$F_{2n}=(2F_{n-1}+F_n)F_n\\F_{2n+1}=F_n^2+F_{n+1}^2$

これをF_n=F_n-1+F_n-2の代わりに使ってみましょう。

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0  # for negative n
  fib_hash = Hash.new {|h, n|
    h[n] = case n
      when 0, 1
        n
      else
        n.even? ? (h[n/2-1]*2 + h[n/2]) * h[n/2] : h[n/2]**2 + h[n/2+1]**2
    end
  }
  fib_hash[n]
end

ベンチマークは以下の通り。

                      user     system      total        real
fib(50)x1000      0.020000   0.000000   0.020000 (  0.019347)
fib(100)x1000     0.020000   0.000000   0.020000 (  0.019735)
fib(500)x1000     0.030000   0.000000   0.030000 (  0.032989)
fib(1000)x1000    0.050000   0.000000   0.050000 (  0.043373)
fib(5000)x1000    0.080000   0.000000   0.080000 (  0.079305)
fib(10000)x1000   0.130000   0.000000   0.130000 (  0.135395)

n=10000でもSystemStackErrorが発生せず、しかもnが100を超えたあたりからは手続き指向の解答例よりも明らかに高速になっています。n=10000だと、「メモ化＋倍数公式」は「手続き指向」の約60倍速！
これは、今までの方法は結局「n番目以下の全てのフィボナッチ数も算出」していたのに対して、今回の方法が「計算に必要なフィボナッチ数だけを算出」しているから、というのが最も大きな理由。
算出する個数は、最初にn/2番目あたりを探し、次にそのさらにn/2あたりを探し…ということで、つまりlog₂n に大体比例することになります。つまり計算量はO(logn)。かなり理想的です。
ただ、nが十分に小さい場合は2乗や掛け算の分だけ余分に計算することになるので、実時間は手続き指向の場合に敵わないこともある、というわけです。

おまけ：フィボナッチ数（＝フィボナッチ数列の一般項）の公式

n番目のフィボナッチ数は、nのみに依る式で以下のように表せます。

$F_n=\frac{1}{\sqrt 5}\left\{\left(\frac{1+\sqrt{5}}{2}\right)^n-\left(\frac{1-\sqrt{5}}{2}\right)^n\right\}$

証明はめんどくさいしブログ1ページじゃHeavyなのでしませんが、興味のある方は、「数学ガール」（ISBN: 4797341378）でも読んでみてください。

これをうまく利用して、「途中のフィボナッチ数を一切計算せず直接算出する」方法ももちろん考えられます。
ただ、そのまま実数（浮動小数点数）を用いてやろうとすると誤差の問題が、それを無理矢理整数範囲で収めようとすると計算量の問題が、結局ついて回ります。
世の中うまくいかないものですね。

まとめ

計算量の考え方は大事。
あと、アルゴリズムの工夫も大事だけど、数学的な裏付けとか、「数学的な考察による数式変形（公式導出）によるさらなる効率化」が、とっても大事。
ちなみに「倍数公式」という呼び方は正式な用語かどうか分かりませんが、同様にF_3nやF_4nに対する公式も存在し（というか探せばいくらでも作り出すことができ）ますが、それらを利用してプログラムを書き換えてみても、必ずしもこれ以上効率的な結果にはなりません。
なぜでしょう？興味があればご自分で考察してみてください。

*1:Rubyの場合正式には「メソッド」と言いますが、言語はここでは重要ではないので一般的名称として「関数」という用語で通します。