15.7.4 Fibonacci数

Fibonacci関数mpz_fib_uiとmpz_fib2_uiは， F[n] もしくは F[n]とF[n-1]を個別に効率的に計算するためのものです。

小さいnの場合，__gmp_fib_tableテーブルに入っている1リム長分の値が使われます。32-bitリムの場合はF[47]まで，64-bitリムの場合はF[93]までが格納されているものです。利便性を上げるため，この表はF[-1]から始まっています。

このテーブルを超える範囲の値の場合は，2進べき乗アルゴリズムを用いてnの2進表現を高い方から低い方へ，F[n] と F[n-1]をペアで求めます。漸化式は下記の通りです。

F[2k+1] = 4*F[k]^2 - F[k-1]^2 + 2*(-1)^k
F[2k-1] =   F[k]^2 + F[k-1]^2

F[2k] = F[2k+1] - F[2k-1]

各ステップにおいて，kはnの高い方のbビット分にあたります。nの次のビットがゼロの時は， F[2k]とF[2k-1]が使用され，1の時はF[2k+1]とF[2k]が使用されます。このプロセスが繰り返され，nの全てのビットに対して行われます。これらの漸化式は，nの1ビットに対して2回の2乗が実施されます。

最初の２，３のnを単一リムテーブル上で加算を使うだけで扱うことができます。この時にはFibonacci数列の漸化式F[k+1]=F[k]+F[k-1]を使います。しかし，nは10～20個程度まとめて扱う方が高速なので，このためのコードは意味がないと思われます。もしそうしたければ，データテーブルを拡大した方がいいでしょう。

テーブルを使うことで，小さい数に対する計算を大幅にカットできますし，小さい nに対しては高速にできます。テーブルを大きくしてやれば更に小さい nに対しては高速にできますが，サイズと計算速度のバランスをどうとるかが問題になってきます。GMPではなるべくコードはコンパクトにしたいので，べき乗アルゴリズムを基本的には使うようにしています。

mpz_fib2_ui関数はF[n]とF[n-1]の両方を返しますが，mpz_fib_ui関数はF[n]だけを計算します。この場合はこのアルゴリズムの最後のステップで，2回の2乗の代わりに乗算を1回だけ実行するようにしています。下記の２式のうち一つがそれで，nの偶奇によって使い分けます。

F[2k]   = F[k]*(F[k]+2F[k-1])

F[2k+1] = (2F[k]+F[k-1])*(2F[k]-F[k-1]) + 2*(-1)^k

ここで，F[2k+1]は上記の場合は同じもので，乗算のためにアレンジしているだけです。 2*(-1)^kの項は両方，桁の借り上げ・桁上がりが発生しないので，低いリムの方向に計算を進めることができるようになっており，コード量の削減ができています。どうやっているかは， mpz_fib_ui関数や基盤関数mpn_fib2_ui関数のコードのコメントをご覧下さい。