はじめに

学部二年の時に動的計画法を学びました。
ただ、その知識も使わないとどんどん忘れていきます。
そんな基本的なところを踏まえているフィボナッチ数列に触れ、思い出すためにこちらの記事を書かせていただきます。

フィボナッチ数列

フィボナッチ数列とは

\begin{align}
F_0 = 0 \quad \quad \quad \quad \quad \quad \quad \\
F_1 = 1 \quad \quad \quad \quad \quad \quad \quad \\
 F_n = F_{n-1} + F_{n-2} \quad (n > 2)\\
\end{align}

で表される数式です。

こちらプログラミングで特にはどのようにすればよろしいのでしょうか。

愚直に考える

fib.py

import time

def fib_stupid(n):
    if(n == 0):
        return 0
    elif(n == 1):
        return 1
    else:
        return fib_stupid(n-1) + fib_stupid(n-2)
        #fib_stupid()を二回呼び出す

start = time.time()
print(fib_stupid(40))
process_time = time.time() - start
print("Time :" + str(process_time))

102334155
Time :75.65844893455505[s]

こちら定義通りですね。
ただこちらの計算のオーダー量はとても大きくなっています。
というのも、fib_stupid(n)を求める際にfib_stupid(n-1),fib_stupid(n-2)を求める必要になり、一回のfib_stupid()で再帰的に二回呼び出すため、結局のところ$O(2^{(n-1)} - 2)$だけ関数を呼び出す羽目になります。現に$F_9$を求めようとすると$F_8$や$F_7$を求める必要があり$F_8$を求めるには$F_7$も求める必要があり$F_7$に関しては二回も計算する必要が出てくるため非常に非効率であります。

ここで効率よく考えるため動的計画法を考えます。

動的計画法 (Dynamic Programing : DP)

動的計画法とは、対象となる問題から複数の小さい問題を見つけ、それらを解くことで元の問題を解く再帰的な構造を持つアルゴリズムである。

以下では、$O(n),O(log(n))$で計算できるアルゴリズムを紹介します。

トップダウン法

fib.py

import sys
sys.setrecursionlimit(30000)
#再帰処理の上限設定

def make_list_for_topdown(n):
        l = [0]
        l.append(1)
        for i in range(n-1):
            l.append(-1)
        return l

lst = make_list_for_topdown(30000)
#今回は3000まで求められるようにlistの形を設定した

def fib_topdown(n):
    if(lst[n] != -1):
        return lst[n]
    else:
        lst[n] = fib_topdown(n-1) + fib_topdown(n-2)
        return lst[n]


start = time.time()
print(fib_topdown(40))
fib_topdown(30000)
process_time = time.time() - start
print("Time :" + str(process_time) + '[s]')

102334155
Time :0.10364508628845215[s]

こちらに関してはlstという配列を作ることで、フィボナッチ数列の値を保存してあります。このことで先ほどのような重複した計算を避けることができます。
関数の呼び出し回数が先ほどのfib_stupid()と変わらないようにも見えますが、配列に結果が残っているため、計算量は$O(2(n-2))$になり先ほどと比べ格段に速くなりました。

ボトムアップ法

fib.py

def fib_bottomup(n):
    if(n < 2):
        return n
    else:
        num_0 = 0
        num_1 = 1
        num_2 = -1
        for i in range(n-1):
            num_2 = num_0 + num_1
            num_0 , num_1 = num_1, num_2
        return num_2

start = time.time()
print(fib_bottomup(40))
fib_bottomup(30000)
process_time = time.time() - start
print("Time :" + str(process_time) + '[s]')

102334155
Time :0.023389101028442383[s]

こちらでは最新の二つの状態をnum_0,num_1で覚えさせて足し合わせていけばいいだけで、再帰構造も持ちません。そのため$O(n-2)$で計算できます。

考察

オーダー的にボトムアップからに比べ、2倍程度の速度を期待していましたが、何故か結果は5倍ほどになりました。listを生成したり読み込んだりする処理の結果がこのような結果になっていると考えれれます(正確に検証するのであれば、処理時間を数回行い平均や分散をとり行うべきですが、話が逸れるため今回は割愛させていてだきます)。

漸化式をといて考える

フィボナッチ数列の一般項は漸化式を解くことで,以下のように表せます。

F_n = \frac{1}{\sqrt{5}}\Biggl\{\biggl( \frac{1+\sqrt{5}}{2}\biggl)^n- \biggl( \frac{1-\sqrt{5}}{2}\biggl)^n\Biggl\}

よって、以下のようにかけます。

fib.py

import math
def fib_general_term(n):
    return math.floor((1/math.sqrt(5))*(((1+math.sqrt(5))/2)**n-((1-math.sqrt(5))/2)**n))
start = time.time()
print(fib_general_term(40))
try:
    fib_general_term(30000)
except OverflowError:
    print("OverflowError")
process_time = time.time() - start
print("Time :" + str(process_time) + '[s]')

102334155
OverflowError
Time :0.00013494491577148438[s]

今回、ルート演算を使用した為小数点が出てきます。その為floor関数で調節いたしました。また、floatの上限は1.7976931348623157e+308であり$F_{30000}$がそれより大きな値のため例外処理を適応いたしました(intには上限がないため計算可能でした)。ただ、それでも上の二つの動的計画法より速いことが考えられます。それは累乗計算のオーダーが$O(log(n))$であることに起因しています。偶数奇数を考えることでn回掛け合わせる処理を短縮しています。似たような例が行列式でも考えられるため、そちらの例の際に説明いたします(勿論例外処理を適応させているため万能なプログラミングとは言えないです)。

行列計算を考える

フィボナッチ数列の漸化式は以下ような行列の積によっても表せます。

\begin{pmatrix}
F_{n} \\
F_{n-1}
\end{pmatrix}
=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}
\begin{pmatrix}
F_{n-1} \\
F_{n-2}
\end{pmatrix}

つまり、

A
:=
\begin{pmatrix}
1 & 1 \\
1 & 0
\end{pmatrix}\\
\begin{align}
\begin{pmatrix}
F_{n} \\
F_{n-1}
\end{pmatrix}
=A
\begin{pmatrix}
F_{n-1} \\
F_{n-2}
\end{pmatrix}
=A^2
\begin{pmatrix}
F_{n-2} \\
F_{n-3}
\end{pmatrix}
=A^{n-1}
\begin{pmatrix}
F_{n-(n-1)} \\
F_{n-1-(n-1)}
\end{pmatrix} \\
=A^{n-1}
\begin{pmatrix}
1 \\
0
\end{pmatrix}
\end{align}

と書けます。(n>0)
$A^n$は$n=2^k$であれば以下のように2乗をk回繰り替えるだけで求まります。
例えば$n = 35$であるならば,
$n = 2^5 + 2^1 + 2^0$のため
$A^{35} = A^{2^5}A^{2^1}A^{2^0} = ((((A^2)^2)^2)^2)^2A^2A$となります
（別の考え方としては、
$35 = 2*17 + 1 = 2*(2*8+1) + 1 = ...$という考え方です。
個人的にはこちらの方が理解しやすいかと思いました）。
以上から分かるように今回は６回すなわち$log(35)$回の計算で終了します。
つまり、今回のアルゴリズム、行列のべき乗によるフィボナッチ数列の計算のオーダーは$log(n)$となります。具体的なコードは以下の通りです(pythonには勿論行列式の演算はnumpyにありますが、それでは先ほどの一般項による求め方と変わらないためあえてnumpyは使用しないものとします)。

fib.py

def make_matrix(row_len,column_len):
    matrix_c = [[0]*column_len for _ in range(row_len)]
    return matrix_c

def matrix_mul(a,b):
    c = make_matrix(len(a),len(b[0]))
    for i in range(len(a)):
        for k in range(len(b)):
            for j in range(len(b[0])):
                c[i][j] +=  a[i][k] * b[k][j]
    return c

def matrix_pow(a,n):
    b = make_matrix(len(a), len(a))
    for i in range(len(a)):
        b[i][i] = 1
    while(n>0):
        if(n & 1):
            b = matrix_mul(a,b)
        else:
            pass
        a = matrix_mul(a,a)
        n >>= 1
    return b

def fib_matrix(n):
    if n == 0:
        return 0
    a = make_matrix(2,2)
    a[0][0] = 1 
    a[0][1] = 1
    a[1][0] = 1 
    a[1][1] = 0
    a = matrix_pow(a, n)
    return a[1][0]

start = time.time()
print(fib_matrix(40))
fib_matrix(30000)
process_time = time.time() - start
print("Time :" + str(process_time) + '[s]')

102334155
Time :0.002933025360107422[s]

こちらでは、ボトムアップ法に比べて約10倍ほど速くなっていることがわかります。

まとめ

以上のような動的計画法を用いることで、計算量は大幅に削減されます。こちらは基本的な事項ですので、後日様々なDPについてまとめたいと思います。

参考

・動的計画法でフィボナッチ数列の計算を速くする。
https://yosuke-furukawa.hatenablog.com/entry/20120120/1327056359