はじめに

今回は、フィボナッチ数列¹の任意の項を得る関数について調査・考察してみました。

そもそも「フィボナッチ数列」とは？

以下の漸化式で表される数列です。

\[\eqalign{ f(0) &= 0 \\ f(1) &= 1 \\ f(n+2) &= f(n+1) + f(n) \\ }\]

ここでは一般項を “0番目” から順に 0, 1, 1, 2, 3, 5, 8, 13, 21, 34, … と考えます。この数列の第 n 項の値を得るプログラムを考えていきたいと思います。

フィボナッチ数列の第n項を得るプログラムの例

(1) シンプルな再帰呼び出し

再帰呼び出しを用いる、最もシンプルな例を書いてみました。

#
# 01_simple.rb
#
def fibonacci_simple( n )
  case n
  when 0, 1
    return n
  else
    return fibonacci_simple( n - 1 ) + fibonacci_simple( n - 2 )
  end
end

if $0 == __FILE__
  p fibonacci_simple( ARGV[0].to_i )
end

このプログラムでは、

f(5) -> f(4) + f(3) -> (f(3) + f(2)) + (f(2) + f(1)) -> ....

と次々に元の関数が呼ばれるため、n = 20 でも fibonacci_simple の呼び出し回数は 20000 回以上！になります。

この時の計算量²は $O(2^n)$ です。

$ time ruby -r profile 01_simple.rb 20
6765
  %   cumulative   self              self     total
 time   seconds   seconds    calls  ms/call  ms/call  name
 80.00     0.68      0.68    21891     0.03     0.51  Object#fibonacci_simple
(略)

real    0m1.086s

$ time ruby -r profile 01_simple.rb 25
75025
  %   cumulative   self              self     total
 time   seconds   seconds    calls  ms/call  ms/call  name
 84.08     9.98      9.98   242785     0.04     0.82  Object#fibonacci_simple

real    0m12.996s

本記事では時間計測とコール回数測定のため、 “time ruby -r profile XXXX.rb” というように time コマンド使用、プロファイルモードで実行します。

n = 20 では 1 秒程度で計算できていますが、n = 25 ではコール回数、実行時間ともに n = 20 のときの 10 倍以上になっています。このように文字通り指数関数的に実行時間が増えていくため、 n の値がごく小さい場合しか使えない方法ともいえます。

(2) 動的計画法を利用する

プログラミングで有用なアルゴリズムの一つとして「動的計画法」³があります。これは「対象となる問題を複数の部分問題に分割し、部分問題の計算結果を記録しながら解いていく手法」の総称で、ざっくり言うと「同じ計算は何度もせずに記録しておいて使う方法」です。

動的計画法を用いたプログラムの実現方法として、トップダウン方式とボトムアップ方式があります。ここではそれぞれについてプログラムを作成します。

(2a) メモ化再帰の利用（トップダウン）

まずは「メモ化再帰」を使ったコードです。こちらは @memo という Array に結果を格納していきます。

すでに @memo[k] が存在する場合は保存されている値をそのまま返し、存在しない場合は再帰呼び出しで @memo[k-1] と @memo[k-2] を参照します。以下 @memo[k-1] に対しても同様にして参照していくため、これはトップダウンの方法といえます。

#
# 02a_memoize.rb
#
@memo = [ 0, 1 ]
def fibonacci_memoize( n )
  @memo[n] ||= fibonacci_memoize( n - 1 ) + fibonacci_memoize( n - 2 )
end

if $0 == __FILE__
  p fibonacci_memoize( ARGV[0].to_i )
end

$ time ruby -r profile 02a_memoize.rb 20
6765
  %   cumulative   self              self     total
 time   seconds   seconds    calls  ms/call  ms/call  name
100.00     0.01      0.01       39     0.26     1.03  Object#fibonacci_memoize
(略)
real    0m0.075s

$ time ruby -r profile 02a_memoize.rb 100
354224848179261915075
(略)
real    0m0.130s

$ time ruby -r profile 02a_memoize.rb 9990

(1) の方法と比べると n = 20 で fibonacci_memoize のコール回数が 39 回と劇的に小さいコール回数になっています。 また n = 100 でも高速に計算できていることがわかります。

計算量は$O(n)$です。ただ、大きな n (私の環境だと n = 10000 前後) になるとコールスタックが深くなりすぎて “stack level too depp” エラーで落ちて計算できませんでした。

(2b) ループの利用（ボトムアップ）

続いてループを利用する方法です。 n1 と n2 を順次使いまわしながらボトムアップで値を求めます。

#
# 02b_iterate.rb
#
def fibonacci_iterate( n )
  n1 = 0; n2 = 1

  (n-1).times{
    n1, n2 = [ n2, n1 + n2 ]
  }

  n == 0 ? 0 : n2
end

if $0 == __FILE__
  p fibonacci_iterate( ARGV[0].to_i )
end

$ time ruby -r profile 02b_iterate.rb 25
75025

real    0m0.081s

$ time ruby -r profile 02b_iterate.rb 100000
(略)

real    0m4.785s

コールが深くなりすぎるということはなく、100000 という値でもしっかり計算できています。

こちらも計算量は$O(n)$ですが、プログラム自体はシンプルなものやメモ化再帰を用いたものより若干読みにくい（コードの意図が把握しづらい）実装かな、というのが率直な感想です。

(3) 一般項の公式を利用する

次に紹介するのは、フィボナッチ数#一般項 に記載のある Binet の公式を使って実装したものです。

\[F_{n}={\frac {1}{\sqrt {5}}}\left\{\left({\frac {1+{\sqrt {5}}}{2}}\right)^{n}-\left({\frac {1-{\sqrt {5}}}{2}}\right)^{n}\right\}\]

#
# 03_binet.rb
#
def fibonacci_general( n )
  r5 = Math.sqrt(5)
  ((((1 + r5) * 0.5) ** n - ((1 - r5) * 0.5) ** n) / r5).round
end

if $0 == __FILE__
  p fibonacci_general( ARGV[0].to_i )
end

$ time ruby -r profile 03_general.rb 10
55
(略)
real    0m0.073s

$ time ruby -r profile 03_general.rb 100
354224848179263111168 # 誤り！！
(略)
real    0m0.074s

やや複雑な計算式ですが、n = 10、 n = 100 ともに動的計画法と同等の時間で計算できているようです。

計算量は $O(\log n)$ になります。詳細は省略しますが、これはべき乗の計算が \(x^{100} = x^{64}\cdot x^{32}\cdot x^4\) のように表せることを利用して、$O(n)$ ではない回数の演算で実装されているからです。

が、致命的なことに n = 100 の場合の値が間違っています（正しくは 354224848179261915075 ）。さらにn > 1000 くらいになると Infinity (FloatDomainError) が発生し計算できなくなります。

これは浮動小数点を使用しているために生じる現象ですね。式変形で $\sqrt{5}$ を消し整数演算に持ち込む方法 ⁴ もあるみたいですが、そうすると結局 $O(n)$ になってしまったりするようで、あまりメリットはないようです。

(4) 行列のべき乗を利用する

最後に行列のべき乗を利用する方法を紹介します。

#
# 04_matrix.rb
#
require 'matrix'

def fibonacci_matrix( n )
  base   = Matrix[[1, 1], [1, 0]]
  result = base ** n
  result[0,1]
end

if $0 == __FILE__
  p fibonacci_matrix( ARGV[0].to_i )
end

こちらも Binet の公式と同じように Wikipedia に式があります。

\[\biggl( \matrix{ F_{n+1} & F_n \\ F_n & F_{n-1} \\ } \biggr) = \biggl( \matrix{ 1 & 1 \\ 1 & 0 \\ } \biggr)^n\]

行列のべき乗に関しても \(B^{100} = B^{64}\cdot B^{32}\cdot B^4\) のように表せることを利用した高速な計算方法があり、これにより計算量は $O(\log n)$ になります。

$ time ruby -r profile 04_matrix.rb 100
354224848179261915075 # しっかり正解の値になってます
(略)
real    0m0.111s

$ time ruby -r profile 04_matrix.rb 100000
(略：20899桁数になりました)
real    0m0.158s

$ time ruby -r profile 02b_iterate.rb 100000
# (2b)との比較。(4)が圧倒的に速いことがわかる。
(略)
real    0m4.031s

(2b) の方法はこれまで紹介した中では唯一、現実的な時間で大きな値 ( n = 100000 ) に対応できるものでしたが、これと比較しても、行列のべき乗を利用した方法はかなり速いということがわかります。

まとめ

以下の 5 つの方法でフィボナッチ数列の第 n 項を求めるプログラムを書いてみました。

• (1) シンプルな再帰呼び出し
• (2a) 動的計画法：メモ化再帰の利用（トップダウン）
• (2b) 動的計画法：ループの利用（ボトムアップ）
• (3) 一般項の公式の利用（Binetの公式）
• (4) 行列のべき乗の利用

主観的ですが、以下に各手法の性格をまとめてみました。

	(1)	(2a)	(2b)	(3)	(4)
安定性	○	△ (Stack level too deep)	○	△ (Float Domain Error)	○
精度	○	○	○	△	○
計算速度	×	○	○	◎	◎
考え方・手法の汎用性	○	○	◎	△	△
コードの理解しやすさ	◎	◎	○	△	△

私が日ごろプログラミングをするときは「性能Max」よりも「できるだけシンプルに」「コードは理解しやすく」を意識しているので、同じようなケースがあった場合は今回のようにまずは (1) で実装してみて厳しければ (2a)or(2b)、さらに厳しければ (3)or(4) のような順番で採用を検討していってもいいのかなと思いました。シンプルな方法で性能がでなければ、メモ化やループを検討する。それでもダメなら数式や根本的なアルゴリズムに踏み込んで、、というように。

アルゴリズムは大事ですね。

参考

1. フィボナッチ数 - Wikipedia ↩

2. ランダウの記号 - Wikipedia ↩

3. 動的計画法 - Wikipedia ↩

4. フィボナッチ数列の第k項をO(log n)で計算しようとしたけど難しかった話 - shogonir blog ↩