フィボナッチ数(2) ver.末尾再帰 - 名古屋で数学するプログラマ(仮)

暑中残暑お見舞い申し上げます。
この記事書き始めた頃は暑中だったのに色々あって立秋過ぎちゃった。

さておき、前回予告した、「フィボナッチ数を求めるための末尾再帰」について。

おさらい …フィボナッチ数

詳細は、前回の記事とか、適当にググってみたりとかしてください。

簡単に言うと、

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

で定義される数列で、例えば

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

こんな関数（メソッド）を書くと一般の n 番目のフィボナッチ数を求められます。
ちなみに最初の方を列挙すると以下の通り。

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

末尾再帰について少しだけ

一般に関数の「再帰呼び出し」というのは、関数の実装の中で自分自身を呼び出す手法のこと。
「末尾再帰」というのは、その再帰呼び出しを関数の最後のステップで（のみ）行う手法のこと。

例えば、再帰呼び出しを利用した階乗の定義を考えると分かりやすいですね。
↓

def factorial n
  return 1 if n < 2
  n * factorial(n-1)
end

最後に自分自身を呼び出しています。逆にそれ以外の部分では余分な再帰呼び出しは行っていません。
これも末尾再帰と考えることも出来ますが、さらに推し進めて「最後の行は本当に再帰呼び出ししか行わない」ようにしてみましょう。この場合↓こうなります。

def factorial n, r=1
  return r if n < 2
  r *= n
  n -= 1
  factorial n, r
end

本当に最後の行が、元の関数の定義式とほぼ同じになるように、少し冗長に書きました。とにかく最後の行は、パラメータを少し変えて自分自身を再帰呼び出ししている、ただそれだけのことをしているのが分かると思います。
これが厳密な「末尾再帰」のカタチ。

この形になっていると、再帰呼び出しをジャンプに変換することで実行効率を上げる、所謂「末尾最適化」が可能になる*1、という特徴もあるのですが、今回はその話はパス。本題はそっちじゃないので。

フィボナッチ数を末尾再帰で(1) O(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

この例も、再帰呼び出しになっています。しかし最後のステップは、自分自身を2回呼び出しています。
これは末尾再帰のカタチではありません。

これをなんとか末尾再帰のカタチに持って行きたいと思います。どうすれば良いでしょう？

ここで、前節の階乗の例で、2つの実装がどう違っていたのかを思い出してみましょう。

最初の例の最後のステップは、詳しく見るとこうなっています。

再帰呼出しで、nが一つ少ない場合の階乗を実際に計算している。
それを計算した上で、階乗の定義に合わせて今回のnの場合の答えを算出している。

この仕組み、先ほどのフィボナッチ数の関数と同じ構造になっていますよね。

再帰呼出しで、nが1少ない場合と2少ない場合のフィボナッチ数を実際に計算している。
それを計算した上で、フィボナッチ数の定義に合わせて今回のnの場合の答えを算出している。

それに対して、階乗の2つめの例は、最後のステップがこんな構造になっています。

最後のステップに行く前に、カウンタとなる変数と、実際の値を計算するパラメータを変更している。
新しいパラメータで再帰呼出している。

この「実際の値を計算するパラメータ」が重要で、階乗の例の場合、実際に階乗をステップごとに計算しています。
ただし、先の例は「1に、1〜nをこの順に掛けている」形になっているのに対して、後の例は「1に、n〜1をこの順に掛けている」形になっています。
でもこれは大きな問題ではありません。この場合、結果はどちらも同じになりますから。

それを踏まえて。
フィボナッチ数の場合は、どうすれば良いのかというと。
まず、カウンタ以外のパラメータは少なくとも2つ必要です。なぜなら、フィボナッチ数の定義に則った計算を行う以上、少なくとも2つのフィボナッチ数からしか算出できない構造になっているのですから。
そして、その2つの値を使って計算して、その結果をまたパラメータとして次の再帰呼出に与える、ということを繰り返していくことを考えればOK、というわけです。

その2つをこのようにすれば、実はうまくいきます。

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

つまり。
最初はa=0(=F₀),b=1(=F₁)とし、再帰呼出のたびに、aにb、bにa+bを代入することでF_n-2とF_n-1をF_n-1とF_nに変換している、というわけ。
これを所定の回数繰り返せば、aにF_nの計算結果が入っているはずだから、それを返して終わり。

もちろんすぐに分かるように、この関数の場合の計算量は O(n) です。単純にn回繰り返しているだけですからね。

改良(?)バージョン

少し考えると、この方がもっと効率よいんじゃないか？という以下の実装も思いつきます。

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

つまり1回の再帰呼出で、カウンタを1ずつ減らす代わりに2ずつ減らして、与えるパラメータも「F_n-2とF_n-1をF_nとF_n+1に変換」する、というわけ。
計算量はやはり O(n) ですが、実際は、というと。

実は実行速度の面ではむしろ遅くなります。1回繰り返すときに、先の例では足し算1回だけだったのに対して、今回は1回繰り返すとき（先の2回分に相当）に足し算3回。つまり1.5倍になってしまうからです。
ただし再帰呼出の回数は減るので、関数呼び出しのオーバーヘッドとか、スタックオーバーフローの問題は軽減できます。ただこの話も今回の本題ではないので言葉だけ。

フィボナッチ数を末尾再帰で(2) O(logn) バージョン

次は、倍数公式をうまく活用して計算量が O(logn) になるようにできないか考えてみましょう。
「倍数公式」とは、以下のようなモノでした。

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

実はこれをもっと一般化した、以下のような「和の公式」が存在します。

$F_{n+m}=F_nF_m+F_{n-1}F_m+F_nF_{m-1}\\F_{n+m-1}=F_nF_m+F_{n-1}F_{m-1}$

これを利用するとさらにこんな風に書けます：

$F_{2^k+m}=F_{2^k}F_m+F_{2^k-1}F_m+F_{2^k}F_{m-1}$

ここで0≦m＜2^kと考えれば、nを2進数表記して、それを利用してまずF_nをF_2^kとF_m（およびその周辺のフィボナッチ数）で表して、F_mをまた同様に分解して…という繰り返しで計算していくことが出来そうです。

実際にやってみたのが、こんな例になります：

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0
  fibFastIter = lambda do |a, b, p, q, c|
    return a if c < 1
    if c.even?
      a1, b1, p1, q1, c1 = a, b, p**2 + q**2, (2*p+q)*q, c.div(2)
    else
      a1, b1, p1, q1, c1 = a*p + b*q, a*q + b*q + b*p, p, q, c-1
    end
    fibFastIter[a1, b1, p1, q1, c1]
  end
  fibFastIter[0, 1, 0, 1, n]
end

ちょっと先ほどまでと書き方が変わってますが、自分自身を再帰関数としていたか、別の再帰関数を（ラムダ式で）用意したかの違いだけなのであまり気にしないでください。
p, q の方が、倍数公式を利用した、F_2^kを表すフィボナッチ数となっており、
a, b の方が、和の公式を利用した、F_2^k+mを表すフィボナッチ数となっています。
で、繰り返していって最後に、aに求めるフィボナッチ数が入っている、という仕組み。

これ、結構速いです。前回の記事の『メモ化＋倍数公式』に肉薄するくらい。どちらも O(logn) ですしね。
ベンチマークを取ってみました。

                      user     system      total        real
fib(50)x1000      0.010000   0.000000   0.010000 (  0.008081)
fib(100)x1000     0.010000   0.000000   0.010000 (  0.010353)
fib(500)x1000     0.020000   0.000000   0.020000 (  0.025588)
fib(1000)x1000    0.040000   0.000000   0.040000 (  0.033170)
fib(5000)x1000    0.080000   0.000000   0.080000 (  0.082407)
fib(10000)x1000   0.190000   0.000000   0.190000 (  0.192046)

与える n が十分に小さいときは、メモ化＋倍数公式の時よりも速くなっているのが分かります。
でも n=5000 の時に同じ数値となり、n=10000 のときは、逆に少し時間かかってますね。
それは、この方式にまだ少し無駄があるから。

改良(!)バージョン

実は先のコードは、ググって見付けたHaskelか何かの関数型言語の実装例をRubyに翻訳したモノです。
でも、明らかな無駄があるの、分かりますよね？
カウンタ変数 c を、どんどん2で割れば良いのに、奇数の時は1引いて再帰呼び出しして、偶数の時に2で割っています。これをまとめちゃえば、効率よくなるはず。
ということでまとめちゃいました。

def fib n
  return n.even? ? -fib(-n) : fib(-n) if n < 0
  return 0 if n < 1
  fibFastIter = lambda do |a, b, p, q, c|
    return a*p + b*q if c < 2
    a, b = a*p + b*q, a*q + b*q + b*p if c.odd?
    fibFastIter[a, b, p**2 + q**2, (2*p+q)*q, c.div(2)]
  end
  fibFastIter[0, 1, 0, 1, n]
end

ベンチマーク。

                      user     system      total        real
fib(50)x1000      0.000000   0.000000   0.000000 (  0.005049)
fib(100)x1000     0.010000   0.000000   0.010000 (  0.006286)
fib(500)x1000     0.020000   0.000000   0.020000 (  0.016950)
fib(1000)x1000    0.030000   0.000000   0.030000 (  0.029278)
fib(5000)x1000    0.060000   0.000000   0.060000 (  0.062313)
fib(10000)x1000   0.130000   0.000000   0.130000 (  0.132980)

うん。やっぱり先のモノより速くなりました(^_^) n=10000 の時も『メモ化＋倍数公式』の時と同じ数値を出しています。
ていうか、これ以上大きな n のことを考えると、やっぱりメモ化を利用した場合の方が速くなる傾向にあるってことですよね、これ。
メモ化利用の方が再帰呼出の回数は多いけれど、演算（足し算や掛け算）の総回数は下回っていくから、なんでしょうね。

*1:実際には厳密バージョンではなくても末尾最適化は適用可能です。