フィボナッチ数列の計算量について

2022 1/6

2022-03-4

マージソートは、与えられたデータを2分割し、その2つの山をそれぞれマージソートを行う。この結果の2つの山の頂上から、大きい方を取り出す…という処理を繰り返すことで、ソートを行う。

第1章フィボナッチ数と行列

まず、行列ではなく素直に正整数 $ a,n $ フィボナッチ数列の計算量についてについて $ a^n $ を求めることを考えます(なお、実際はnが大きくなるにつれ $ a \neq 1 $ のとき $ a^n $ の値もかなり大きくなってしまうので実際はmodを取ることを考えます)。そして、この計算は $ O(\log n) $ で求めることができます。その方法について以下に示していきます( $ a=3 $ , $ n=22 $ とします)。

まず、指数法則より $ 3^22 = (3^11)^2 $ となります。そして、 $ 3^11 $ を指数的に分解していきたいのですが、半分に割ることはできません(かなしい)。ではどうするかというと単純で、 $ 3^11 = 3^1 * 3^10 $ とすれば片方の指数を偶数にすることができました。

同様にして、 $ 3^22 = (3^1 * (3^1 * ((3^1)^2)^2)^2)^2 $ と表すことができます。整理すると、フィボナッチ数列の計算量について $ 3^22 = 3^2 * 3^4 * 3^ $ と表すことができます。これはつまり、指数を2のべき乗の数で分解したということです。さて、ではこの分解した2のべき乗(ここで $ 2^k $ とします)の指数の累乗についてはどのように求めればよいでしょうか？そう、例えば $ 3^2 = (3^1)^2 $ ですし、 $ 3^4 = (3^2)^2 $ です。つまり、 $ 3^(2^k) $ は $ 3 $ を $ \log_ k $ 回2乗した数と考えることができます。このことから $ O(\log k) $ で $ a^n $ の値を求めることができます。これを繰り返し二乗法といいます。コードはこ↑こ↓に書いてありますので是非見てみてください。決して言葉で説明するのがめんどくさくなったというわけではないですよ。

1.5 行列で応用

整数ができるのならば行列でもできそうなことは直感的にわかりますが(本当？)、実際にできます。なので行列累乗のの計算は行列のサイズを無視すると $ O(\log k) $ で求められます。いい話。

1.6 フィボナッチ数で応用

では、このような行列 $ A $ が存在すると嬉しいと思いませんか？

ここで行列 $ A $ について考えましょう。ヒントとして、 $ A $ は2×2の行列です。そして、 $ f_ = f_n + f_,f_n = f_n $ であり、行列積の計算方法を考えると因数分解みたいなイメージで考えることができます。

整数の公式でフィボナッチ数列を求める

よって、もし $m = \left(\begin 1 & 1 \\ フィボナッチ数列の計算量について 1 & 0 \end \right)^n$ なら、 $bn = m\$ になります（Pythonと違って、行列の添え字は通常1が基準になることに注意してください）。 NumPy行列のべき乗が繰り返し二乗法のような振る舞いをすると想定すると、計算量は $O(\mathrm\ n)$ になります。さらに、漸化式を解くために、閉じた式を見つける方法もあります。これにより、次の実数値の公式が導かれます： $\phi = (1 + \sqrt) / 2$ 、 $\psi = (フィボナッチ数列の計算量について 1 - \sqrt) / 2$ とすると、 $\mathrm(n) = (\phi^ フィボナッチ数列の計算量について - \psi^) / \sqrt)$ 。この手法には、任意精度の実数計算を要するという実用上の欠点がありますが、 $n$ の値が小さければ問題はありません。

任意の数列 $an$ の母関数は、 $\Sigman anx^n$ の無限和です。フィボナッチ数列の場合は $\Sigman \mathrm(n)x^n$ になります。つまり、これは無限に続くべき級数であり、 $x^n$ の係数は $n$ 番目のフィボナッチ数に相当します。

この式に $x^$ をかけて $n$ 全体で和をとると、以下の式が得られます。

$F(x)$ を $\mathrm$ の母関数として、それを $フィボナッチ数列の計算量について \Sigma_n\mathrm(n)x^n$ と定義すると、上の式は次のように簡略化できます。

\[F(x) - x - 1 = x(F(x) - 1) + x^2F(x)\]

\[F(x) = xF(x) + x^2F(x) + 1\]

これを $F$ について解くと以下の式が得られます。

整数の公式

まずは、この公式を直感的にとらえるため、 $10^$ で母関数 $F$ を評価してみましょう。

興味深いことに、小数展開した部分に $1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89$ と、フィボナッチ数列が現れています。魔法のような結果に驚いてしまいますが、その理由は次の式から分かります。

$F(10^) = \mathrm(0) + \mathrm(1)/10^3 + \mathrm(2)/10^6 + \mathrm(3)/10^9 + \ldots$

この例では、フィボナッチ数列が次々と $1/1000$ 倍されて並んでいきます。つまり、その値が一旦1000を超えると、隣り合う数に影響を及ぼし始めるということです。この現象は、上記の $F(10^)フィボナッチ数列の計算量について $ の計算で988から確認できます。正しいフィボナッチ数は987ですが、数列の次の数から1だけオーバーフローが発生しています。その結果Off-by-oneエラーが発生し、以降はパターンが崩れてしまうのです。

しかし、いかなる $n$ の値に対しても、10の負の指数を十分大きく取れば、たとえオーバーフローが発生したとしても、 $n$ 番目のフィボナッチ数に悪影響が出ることはありません。ここでは、ある $k$ という値について $10^$ が妥当な値になると仮定しましょう。この値は後ほど選定します。

さらに整数で計算をしたいので（その方がコーディングしやすいので）、全体を $10^$ 倍して $n$ 番目のフィボナッチ数が整数の範囲にくるようにして、式を整理します。

この結果を $10^k$ を法として見ると、 $n$ 番目のフィボナッチ数が得られます（先ほども書きましたが、 $k$ には十分大きな値を選んだものと想定しています）。

あとは $\mathrm(n+1)$ $2^k$ になるように、 $k$ を十分大きく取るだけです。フィボナッチ数列は $\phi^n$ のように増大して、 $\phi$ < $2$ なので、 $k = n+1$ とすれば安全です。

非反復型で閉じた解が得られたのは興味深いですが、これは全く実用的な手法ではありません。ここではサイズが $ O(n^2) $ ビットの整数を用いて、整数演算を実行しています。でも実際には、最終的にビット単位で論理積を取る前に、最初の $ n $ 個のフィボナッチ数が全て連結した整数値を取得しているのです。

最速のフィボナッチ数計算を考える

さらに、フィボナッチ数列の漸化式が $フィボナッチ数列の計算量について F_=F_n+F_$ だったことを思い出せば、\begin
F_&=F_mF_+F_F_n-F_mF_n,\\
F_&=F_mF_n+F_F_
\endという2本の式が得られる。つまり、行列の乗算の代わりに\[(F_m,F_)\cdot(F_n,F_)\mapsto(F_mF_+F_F_n-F_mF_n,F_mF_n+F_F_)\]という演算を使って計算できる。

$(F_n,F_)$ の形の組を FibPair と呼ぶことにすれば、 FibPair フィボナッチ数列の計算量についてには\[(a,b)\cdot(a’,b’)\mapsto(a b’+ba’-aa’,aa’+bb’)\]によるモノイド演算（単位元は $(F_0,F_1)=(0,1)$）が定まる（結合法則の確認は読者の演習問題とする）。Haskellでの実装例は次のようになる：

ここで、 stimesMonoid は、モノイドの $フィボナッチ数列の計算量についてフィボナッチ数列の計算量について n$ 乗を $n$ の二進展開を利用して高速に（O(log n)で）計算してくれる関数である（同様のアルゴリズムは、行列の $n$ フィボナッチ数列の計算量についてフィボナッチ数列の計算量について乗の計算にも利用していた）。

このFibPairを使って 123456789 （およそいちおく）番目のフィボナッチ数を計算した場合と、一般の2×2行列の方法 (MatFib) での計算時間を比較してみる。

整数 n に対して n 番目（および n+1 番目）のフィボナッチ数 $(F_n,F_)$ を対応させる写像は、整数の加法に関するモノイドから、FibPairへのモノイド準同型だとみなせる。

Fast doubling

Haskellのべき乗 (^) やモノイドのn倍 stimesMonoid では、nを二進展開したものを「右から左に」辿って計算する、ということをやっている。

例えば、 a の19乗（二進法で10011）を計算したいというときは、アルゴリズム的には次のようなことを行なっている：

Y := 1, Z := a とおく。
n=19 の最下位ビット（1の位）は 1 なので、Y に Z をかける（Y := Y * Z）。
Z を自乗する（Z := Z * Z）。
n=19 の下から2番目のビット（2の位）は 1 なので、 Y フィボナッチ数列の計算量についてに Z をかける（Y := Y * Z）。
Z を自乗する（Z := Z フィボナッチ数列の計算量について * Z）。
n=19 の下から3番目のビット（4の位）は 0 なので、 Y には何もしない。
Z を自乗する（Z := フィボナッチ数列の計算量について Z * Z）。
n=19 の下から4番目のビット（8の位）は 0 なので、 Y には何もしない。
Z を自乗する（Z :フィボナッチ数列の計算量について = Z * Z）。
n=19 の下から5番目のビット（16の位）は 1 なので、 Y に Z をかける（Y := Y * Z）。
n は二進法で5桁なので、これでアルゴリズムを終了する。 Y は a フィボナッチ数列の計算量について の19乗である。

一つの式で書けば、\[a^=a\cdot a^2\cdot (((a^2)^2)^2)^2\]と計算していることになるだろう。

Y := 1 とおく。
n=19 の最上位ビット（16の位）は 1 なので、 Y に a をかける（Y := Y * a）。
Y を自乗する（Y := Y * Y）。
n=19 の上から2番目のビット（8の位）は 0 なので、 Y フィボナッチ数列の計算量についてには何もしない。
Y を自乗する（Y := Y * Y）。
n=19 の上から3番目のビット（4の位）は 0 なので、フィボナッチ数列の計算量についてフィボナッチ数列の計算量について Y には何もしない。
Y を自乗する（Y := Y * Y）。
n=19 の上から4番目のビット（2の位）は 1 フィボナッチ数列の計算量についてなので、 Y に a をかける（Y := Y * a）。
Y を自乗する（Y := Y * Y）。
n=19 の上から5番目のビット（2の位）は 1 なので、 Y に a をかける（Y := Y * a）。
n は二進法で5桁なので、これでアルゴリズムを終了する。 Y は a の19乗である。

一つの式で書けば、\[a^=(((a^2)^2)^2\cdot a)^2\cdot a\]と計算していることになるだろう。

「右から左」と「左から右」の比較だが、計算機上では整数の二進表記を「右から左」に辿る方が実装しやすい（ひたすら2で割って、あまりを見れば良い）。そのため、「右から左」が使われることが多い。実際、Haskellのべき乗 (^) や stimesMonoid もそうなっている。

（二進表記の「右から左」と「左から右」のアルゴリズムについては、The Art of Computer Programming Vol. 2 に記載がある）

さて、「左から右」には、「右から左」にはない特徴がある。それは、アルゴリズム中で使う乗算が「自乗」と「a をかける」の2種類だけ、ということだ。この特徴と、フィボナッチ数計算 (FibPair) の事情を組み合わせるとどうなるか。

n乗のアルゴリズムをフィボナッチ数計算に使う場合、 a としては組 $\operatorname(1)=(F_1,F_2)=(1,1)$ フィボナッチ数列の計算量についてを用いる。そして、「a をかける」という操作は、「次のフィボナッチ数を計算する」ということであり、足し算1回でできてしまう：\[(F_n,F_)\cdot (1,1)=(F_,F_n+F_)\]「右から左」の場合は（自乗のほかに）フィボナッチ数列の計算量について FibPair の演算が3回必要だったところ、「左から右」なら多倍長整数の加算3回で済んでしまうのだから、「左から右」の方が有利である。

この「左から右」に辿る方法をフィボナッチ数に適用したものは、fast doubling と呼ばれているようだ。二進表記を「左から右」に辿るため、実装の際に末尾ではない再帰呼び出しを使うことになる。

Fast doubling をHaskellで実装すると次のようになる：

この時点ですでに FibPair + stimesMonoid よりも早くなっているが、いくつかの小細工を加えるとさらに早くなる。

まず、ここまできたらもはや FibPair フィボナッチ数列の計算量についての汎用的なモノイド演算は必要ない。モノイドとしての自乗 p <> p さえ計算できれば十分で、これは\[(a,b)\cdot(a,b)=(2ab-a^2,a^2+b^2)=(a(2b-a),a^2+b^2)\]で計算できる。モノイド演算では多倍長整数の乗算が4回、加減算が3回だったのが、自乗であれば多倍長整数の乗算が4回、加減算が2回となる。

そして、「FibPair を自乗してから次のフィボナッチ数を計算する」部分はひとまとめにできる。つまり、コード上は p <> p と FibPair b (a + b) に分かれていたのを、

実験（実行時間の計測）

筆者の環境（MacBook Pro (Late 2013), GHC 8.6.3）では、行列を使った版で123456789番目のフィボナッチ数を計算したところ11秒程度、FibPairとstimesMonoidを使った版で同じ計算をしたところ5秒程度、fast doublingを使った版（最後に載せた、タプルを使ったコード）では1.7秒程度かかった。

なお、記事の最初の方で書いた、一般項を $\mathbf(\sqrt)$ で計算するやつ（$\left(\frac<1+\sqrt>\right)^n$ の計算だけで済ませる）は、5秒程度だった。こちらは有理数計算を伴うため、整数演算のみのFastDoublingには勝てないのだろう（推測）。

結論として、この記事で触れたアルゴリズムの中では fast doubling が一番早い。

おまけ：多倍長計算の特性を考慮する

演算回数でいうと、A, B どちらも乗算が4回、加減算が3回ずつである。では、どちらを使っても同じなのだろうか。この辺の話は、演算対象となる型によって変わってくる。

仮に、変数 a , a' , b , b' のいずれもn桁の整数だとしよう。すると、Aの方は

n桁の整数同士の乗算が4回
2n桁の整数同士の加減算が3回

n桁の整数同士の乗算が4回
2n桁の整数同士の加減算が2回
n桁の整数同士の加減算が1回

余談：HaskellのSemigroup/Monoidには自乗するためのメソッドがない

モノイドの元のn乗を計算する際、「左から右」「右から左」のいずれも、モノイドの元の自乗 $x^2$ の計算を利用した。Data.SemigroupやData.Monoidには積を計算するメソッド <> はあるが、自乗に特化されたメソッドはないので、自乗の計算には積演算 <> が使われることになる。

自乗に特化されたメソッドがあると何が嬉しいかというと、インスタンスによっては一般の積演算 <> よりも高速な実装を提供できる可能性があり、 stimesMonoid のような関数がそれを利用できるようになることである。 sconcat や stimes がクラス外の関数ではなくて Semigroup クラスの中で定義されているのと同じ理由である。もちろん「自乗するメソッド」のデフォルト実装は \x -> x <> x とする。

コードで書けば、 Semigroup クラスがこういう風に定義されていてほしい：

FibPairの場合に、この「自乗する」 twice メソッドがあるとどの程度嬉しいか。FibPairのモノイド演算は次のように定義されていた：

（fast doublingの説明で同じことを書いたが）こちらは多倍長整数の乗算が4回、多倍長整数の加減算が2回で、加減算が1回減った上に乗算1回分が定数 2 の掛け算に変わっている。

先ほどリンクを貼った実験プログラムでは、 FibPair の自乗で特化したコードを使うようにしたものを FibPairX という名前（名付けが雑だ）で実装しており、それぞれ234567890番目のフィボナッチ数を計算させてみると

という風に FibPairX の方が若干早い。

（このことから、GHCの最適化器は \x -> x <> x というコードを式変形して上記の twice の実装に変えるようなことは行わない、ということがわかる。数学的に同値な式でも多倍長整数のコストを考えると優劣があるのは「おまけ：多倍長計算の特性を考慮する」で見た通りなので、コンパイラーが勝手にそういう式変形を行わないことはプログラマーにとっては生成コードを予測できるということであり、良いことなのだが。）

線形漸化式と母関数

体 $\mathbb$ フィボナッチ数列の計算量について 上の数列 $a_0,\,a_1,\,\dotsc,$ が $d+1$ 項間線形漸化式を満たすとは、ある $(c_i\in\mathbb)_$ が存在して $$ a_n = \sum_^d c_i a_ $$ が任意の $n\ge d$ について成り立つことを言う。数列の母関数は形式的冪級数として $$ G(x) := \sum_ a_i x^i フィボナッチ数列の計算量について $$ と定義される。母関数 $G(x)$ に対して $Q(x) := 1-\sum_^d c_i x^i$ を掛けると、$x^d$ 以上の項はすべて消えるので、高々 $d-1$ 次の多項式 $P(x)$ となる。よって $$ G(x) = \frac $$ となる。また、逆に高々 $d-1$ 次の多項式 $P(x)フィボナッチ数列の計算量について $ と$d$ 次の多項式 $Q(x)$ が $Q(0)=1$ フィボナッチ数列の計算量についてを満たし、$G(x) = P(x)/Q(x)$ となるとき、 $G(x)フィボナッチ数列の計算量について $ が $d+1$ 項間線形漸化式を満たす数列の母関数となることは $G(x)Q(x)$ の $d$ 次以上が0になることから簡単に確認できる。

アルゴリズム

アルゴリズムの実装上の工夫

高速フーリエ変換が使える体 $\mathbb$ の場合、$\mathsf(d) = O(d\log d)$ となる。時間間引き、周波数間引きの高速フーリエ変換の式を考慮することにより、漸化式の1回の適用の中で計算する高速フーリエ変換(計4回)のサイズを半分にできる。その結果、計算量は$(2/3)\mathsf(d)\lceil\log(N+1)\rceil$ となる。実装例はこちら。ただし、周波数領域ではビットリバーサル置換がかかっているので少し読み辛い。明文化されているのを見たことがないが、FFTダブリング(長さ $k$ のDFTを長さ $2k$ のDFTに変換する手続き)をする際には、周波数領域でビットリバーサル置換をかかっていると DFTの偶数番目が前半に連続して並んでいるので都合がよい。 © 2018 Ryuhei Mori このページについての連絡は email: mori[at]c.titech.ac.jp

タグアーカイブ: 再帰方程式

ハノイの塔は、３本の塔にN枚のディスクを積み、(1)1回の移動ではフィボナッチ数列の計算量についてディスクを1枚しか動かせない、(2)ディスクの上により大きいディスクを積まない …という条件で、山積みのディスクを目的の山に移動させるパズル。

一般解の予想

ハノイの塔の移動回数をとした場合、少ない枚数での回数の考察から、以下の一般式で表せることが予想できる。

再帰方程式

ということが言える。(フィボナッチ数列の計算量についてこれがハノイの塔の移動回数の再帰方程式)
ディスクが枚の時、予想が正しいのは明らか①,②。
ディスクが枚で、予想が正しいと仮定すると、枚では、

となり、枚でも、予想が正しいことが証明された。よって 数学的帰納法 により、1枚以上で予想が常に成り立つことが証明できた。

理解度確認

の「ピラミッドの体積」pyra() を、ループにより計算するプログラムを記述せよ。での2分探索法のプログラムを、再帰によって記述せよ。(以下のプログラムを参考に)。また、このプログラムの処理時間にふさわしい再帰方程式を示せ。

再帰を使ったソートアルゴリズムの分析

参考：ソートアルゴリズム12種を可視化してみた

この中で、高速なソートアルゴリズムは、クイックソート(最速のアルゴリズム)とマージソート(オーダでは同程度だが若干効率が悪い)であるが、ここでは、再帰方程式で処理時間をイメージしやすい、マージソートにて説明を行う。

マージソートの分析

この再帰方程式を、N=1,2,4,8…と代入を繰り返していくと、最終的に処理時間のオーダがとなる。

選択法とクイックソートの処理時間の比較

データ数 N = 20 件でソート処理の時間を計測したら、選択法で 10msec 、クイックソートで 20msec であった。

データ件数 N = 100 件では、選択法,クイックソートは、それぞれどの程度の時間がかかるか答えよ。
データ件数何件以上なら、クイックソートの方が高速になるか答えよ。

再帰呼び出しの処理時間の見積もり

前回の授業の復習と練習問題

前回の授業では、for ループによる繰り返し処理のプログラムについて、処理時間を T(N) の一般式で表現することを説明し、それを用いたオーダー記法について説明を行った。理解を確認するための練習問題を以下に示す。

ある処理のデータ数Nに対する処理時間が、であった場合、オーダー記法で書くとどうなるか？
の処理時間を要するアルゴリズムを、オーダー記法で書くとどうなるか？また、このような処理時間となるアルゴリズムの例を答えよ。
の処理時間を要するアルゴリズムを、オーダー記法で書くとどうなるか？
(ヒント: ロピタルの定理)

1は、N→∞において、N 2 ≪ 2 N なので、O(2 N ) 。厳密に回答するなら、練習問題3と同様の証明を行うべき。
2は、O(1)。誤答の例：O(0)と書いちゃうと、T(N)=T_α×0=0になってしまう。事例は、電話番号を、巨大配列の”電話番号”番目の場所に記憶するといった方法。(これはハッシュ法で改めて講義予定)

再帰呼び出しの基本

再帰関数は、自分自身の処理の中に「問題を小さくした」自分自身の呼び出しを含む関数。プログラムには問題が最小となった時の処理があることで、再帰の繰り返しが止まる。

階乗 fact(N) を求める処理は、以下の様に再帰が進む。

また、フィボナッチ数列 fib(N) を求める処理は以下の様に再帰が進む。

再帰呼び出しの処理時間

> 再帰方程式

このような、式の定義自体を再帰を使って表した式は再帰方程式と呼ばれる。これを以下のような代入の繰り返しによって解けば、一般式 が得られる。

一般的に、再帰呼び出しプログラムは(考え方に慣れれば) 分かりやすくプログラムが書けるが、プログラムを実行する時には、局所変数や関数の戻り先を覚える必要があり、深い再帰ではメモリ使用量が多くなる。
ただし、fact() や pyra() のような関数は、プログラムの末端で再帰が行われている。(fib()は、再帰の一方が末尾ではない)
このような再帰は、末尾再帰(tail recursion) と呼ばれ、関数呼び出しの return を、再帰処理の先頭への goto 文に書き換えるといった最適化が可能である。言い換えるならば、末尾再帰の処理は繰り返し処理に書き換えが可能である。このため、末尾再帰の処理をループにすれば再帰のメモリ使用量の問題を克服できる。