projectpn
projectpn

フィボナッチ数の定義

フィボナッチ数は、
$\displaystyle F_{n+2}=F_{n+1}+F_n,F_0=0,F_1=1$
という漸化式で定義され、
0,1,1,2,3,5,8,13,21,34,…
という数列になる。
この数列の値を、コンピュータを用いてなるべく速く求めることを考える。
まずは計算法とその時間計算量を示す。(時間計算量の証明は後でやったほうが楽なので後でまとめて行う)

漸化式に従った計算

+ ...
まずは、漸化式に従って計算することを考える。
$\displaystyle F_n=F_{n-1}+F_{n-2} \\ =(F_{n-2}+F_{n-3})+(F_{n-3}+F_{n-4}) \\ \vdots $
と計算する。時間計算量は、
足し算 全体
{\scriptstyle O(1)} {\scriptstyle O\Bigl(\bigl(\frac{1+\sqrt 5}{2}\bigr)^n\Bigr)}
{\scriptstyle O(\log(n))} {\scriptstyle O\Bigl(\bigl(\frac{1+\sqrt 5}{2}\bigr)^n\Bigr)}以上
となるので、この方法ではnが大きくなると飛躍的に計算量が増大する。

動的計画法

+ ...
計算を高速化する方法の一つに、動的計画法がある。
具体的には、
$\displaystyle F_2=F_1+F_0 \\ F_3=F_2+F_1 \\ F_4=F_3+F_2 \\ \vdots \\ F_n=F_{n-1}+F_{n-2} $
というように、最初から順番に計算してゆく。このとき時間計算量は、
足し算 全体
{\scriptstyle O(1)} {\scriptstyle O(n)}
{\scriptstyle O(\log(n))} {\scriptstyle O(n^2)}
となりかなり改善される。

漸化式の改良

+ ...
動的計画法でかなり高速化したとはいえ、nが非常に大きくなると大変である。そこで、別の形の漸化式を用いることで高速化を図る。
このように漸化式を変形し、
$\displaystyle F_n=F_{n-1}+F_{n-2} \\ =2F_{n-2}+F_{n-3} \\ =3F_{n-3}+2F_{n-4} \\ =5F_{n-4}+3F_{n-5} \\ \vdots \\ =F_{k+1}F_{n-k}+F_{k}F_{n-k-1} $
この式に{\scriptstyle n=2m,k=m}や、{\scriptstyle n=2m+1,k=m}を代入すると、
$\displaystyle F_{2m}=F_{m+1}F_m+F_mF_{m-1} \\ =F_m(F_{m+1}+F_{m-1}) \\ =F_m(F_m+2F_{m-1}) \\ F_{2m+1}=F_{m+1}F_{m+1}+F_mF_m \\ =F_{m+1}^2+F_m^2 $
この漸化式を用いてそのまま計算すると、時間計算量は以下のようになる。
掛け算 足し算 全体
{\scriptstyle O(1)} {\scriptstyle O(1)} {\scriptstyle O(n)}
{\scriptstyle O((\log(n))^2)} {\scriptstyle O(\log(n))} {\scriptstyle O(n^3)}
{\scriptstyle O(\log(n)\log(\log(n)))} {\scriptstyle O(\log(n))} {\scriptstyle O(n^2\log(n))}
これは先程の動的計画法より遅いが、この漸化式にも動的計画法が適用できて、例えば、n=100の時の計算を行おうと思えば、
$\displaystyle F_2=F_1(F_1+2F_0) \\ F_3=F_2^2+F_1^2 \\ F_4=F_2(F_2+2F_1) \\ F_5=F_3^2+F_2^2 \\ F_6=F_3(F_3+F_2) \\ F_7=F_4^2+F_3^2 \\ F_{12}=F_6(F_6+2F_5) \\ F_{13}=F_7^2+F_6^2 \\ F_{24}=F_{12}(F_{12}+2F_{11}) \\ F_{25}=F_{13}^2+F_{12}^2 \\ F_{49}=F_{25}^2+F_{24}^2 \\ F_{50}=F_{25}(F_{25}+2F_{24}) \\ F_{100}=F_{50}(F_{50}+2F_{49}) $
というようにできる。細々とした最適化になるが、n<10くらいは先に元の漸化式で計算したほうがいいかもしれない。また、n=24の式を{\scriptstyle F_{24}=F_{12}(2F_{13}-F_{12})}にすることでn=5の時の計算をする必要がなくなる。時間計算量は以下のとおりである。
掛け算 足し算 全体
{\scriptstyle O(1)} {\scriptstyle O(1)} {\scriptstyle O(\log(n))}
{\scriptstyle O((\log(n)^2)} {\scriptstyle O(\log(n))} {\scriptstyle O(n^2)}
{\scriptstyle O(\log(n)\log(\log(n)))} {\scriptstyle O(\log(n))} {\scriptstyle O(n\log(n))}

一般項を用いる・・・しかし

+ ...
フィボナッチ数列は簡単な三項間漸化式であるから、一般項を求めることができて、
$\displaystyle F_n=\frac{\phi ^n-(-\phi)^{-n}}{\sqrt 5},\phi=\frac{1+\sqrt 5}{2} $
となる。分子の第二項はnが十分大きい時、ほぼ0となるので、実用的には
$\displaystyle F_n=\lfloor \frac{\phi ^n}{\sqrt 5}+\frac{1}{2} \rfloor $
で計算できる。
累乗の計算は
$\displaystyle a^{2n}=(a^n)^2,a^{2n+1}=a^na^{n+1} $
と計算できる。これをバイナリ法(または繰り返し二乗法)と呼ぶ。
動的計画法を用いて計算することで時間計算量は、
掛け算 累乗
{\scriptstyle O(1)} {\scriptstyle O(\log(n))}
{\scriptstyle O((\log(n)^2)} {\scriptstyle O(n^2)}
{\scriptstyle O(\log(n)\log(\log(n)))} {\scriptstyle O(n\log(n))}
この時間計算量は先程の改良した漸化式で動的計画法を用いた時と同じなので、問題は
$\displaystyle \phi=\frac{1+\sqrt 5}{2} $
を求めるための時間計算量である。√5を求めるために筆算の開平のようなアルゴリズムを用いていては、十分な精度を得るのに時間がかかりすぎるので、二次収束する近似アルゴリズムであるニュートン法を使う。(ニュートン法の詳しい説明は割愛する)
$\displaystyle a_0=1,a_{n+1}=\frac{a_n^2+1}{2a_n-1} $
とおくと、
$\displaystyle \lim_{n \to \infty}a_n=\phi
となる。ところが、{\scriptstyle a_n=\frac{b_n}{c_n}(b_n\perp c_n)}とおくと、
$\displaystyle \frac{b_{n+1}}{c_{n+1}}=\frac{(\frac{b_n}{c_n})^2+1}{2\frac{b_n}{c_n}-1} \\ =\frac{c_n^2+b_n^2}{2b_nc_n-c_n^2} \\ =\frac{c_n^2+b_n^2}{c_n(2b_n-c_n)} $
であるから、{\scriptstyle b_n=F_{m+1},c_n=F_m}とおくと、
$\displaystyle \frac{b_{n+1}}{c_{n+1}}=\frac{F_m^2+F_{m+1}^2}{F_m(2F_{m+1}-F_m)} \\ =\frac{F_{2m+1}}{F_{2m}} $
となり、実は黄金比の近似計算と先ほどの改良した漸化式と同じ計算をしていることが分かる。ここで、
$\displaystyle b_0=1=F_2,c_0=1=F_1$
であったから、
$\displaystyle b_n=F_{2^n+1},c_n=F_{2^n},a_n=\frac{F_{2^n+1}}{F_{2^n}} $
といえる。この近似法で得た
$\displaystyle \frac{F_{2^n+1}}{F_{2^n}} $
という近似値では、{\scriptstyle m&lt;2^{n+1}}程度でないと正しい値がでないため、あまり意味がない。
よって、一般項を使った方法では、ニュートン法より収束が速いか、時間計算量の少ない近似アルゴリズムを用いるか、累乗と無理数の処理法を考える必要がある。(僕は思いつきませんでした)

行列累乗

+ ...
この記事を公開したあとにいろいろな人から行列累乗について言われたので追記。
フィボナッチ数の漸化式は、行列を用いて次のように表せる。
$\displaystyle \left(\begin{array}{c} F_{n+1}\\F_n \end{array}\right) = \left( \begin{array}{cc} 1&1\\ 1&0 \\ \end{array}\right) \left(\begin{array}{c} F_n\\F_{n-1}\\ \end{array}\right) $
これより、
$\displaystyle \left( \begin{array}{cc} 1&1\\ 1&0\\ \end{array}\right)^n = \left( \begin{array}{cc} F_{n+1}&F_n\\ F_n&F_{n-1} \\ \end{array}\right) $
が得られる。行列の累乗は先に述べたバイナリ法が使えるので、
この方法の時間計算量は
累乗 全体
{\scriptstyle O(1)} {\scriptstyle O(\log(n))}
{\scriptstyle O((\log(n)^2)} {\scriptstyle O(n^2)}
{\scriptstyle O(\log(n)\log(\log(n)))} {\scriptstyle O(n\log(n))}
となる(改良した漸化式で動的計画法を使うのと同じ)

さらなる改良

+ ...
フィボナッチ数とよく似た数列に、リュカ数がある。
$\displaystyle L_{n+2}=L_{n+1}+L_{n},L_0=2,L_1=1 $
というように定義されていて、漸化式はフィボナッチ数と同じで、初期値だけが異なっている。
ここで、次の漸化式を用いることで、改良した漸化式や行列累乗と同じ時間計算量で計算できるうえに、どの漸化式も、大きな数同士の掛け算もしくは足し算は一度しか現れないため、定数倍高速化されることが期待できる。
また、同時にリュカ数を得たい場合には非常に有用である。
$\displaystyle L_{2n}=L_n^2-2(-1)^n \\ F_{2n}=F_nL_n \\ L_{2n+1}=L_nL_{n+1}+(-1)^n \\ F_{2n+1}=L_nF_{n+1}-(-1)^n \\ L_{n+1}=\frac{L_n+5F_n}{2} \\ F_{n+1}=\frac{L_n+F_n}{2} $
2で割る演算は、右ビットシフトを用いて高速に計算できる。
さらに、次のような漸化式も存在し、場合によっては定数倍オーダーで高速になる。(その逆もありうる)
$\displaystyle L_{3n}=L_n(L_n^2-3(-1)^n) \\ F_{3n}=F_n(L_n^2-(-1)^n) $

特殊目的での計算

+ ...
フィボナッチ数の計算結果を利用する際、多くの場合は上何桁、下何桁、あるいはそのフィボナッチ数をある数で割った余りを求められれば十分である。
上の方の数値を利用するとき
上で述べた一般項を、ある程度の精度の{\phi}を求めて計算すれば良い。{\phi}の有効桁数が求める近似値の有効桁数とほぼ等しい(若干小さくなる)。
下の方の数値、余りを利用するとき
改良した漸化式、行列累乗、リュカ数を用いた漸化式のどれかを、各ステップで余りを求めつつ計算すれば良い。このとき、２で割る演算の代わりに割る数を法とした時の2の逆数を掛ける操作をする。その法で2の逆数がないときは、割られる数は偶数になっているので、そのまま2で割る(証明は省略)。
素数で割った余りを求めるとき
  • 2で割るとき
フィボナッチ数を2で割った余りは、{\scriptstyle F_0}から順に{\scriptstyle 0,1,1,0,1,1,0,\ldots}となっているので、n番目のフィボナッチ数を2で割った余りは、(nを3で割った余り)番目のフィボナッチ数と等しい。
  • 5で割るとき
2の場合と同様に、余りの取る値は20回でループするので、n番目のフィボナッチ数を5で割った余りは、(nを20で割った余り)番目のフィボナッチ数を5で割った余りと等しい。
  • その他の素数で割るとき
pを5でない奇素数とするとき、
$\displaystyle F_{mp+n}\equiv\left(\frac{p}{5}\right)F_{m+\left(\frac{p}{5}\right)n} \pmod{a} $
が成り立つ。ここで、{\scriptstyle \left(\frac{p}{5}\right)}は平方剰余記号で、具体的には
$\displaystyle \left(\frac{p}{5}\right)= \begin{cases} 1&(p\equiv1,9\pmod(p)) \\ -1&(p\equiv3,7\pmod(p)) \end{cases} $
という値を取る。
さらに、変形して、
$\displaystyle F_{m\bigl(p-\left(\frac{p}{5}\right)\bigr)+n}\equiv\left(\frac{p}{5}\right)^mF_{n} \pmod{a} $
ここで、{\scriptstyle F_k \bmod p}を求めるとき、{{\scriptstyle k=m\bigl(p-\left(\frac{p}{5}\right)\bigr)+n}}となるm,nを求めれば良い。
さらにmは偶数か奇数かさえ分かればよいので、最初に{{\scriptstyle 2\bigl(p-\left(\frac{p}{5}\right)\bigr)}で割った余りrを求めると、nはすぐに求まり、mは{{\scriptstyle r&lt;p-\left(\frac{p}{5}\right)}であれば偶数、そうでなければ奇数と判定すれば良い。
さらに、kがある数の累乗だとわかっているとき、割った余りを求めるのにもバイナリ法が使える。
これにより、{F_{10^{10^{24}}} \bmod 10^{24}+7}({\scriptstyle10^{24}+7}は素数)もコンピュータを用いればすぐに(多くの場合1msもかからずに)求められることになる。実際に10の(10の24乗)乗番目のフィボナッチ数を計算したとすると、「そのフィボナッチ数の桁数」の桁数でさえ1,000,000,000,000,000,000,000,000桁にもなるので、とても普通の方法では一生かかっても計算できないが、これまで述べた方法を用いることで現実的な計算量にまでなるのである。

計算量の概算

+ ...

漸化式に従った計算

足し算が{\scriptstyle O(1)}
{\scriptstyle F_n}の計算にかかる時間を{\scriptstyle T_n}とすると、詳しい計算は省くが、
$\displaystyle T_1=T_2=O(1),T_{n+2}=T_{n+1}+T_n+O(1) \\ \Rightarrow O(T_n) = O(F_n) $
修正しました(2012/02/18)
となり、フィボナッチ数のオーダーと同じ{\scriptstyle O\Bigl(\bigl(\frac{1+\sqrt 5}{2}\bigr)^n\Bigr)}の時間計算量となる。
足し算が定数時間で終わらない場合は当然時間計算量のオーダーは大きくなる。(はず。求めるのめんどくさいんで誰かやってください、の意)

動的計画法

足し算が定数時間の場合は当然{\scriptstyle O(n)}である。
足し算が桁数に比例した時間がかかる場合(オーダー表記では{\scriptstyle O(\log(n))})、フィボナッチ数は桁数がnにおよそ比例するので、{\scriptstyle O(n^2)}となる。

漸化式の改良以降のアルゴリズム

それぞれの漸化式を使うとこれまでに説明したことで普通に出せる(はず)。
「フィボナッチ数の高速計算」をウィキ内検索
最終更新：2012年04月30日 12:40