まずは作成したプログラムをペタリ
import time
import numpy as np
from functools import lru_cache
# 関数の実行時間を返す
def measure_time(func, n):
start = time.perf_counter_ns()
func(n)
end = time.perf_counter_ns()
dt = end - start
if dt >= 1000000000:
# ns -> s
dt = str(dt / 1000000000) + " s"
elif dt >= 1000000:
# ns -> ms
dt = str(dt / 1000000) + " ms"
elif dt >= 1000:
# ns -> μs
dt = str(dt / 1000) + " μs"
else:
dt = str(dt) + " ns"
return dt
# 一般的な全探索型再帰関数
def fib1(n):
if n <= 2:
return 1
else:
return fib1(n - 1) + fib1(n - 2)
# キャッシュを伴う再帰関数
@lru_cache(maxsize=None)
def fib2(n):
if n <= 2:
return 1
return fib2(n - 1) + fib2(n - 2)
# ボトムアップ型(所謂動的計画法)関数
def fib3(n):
if n <= 1:
return 1
else:
a, b = 1, 1
for _ in range(2, n):
a, b = b, a + b
return a + b
# 行列の冪乗を用いた関数
def fib4(n):
A = np.matrix([[0, 1],
[1, 1]], dtype=np.int64)
return fast_pow(A, n)[0, 1]
# 一般項(ビネの公式)をダイレクトに使用した関数
def fib5(n):
eq = (1 / np.sqrt(5)) * \
(((1 + np.sqrt(5)) / 2) ** n - ((1 - np.sqrt(5)) / 2) ** n)
return round(eq.real)
# 行列の冪乗計算を早めるアルゴリズム
def fast_pow(x, n):
if n == 0:
return np.matrix([[1, 0], [0, 1]], dtype=np.int64)
elif n % 2:
return x * fast_pow(x, n-1)
else:
y = fast_pow(x, n/2)
return y ** 2
一番上の関数は引数に受け取った関数の処理時間を計測し、桁に合わせてナノ秒、マイクロ秒、ミリ秒か秒に変換して文字列として返します。
一番下の関数は行列の冪乗計算用に定義しました。
たとえば $A^{11}$ を考える時、通常であればコンピュータに対して $A \times A \times A \times A \times A \times A \times A \times A \times A \times A \times A$
と計算させるのも手ですが、この愚直なやり方を工夫すると
$A \times A^{10} $
$A \times (A^5)^2$
$A \times (A \times A^4)^2$
$A \times A^2 \times ((A^2)^2)^2$
と変形でき、多くても2乗までの計算さえすればよくなります。
該当する関数はこの計算を定義したものです。
そして下のfib関数は、上から全探索型の再帰関数、キャッシュを利用した再帰関数、初項から計算するDP的関数、行列の冪乗を用いた関数、最後に一般項をそのまま使った関数です。
fib2関数に使われているLRUキャッシュについてはこちらを参照してください。
まずは40項目の数を求めるのにかかる時間を計測
ちなみに第40項の値は102334155です
print(measure_time(fib1, 40))
print(measure_time(fib2, 40))
print(measure_time(fib3, 40))
print(measure_time(fib4, 40))
print(measure_time(fib5, 40))
11.6705681 s
16.4 μs
4.0 μs
163.7 μs
27.8 μs
再帰的なやりかたは断トツでもっとも時間がかかっています。
最速だったボトムアップ型と比較すると、その処理速度に2917倍の差があります。
続いて100万項目を計算しました。
ただし、fib1関数でこれを実行すると途方もない時間がかかり
また、fib2に関しては、RecursionError: maximum recursion depth exceeded in comparison が発生
fib5に関しては、OverflowError: cannot convert float infinity to integer が発生するため、fib3, fib4でのみ実行した。
ちなみにフィボナッチ数列の第100万項の値は208988桁あります。
print(measure_time(fib3, 1000000))
print(measure_time(fib4, 1000000))
5.8629356 s
247.9 μs
40項の時はボトムアップ型の方が行列計算の関数より速度が速かったが、100万項になると早さが逆転している。というより、行列計算の方は項がだいぶ増えようが処理スピードの増加割合がほとんど変わってないことがわかります。
ボトムアップ型のやり方では、行列計算の関数の23650倍も処理に時間がかかりました。
追記
√5の計算に時間を要するため、fib5を変えた結果も掲示します。
# 一般項(ビネの公式)をダイレクトに使用した関数
def fib5(n):
eq = (1 / 5 ** 0.5) * \
(((1 + 5 ** 0.5) / 2) ** n - ((1 - 5 ** 0.5) / 2) ** n)
return round(eq.real)
print(measure_time(fib1, 40))
print(measure_time(fib2, 40))
print(measure_time(fib3, 40))
print(measure_time(fib4, 40))
print(measure_time(fib5, 40))
11.7035785 s
16.0 μs
3.8 μs
166.8 μs
3.7 μs
27.8 μs -> 3.7 μs
めっちゃ早くなった。