フィボナッチ数列の一般項について

# はじめに
これは[Python その2 Advent Calendar 2018](https://qiita.com/advent-calendar/2018/python_python)の22日目の記事となります。

ある機械学習の本を読んでいたら、フィボナッチ数列に一般項があることを知った。
検索してみたら下記サイトを見つけたがHaskellで実装されていたので、これをPythonでやってみたくなった。

* [フィボナッチ数列の一般項を計算する(※ただし有理数に限る)](http://labs.timedia.co.jp/2012/11/fibonacci-general-term-using-rational.html)
* [フィボナッチ数列の一般項をもう少しだけ速く計算する(ライブラリ活用編)](http://labs.timedia.co.jp/2012/12/fibonacci-general-term-using-library.html)
* [フィボナッチ数を一般項から求めてみる。](http://tsumuji.cocolog-nifty.com/tsumuji/2013/11/post-4e7c.html)

# レオナルド・フィボナッチ
名前の由来となったイタリアの数学者[レオナルド・フィボナッチ](https://ja.wikipedia.org/wiki/%E3%83%AC%E3%82%AA%E3%83%8A%E3%83%AB%E3%83%89%E3%83%BB%E3%83%95%E3%82%A3%E3%83%9C%E3%83%8A%E3%83%83%E3%83%81)ですが、32歳の1202年に「算盤の書」を出版、この中のインドの方法としてアラビア数字(0-9)を紹介。これがヨーロッパの知識層へ広く受け入れられ、ローマ数字からアラビア数字へと一般に常用されるには15世紀頃までかかりました。
1450年頃にグーテンベルクにより活版印刷が実用化され印刷物が大量に造られるようになったことが普及に大きく影響している。

* [ローマ数字のⅣとⅥ、ⅨとⅪを混同しなくなる方法](https://qiita.com/yaju/items/cd578d1a7b7bc9602156)
* [記数法のかんたんな歴史 ― アレッサンドラ・キング - TED Ed](https://www.youtube.com/watch?v=cZH0YnFpjwU)

# フィボナッチ数列
フィボナッチ数とは、前の2つの数を足したものが次の数になるという規則に基づいている数列です。
例 1,1,2,3,5,8,13,21,34,55,89144…と続く

プログラミング言語の実装では、再帰的処理の例としてよく紹介されています。

## 漸化式
漸化式(ぜんかしき)とは再帰関係式である。

```math
\begin{eqnarray*}
F_1 &=& 1 \\
F_2 &=& 1 \\
F_{n+2} &=& F_n + F_{n+1} \ \ (1 \leq n)
\end{eqnarray*}
```

## 一般項
整数の数列なのに一般項に無理数入ってる・・・

```math
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}$ で割ればフィボナッチ数列が求まります。
詳しく下記を参照してください。

![Fibonacci.png](https://qiita-image-store.s3.amazonaws.com/0/14053/e5e4e756-b760-5cf5-6275-099913bc0929.png)

* [「ダヴィンチコード」に登場する数学 - pdf](http://www.gifu-nct.ac.jp/sizen/okada/DaVinci.pdf)
* [フィボナッチ数列の漸化式を黄金比で求める方法 - YouTube](https://www.youtube.com/watch?v=-apeSRByomo)
* [ミルカさんとフィボナッチ数列 結城浩 - pdf](https://www.hyuki.com/story/genfunc.pdf)

# 普通に実装

```python
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 で平方根の演算ができます。

```python:結果
1.0
1.0
2.0
3.0000000000000004
5.000000000000001
8.000000000000002
13.000000000000002
21.000000000000004
34.00000000000001
55.000000000000014
```

次は100までの結果を求めます。

```python
print(int(fib(100)))

# 354224848179263111168
```

100 番目の本当のフィボナッチ数は $354224848179261915075$ です。
残念ながら値が違ってしまいますね。
$354224848179263111168$
$354224848179261915075$

**たかが 100 番目くらいで正しい答えを得られないようでは"なってない"といわざるを得ないではない。**

# 計算方法を変える
では、どうすればいいのだろうとネット検索すると下記サイトが参考になりました。

>計算方針を考える
落ち着こう。確かに無理数は避けて通れない。しかし、我々は無理数どころか、"もっとスゴい"数を扱ったことがあるはずだ。そう、二乗すると -1 になっちゃう愛に溢れた"アレ"を含む数、複素数である。
複素数をどう扱ったか思い出そう。次のように表現したはずだ。
$a+bi$
虚数単位 $i$ はもはや無理数ですらない。それでも我々は扱ってきた。ならばたかが $\sqrt{5}$ 程度が扱えないわけがない。
というわけで、複素数の i を $\sqrt{5}$ に置き換えた形で今回の式を取り扱うことにしよう。
[フィボナッチ数列の一般項を計算する(※ただし有理数に限る)](http://labs.timedia.co.jp/2012/11/fibonacci-general-term-using-rational.html)

Haskellでは標準で Rational という有理数を扱う型があるが、Pythonでは標準 fractions を使うと分数(有理数)での計算ができる。

分子と分母をそれぞれ整数で指定する。

```python
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})$ である。

```math
\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*}
```

```haskell
fibMul :: FibNum -> FibNum -> FibNum
fibMul (a, b) (c, d) = (a*c + 5*b*d, a*d + b*c)
```

Pythonに移植

```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)
```

## 演算子を実装する : 累乗
```haskell
fibExp :: FibNum -> Integer -> FibNum
fibExp _ 0 = (1, 0)
fibExp a 1 = a
fibExp a n = fibMul a (fibExp a (n - 1))
```

Pythonに移植

```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乗の値を使用する。

```python
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
```

* [1.2.4 べき乗 - SCIP](https://sicp.iijlab.net/fulltext/x124.html)
* [フィボナッチ数列の任意の項を高速に計算する方法](http://ynomura.dip.jp/archives/2015/08/post_51.html)

普通のべき乗の計算なら pow 関数を使えばいいが、今回はタプル の Fraction を使っているので考え方を理解して作る必要がある。

```haskell
fibExp a n
    | even n    = ex
    | otherwise = ex `fibMul` a
    where ex = fibExp (fibMul a a) (n `div` 2)
```

Pythonに移植

```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));
```

## 演算子を実装する : 減算
```haskell
fibMin :: FibNum -> FibNum -> FibNum
fibMin (a, b) (c, d) = (a - c, b - d)
```

Pythonに移植

```python
def fibMin(fn1, fn2):
  a = fn1[0]
  b = fn1[1]
  c = fn2[0]
  d = fn2[1]
  return (a - c, b - d)
```

## メインを実装する
```haskell
fibMain :: Integer -> FibNum
fibMain n = (fibExp (1 % 2, 1 % 2) n `fibMin` fibExp (1 % 2, -1 % 2) n) `fibMul` (0, 1 % 5)
```

Pythonに移植

```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)
```

# 使ってみよう！

```python
for x in range(1, 10 + 1):
    print(fibMain(x))
```

```python:結果
(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$ になる。余分なコトを出力すると見づらいので、その辺は削ってしまおう。

```python
def fib(n):
  return fibMain(n)[0].numerator

for n in range(1, 10 + 1):
  print(fib(n))
```

表示がすっきりした。

```python:結果
1
1
2
3
5
8
13
21
34
55
```

## 実力を確認する
先程違っていた 100 番目の結果

```python
print(fib(100))

# 354224848179261915075
```

100 番目のフィボナッチ数 $354224848179261915075$ で一致しました。
やったね。

10,000 番目だっていけると思ったらエラーだ。

```python
print(fib(10000))

RecursionError: maximum recursion depth exceeded in comparison
```

どうやら再帰の上限を超えてしまうようだ。
[Python: maximum recursion depth exceeded while calling a Python object - stackoverflow](https://stackoverflow.com/questions/6809402/python-maximum-recursion-depth-exceeded-while-calling-a-python-object)
上限を追加すればいいようだ。

```python
import sys
sys.setrecursionlimit(15000)

print(fib(10000))

# 2090桁ある
# 3364476487643178326662161200510754331030214846068006390656476997468008144216666236815559551363373402...(中略)...366875
```

# 完成品
2018/12/30にfibExpを高速化版に変更

```python
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)でも出来そうな気がしますね。

# 複素数版
出来そうな気がしたのでやってみました。

```python
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番目もいけるかと思ったのですが、残念ながら結果が一致しませんでした。