はじめに
これはPython その2 Advent Calendar 2018の22日目の記事となります。
ある機械学習の本を読んでいたら、フィボナッチ数列に一般項があることを知った。
検索してみたら下記サイトを見つけたがHaskellで実装されていたので、これをPythonでやってみたくなった。
レオナルド・フィボナッチ
名前の由来となったイタリアの数学者レオナルド・フィボナッチですが、32歳の1202年に「算盤の書」を出版、この中のインドの方法としてアラビア数字(0-9)を紹介。これがヨーロッパの知識層へ広く受け入れられ、ローマ数字からアラビア数字へと一般に常用されるには15世紀頃までかかりました。
1450年頃にグーテンベルクにより活版印刷が実用化され印刷物が大量に造られるようになったことが普及に大きく影響している。
フィボナッチ数列
フィボナッチ数とは、前の2つの数を足したものが次の数になるという規則に基づいている数列です。
例 1,1,2,3,5,8,13,21,34,55,89144…と続く
プログラミング言語の実装では、再帰的処理の例としてよく紹介されています。
漸化式
漸化式(ぜんかしき)とは再帰関係式である。
\begin{eqnarray*}
F_1 &=& 1 \\
F_2 &=& 1 \\
F_{n+2} &=& F_n + F_{n+1} \ \ (1 \leq n)
\end{eqnarray*}
一般項
整数の数列なのに一般項に無理数入ってる・・・
F_n = \frac{1}{\sqrt{5}}\biggl\{\biggl(\frac{1 + \sqrt{5}}{2}\biggr)^n - \biggl(\frac{1 - \sqrt{5}}{2}\biggr)^n\biggr\}
一般項にルート5が現れる理由
整数の数列なのに無理数の $\sqrt{5}$ が現れるのは、二次方程式の解の公式が関係しています。
それを説明するには、等比数列の一般式の理解が必要となります。
等比数列の一般式
例えば「3, 6, 12, 24, 48...」の数列で考えていきましょう。初項は3、公比は2です。
一般項、つまりn番目の項は「初項3に公比2をn-1回かけた数」なので
$a_{n} = 3 \times 2^{n-1}$ となります。
$n=2$ の時、$3 \times 2^{2-1}=3 \times 2 = 6$
$n=3$ の時、$3 \times 2^{3-1}=3 \times 4 = 12$
$n=4$ の時、$3 \times 2^{4-1}=3 \times 8 = 24$
これを一般化すると、初項 $a$、公比 $r$ の等比数列における一般項は
$a_{n} = a \times r^{n-1}$ となります。
フィボナッチ数の等比数列
等比数列の一般項は、$a_{n} = ar^{n−1}$ ですから、次項は $a_{n+1} = ar^{n}$、次の次項 $a_{n+2} = ar^{n+1}$ です。
次のフィボナッチ数 $a_{n+2}$ は、1つ前の値 $a_{n+1}$ と 2つ前の値 $a_{n}$ を足した値となるので、
$ar^{n+1} = ar^{n} + ar^{n−1}$ となります。例 $5 = 3 + 2$、$8 = 5 + 3$
一般に $a$、$r \neq 0$ だから両辺を $ar^{n−1}$ で割って
$\displaystyle \frac{ar^{n+1}}{ar^{n−1}} = \frac{ar^{n}}{ar^{n−1}} + \frac{ar^{n−1}}{ar^{n−1}}$
割り算のときは指数部分は引き算になる。
$r^{(n+1)-(n−1)} = r^{n-(n−1)} + r^{(n−1)-(n−1)}$
指数nは消える。
$r^{n-n+1+1} = r^{n-n+1} + r^{n-n−1+1}$
$r^{1+1} = r^{1} + r^{−1+1}$
$r^{2} = r^{1} + r^{0}$
$r^{2} = r + 1$ が求まります。
これを二次方程式に変形すると
$r^{2}-r-1 = 0$
二次方程式の解の公式に当てはめる
$ax^{2}+bx+c=0$
二次方程式の解の公式
$\displaystyle x=\frac{-b\pm \sqrt{b^{2}-4ac}}{2a}$
上記で求めたフィボナッチ数の等比数列 $r^{2}-r-1 = 0$ を
二次方程式の解の公式に当てはめます。
$a=1$、$b=-1$、$c=-1$ なので式にセットします。
$\displaystyle r=\frac{-(-1)\pm \sqrt{(-1)^{2}-4\times1\times(-1)}}{2\times1}$
$\displaystyle r=\frac{1\pm \sqrt{1+4}}{2}$
$\displaystyle r=\frac{1\pm \sqrt{5}}{2}$
ただ、これだけだと先頭の $1$ と $1$ が満たせないので、もう一工夫必要です。
解の公式により解が2つ求まったので、$α$ と $β$ にして、$α-β$ とすると $\sqrt{5}$ がそれぞれに付きます。あとは、$\sqrt{5}$ で割ればフィボナッチ数列が求まります。
詳しく下記を参照してください。
普通に実装
import numpy as np
def fib(n):
f = (((1 + np.sqrt(5)) / 2)**n - ((1 - np.sqrt(5)) / 2)**n ) / np.sqrt(5)
return f
for n in range(1, 10 + 1):
print(fib(n))
10までの結果は55なので、Int型にすると55になります。
※ np.sqrt を使わないでPythonの基本機能で実現する場合、5**0.5 で平方根の演算ができます。
1.0
1.0
2.0
3.0000000000000004
5.000000000000001
8.000000000000002
13.000000000000002
21.000000000000004
34.00000000000001
55.000000000000014
次は100までの結果を求めます。
print(int(fib(100)))
# 354224848179263111168
100 番目の本当のフィボナッチ数は $354224848179261915075$ です。
残念ながら値が違ってしまいますね。
$354224848179263111168$
$354224848179261915075$
たかが 100 番目くらいで正しい答えを得られないようでは"なってない"といわざるを得ないではない。
計算方法を変える
では、どうすればいいのだろうとネット検索すると下記サイトが参考になりました。
計算方針を考える
落ち着こう。確かに無理数は避けて通れない。しかし、我々は無理数どころか、"もっとスゴい"数を扱ったことがあるはずだ。そう、二乗すると -1 になっちゃう愛に溢れた"アレ"を含む数、複素数である。
複素数をどう扱ったか思い出そう。次のように表現したはずだ。
$a+bi$
虚数単位 $i$ はもはや無理数ですらない。それでも我々は扱ってきた。ならばたかが $\sqrt{5}$ 程度が扱えないわけがない。
というわけで、複素数の i を $\sqrt{5}$ に置き換えた形で今回の式を取り扱うことにしよう。
フィボナッチ数列の一般項を計算する(※ただし有理数に限る)
Haskellでは標準で Rational という有理数を扱う型があるが、Pythonでは標準 fractions を使うと分数(有理数)での計算ができる。
分子と分母をそれぞれ整数で指定する。
from fractions import Fraction
print(Fraction(1, 3))
# 1/3
FibNum こと Rational の二要素からなるタプルは、左に $\sqrt{5}$ が付かない項を、右に $\sqrt{5}$ が付く項を格納することしよう。つまり (1, 1) と書けば $1+\sqrt{5}$ のこと。(0, 1) と書けば $\sqrt{5}$ のこと。(1 % 2, 1 % 3) と書けば $\frac{1}{2}+\frac{1}{3}\sqrt{5}$ のことを表す。(中略)
除算は 0 で割れないとか面倒なこともあるので、乗算の形にしておきたい。
まず、(1, 1) `fibDiv` (2, 0) は要するに $\frac{1+\sqrt{5}}{2}$ のことだが、こんなものは $\frac{1}{2}+\frac{1}{2}\sqrt{5}$、つまり (1 % 2, 1 % 2) としてしまえば良い。
後ろ側の $\sqrt{5}$ で割る処理は、逆数であるところの $\frac{1}{\sqrt{5}}$ 、つまり $\frac{\sqrt{5}}{5}$ を掛ければ良い。$\frac{\sqrt{5}}{5}$ ってことは $0+\frac{1}{5}\sqrt{5}$ だから、ここでの表現では (0, 1 % 5) ってことだ。
HaskellとPythonのどちらもおぼつかないが、がんばって移植してみる。
演算子を実装する : 乗算
FibNum の乗算とは何かと言うと、$(a+b\sqrt{5})(c+d\sqrt{5})$ である。
\begin{eqnarray*}
&=&(a+b\sqrt{5})(c+d\sqrt{5}) \\
&=&ac+ad\sqrt{5}+bc\sqrt{5}+5bd \\
&=&(ac+5bd)+(ad+bc)\sqrt{5}
\end{eqnarray*}
fibMul :: FibNum -> FibNum -> FibNum
fibMul (a, b) (c, d) = (a*c + 5*b*d, a*d + b*c)
Pythonに移植
def fibMul(fn1, fn2):
a = fn1[0]
b = fn1[1]
c = fn2[0]
d = fn2[1]
return (a*c + 5*b*d, a*d + b*c)
演算子を実装する : 累乗
fibExp :: FibNum -> Integer -> FibNum
fibExp _ 0 = (1, 0)
fibExp a 1 = a
fibExp a n = fibMul a (fibExp a (n - 1))
Pythonに移植
def fibExp(fn, n):
if n == 0:
return (Fraction(1), Fraction(1))
elif n == 1:
return fn
return fibMul(fn, fibExp(fn, n - 1))
高速化版
時間がなかったので今回は累乗の高速化版は省きました。
2018/12/29 高速化版を追記
逐次平方を使えばより少ないステップでべき乗が計算出来る。
例えば $b^8$ の計算は
$b・(b・(b・(b・(b・(b・(b・b)))))))$
とするよりは乗算三回で
$b2 = b・ b$
$b4 = b2・ b2$
$b8 = b4・ b4$
と計算する。つまり、指数値 n が偶数の時には2乗の値を使用する。
def fastExp(a, n):
if (n == 0):
return 1;
if (n % 2 == 0):
return fastExp(a, n / 2) ** 2;
return a * fastExp(a, n - 1);
print(fastExp(2, 8));
#256
print(pow(2,8));
#256
普通のべき乗の計算なら pow 関数を使えばいいが、今回はタプル の Fraction を使っているので考え方を理解して作る必要がある。
fibExp a n
| even n = ex
| otherwise = ex `fibMul` a
where ex = fibExp (fibMul a a) (n `div` 2)
Pythonに移植
def fibExp(fn, n):
if (n == 0):
return (Fraction(1), Fraction(1));
if (n % 2 == 0):
return fibExp(fibMul(fn, fn), n / 2);
return fibMul(fn, fibExp(fn, n - 1));
演算子を実装する : 減算
fibMin :: FibNum -> FibNum -> FibNum
fibMin (a, b) (c, d) = (a - c, b - d)
Pythonに移植
def fibMin(fn1, fn2):
a = fn1[0]
b = fn1[1]
c = fn2[0]
d = fn2[1]
return (a - c, b - d)
メインを実装する
fibMain :: Integer -> FibNum
fibMain n = (fibExp (1 % 2, 1 % 2) n `fibMin` fibExp (1 % 2, -1 % 2) n) `fibMul` (0, 1 % 5)
Pythonに移植
def fibMain(n):
fibNum1 = (Fraction(1, 2), Fraction( 1, 2))
fibNum2 = (Fraction(1, 2), Fraction(-1, 2))
fibNum3 = fibMin(fibExp(fibNum1, n), fibExp(fibNum2, n))
fibNum4 = (Fraction(0, 1), Fraction(1, 5))
return fibMul(fibNum3, fibNum4)
使ってみよう！
for x in range(1, 10 + 1):
print(fibMain(x))
(Fraction(1, 1), Fraction(0, 1))
(Fraction(1, 1), Fraction(0, 1))
(Fraction(2, 1), Fraction(0, 1))
(Fraction(3, 1), Fraction(0, 1))
(Fraction(5, 1), Fraction(0, 1))
(Fraction(8, 1), Fraction(0, 1))
(Fraction(13, 1), Fraction(0, 1))
(Fraction(21, 1), Fraction(0, 1))
(Fraction(34, 1), Fraction(0, 1))
(Fraction(55, 1), Fraction(0, 1))
フィボナッチ数列は元々整数の数列なので、当然計算結果に無理数であるところの $\sqrt{5}$ は絶対に含まれない。だから、 $a+b\sqrt{5}$ の $b$ は絶対に $0$ になるのだ。それに、整数ということは有理数 $a$ の分母も絶対に $1$ になる。余分なコトを出力すると見づらいので、その辺は削ってしまおう。
def fib(n):
return fibMain(n)[0].numerator
for n in range(1, 10 + 1):
print(fib(n))
表示がすっきりした。
1
1
2
3
5
8
13
21
34
55
実力を確認する
先程違っていた 100 番目の結果
print(fib(100))
# 354224848179261915075
100 番目のフィボナッチ数 $354224848179261915075$ で一致しました。
やったね。
10,000 番目だっていけると思ったらエラーだ。
print(fib(10000))
RecursionError: maximum recursion depth exceeded in comparison
どうやら再帰の上限を超えてしまうようだ。
Python: maximum recursion depth exceeded while calling a Python object - stackoverflow
上限を追加すればいいようだ。
import sys
sys.setrecursionlimit(15000)
print(fib(10000))
# 2090桁ある
# 3364476487643178326662161200510754331030214846068006390656476997468008144216666236815559551363373402...(中略)...366875
完成品
2018/12/30にfibExpを高速化版に変更
from fractions import Fraction
import sys
sys.setrecursionlimit(15000)
def fibMul(fn1, fn2):
a = fn1[0]
b = fn1[1]
c = fn2[0]
d = fn2[1]
return (a*c + 5*b*d, a*d + b*c)
def fibExp(fn, n):
if (n == 0):
return (Fraction(1), Fraction(1));
if (n % 2 == 0):
return fibExp(fibMul(fn, fn), n / 2);
return fibMul(fn, fibExp(fn, n - 1));
def fibMin(fn1, fn2):
a = fn1[0]
b = fn1[1]
c = fn2[0]
d = fn2[1]
return (a - c, b - d)
def fibMain(n):
fibNum1 = (Fraction(1, 2), Fraction( 1, 2))
fibNum2 = (Fraction(1, 2), Fraction(-1, 2))
fibNum3 = fibMin(fibExp(fibNum1, n), fibExp(fibNum2, n))
fibNum4 = (Fraction(0, 1), Fraction(1, 5))
return fibMul(fibNum3, fibNum4)
def fib(n):
return fibMain(n)[0].numerator
print(fib(100))
最後に
HaskellとPythonのどちらもおぼつかない状態で、とりあえず出来て良かった。
考え方だけでいえば、分数(Fraction)ではなく複素数(Complex)でも出来そうな気がしますね。
複素数版
出来そうな気がしたのでやってみました。
from fractions import Fraction
def fibMul(fn1, fn2):
a = fn1.real
b = fn1.imag
c = fn2.real
d = fn2.imag
return complex(a*c + 5*b*d, a*d + b*c)
def fibExp(fn, n):
if n == 0:
return (1 + 1j)
elif n == 1:
return fn
return fibMul(fn, fibExp(fn, n - 1))
def fibMin(fn1, fn2):
return fn1 - fn2
def fibMain(n):
fibNum1 = 1/2 + 1/2j
fibNum2 = 1/2 - 1/2j
fibNum3 = fibMin(fibExp(fibNum1, n), fibExp(fibNum2, n))
fibNum4 = 0 + 1/5j
return fibMul(fibNum3, fibNum4)
def fib(n):
return int(fibMain(n).real)
print(fib(100))
# 354224848179261865984
10番目のフィボナッチ数は55になったので100番目もいけるかと思ったのですが、残念ながら結果が一致しませんでした。