拡張されたフィボナッチ数列をPythonで実装してみるPart3 (初期値の変更)

# はじめに
フィボナッチ数列をPythonで実装するシリーズ記事のPart3です。Part0からPart2を読了したことを前提としているのでまだの方は是非お読みください。

- Part0:[Pythonによるフィボナッチ数列の色々な求め方](https://qiita.com/dovedove/items/3456c4f317a5c680f437)
- Part1:[負の数番を考える](https://qiita.com/dovedove/items/5146718b7500f882c211)
- Part2:[項数の一般化](https://qiita.com/dovedove/items/fcdd945ef7dcd8948530)
- Part3:初期値の一般化 (イマココ)

ではPart3の今回はフィボナッチ数列の初期値をPython3で変更していきましょう。


# 初期値の変更
フィボナッチ数列は二つの初期条件(初項$F_0=0$,第二項$F_1=1$)を持つ漸化式で表される数列であると述べてきました。しかし、例えば$F_0=8$,$F_1=15$という初期条件を持つ数列の数$F_n$というのも考えることができそうです。
以下ではこのような初期値の変更について考えていきます。

## 初期値変更の例
### 0-filと1-fil
まずは簡単な例から考えていきましょう。初期値を0で埋めたものを0-fil、1で埋めたものを1-filといいます。
例えば0-filのフィボナッチ数列の初期値はみなさんご存知の$F_0=0$,$F_1=1$です。
1-filのフィボナッチ数列は$F_0=1$,$F_1=1$の数列です。
一般化すると、0-filのn-ボナッチ数列の初期値は$N_0=N_1=\cdots=N_{n-2}=0$ , $N_{n-1}=1$であり、1-filの場合の初期値は$N_0=N_1=\cdots=N_{n-1}=1$です。

ここでは1-filのフィボナッチ数列を考えてみます。ほぼ普通のフィボナッチ数列と変わらないので再帰法の紹介だけにとどめます。

```.py
def Fib(n):    
    if n in [0,1]:
        return 1
    else:
        return Fib(n-1)+Fib(n-2)

print([Fib(n) for n in range(15)])
# [1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377, 610]
```


### リュカ数
フィボナッチ数列の初期値を$F_0=2$,$F_1=1$に置き換えた数列の項をリュカ数といいます。
ではリュカ数の列[^4] $2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521, 843, 1364, 2207, \cdots$ について考えていきましょう。

[^4]: 「リュカ数列」というと数学的には別の意味になってしまうので、あえて「リュカ数の列」と表記しています。

#### 再帰法
まずは再帰法です。初期条件を変更するだけです。

```.py
def Lucas(n):  
    if n == 0:
        return 2
    if n == 1:
        return 1
    else:
        return Lucas(n-1)+Lucas(n-2)

print([Lucas(n) for n in range(15)])
# [2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521, 843]
```

#### 一般項
次に一般項を求めてみましょう。sympyの漸化式を解く機能は優秀なので(トリボナッチ数列の一般項は導出できませんでしたが)、初期値を辞書型で自由に与えられます。

```.py
import sympy as sym

x = sym.symbol('x')
f = sym.simplify("f(x)") 
s = sym.simplify("f(x)-f(x-1)-f(x-2)")
ini = sym.simplify({0:2,1:1})
Lucas_term = sym.rsolve(s,f,ini) #初期条件iniの漸化式sをfについて解く

print(Lucas_term) 
# (1/2 + sqrt(5)/2)**x + (-sqrt(5)/2 + 1/2)**x
```

求めた`Lucas_term`を使ってリュカ数を求めてみましょう。

```.py
def Lucas(n):
    result=Lucas_term.subs(x,n) #xにnを代入
    result=result.evalf(int(3*n/14)+1) #精度の指定(リュカ数もフィボナッチ数もnに対する桁数はほぼ同じなので変えなくてよい)
    return int(result)

print([Lucas(n) for n in range(15)])
# [2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199, 322, 521, 843]
```


## 任意の整数初期値
ここからが本題です。フィボナッチ数列の漸化式は隣接三項間漸化式というものの一種なのですが、この漸化式は二つの(異なる)初期値を与えると一意に定まります。つまり、任意に初期値を与えることができます。
では、任意の初期値を与えるフィボナッチ数列(の類似数列)について考えていきましょう。

### 実装のルール
実装のルールを考えます。実装する関数は`Fib(n,ini)`で、初期値`ini`のフィボナッチ数$F_n$を求めます。`ini`は辞書型で`{番号1:値1,番号2:値2}`のように与えます。つまり$F_8=21,F_{9}=34$が与えられている際の$F_{11}$を求めたい場合は`Fib(11,{8:21,9:34})`とします。

### 任意の連続する整数番の整数初期値
まずは少し簡単なものを考えます。

#### 再帰法
$F_8=21,F_9=34$ が与えられた時を例にとって説明していきます(ちなみにこれは普通のフィボナッチ数列となります)。
このとき、$F_6$を求めようと$F_n=F_{n-1}+F_{n-2}$を適用していくと、`Fib(6)`は再帰的に`Fib(5)`と`Fib(4)`を呼び出します。`Fib(5)`は`Fib(4)`と`Fib(3)`を呼び出します。`Fib(4)`は...と一向に初期値にたどり着きません。無限再帰に陥ってしまいます。
「あれ？でもこれなんか見たことあるぞ？」と思いませんでしたか？
そう、[負の数番のフィボナッチ数列を考えた時に見た光景](https://qiita.com/dovedove/items/5146718b7500f882c211#やっぱり再帰で)ですね。つまり、同様の解決法が使えるということです。
初期値$F_k,F_{k+1}$が与えられている数列の$F_n$を求めるとき、 $n>k+1$ ならば $F_n=F_{n-1}+F_{n-2}$ を適用し、 $n<k$ ならば$F_{n}=F_{n+2}-F_{n+1}$を適用するようにします。
ただし、引数で与えられる初期値`ini`が必ず$k,k+1$の順番で与えられるわけではない(そもそも辞書型は順番を保証しない)ので、どちらが$k+1$でどちらが$k$なのか判定する必要があります。`max(ini.keys())`で$k+1$、`min(ini.keys())`で$k$を返します。

```.py
def Fib(n,ini):  
    if n in ini:
        return ini[n]
    elif n > max(ini.keys()):
        return Fib(n-1,ini)+Fib(n-2,ini)
    elif n < min(ini.keys()):
        return Fib(n+2,ini)-Fib(n+1,ini)

print([Fib(n, {8:21,9:34}) for n in range(10)])
# [0, 1, 1, 2, 3, 5, 8, 13, 21, 34]
```
正しく出力されているようです。

#### 一般項
前述したようにsympyの漸化式を解く機能は初期値を辞書で自由に与えられます。これを使えば任意初期値のフィボナッチ数列が実装できるというわけです。

```.py
import sympy as sym

def Fib(n,ini):
    x=sym.symbols('x')
    f = sym.simplify("f(x)") 
    s = sym.simplify("f(x)-f(x-1)-f(x-2)")
    ini = sym.simplify(ini)
    Fib_term=sym.rsolve(s,f,ini) #一般項の式
    result=Fib_term.subs(x,n) #xにnを代入
    result=sym.simplify(result) #式の整理
    return result

print([Fib(n,{8:21,9:34}) for n in range(12)]) #フィボナッチ数列
# [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89]

print([Fib(n,{0:2,1:1}) for n in range(12)]) #リュカ数の列
# [2, 1, 3, 4, 7, 11, 18, 29, 47, 76, 123, 199]
```
いちいち一般項の導出をしているため、速度は遅いですが、しっかりと出力できています。

精度の検証もしておきましょう。
$F_{100}=0,F_{101}=1$のときの$F_0$を考えてみましょう。この数列は$F_{100}$を初項とするフィボナッチ数列と同等とみなせるので、$F_0$は負数番フィボナッチ数$F_{-100}$と一致します。すなわち、$F_0=-354224848179261915075$です。
では`Fib(0,{100:0,101:1})`を出力してみます。

```.py
print(Fib(0,{100:0,101:1}))
# -354224848179261915075
```
確かな精度で出力できています。

### 任意の整数番の実数初期値
さきほどは「連続する整数番」「初期値は整数」という縛りを与えて考えました。
次は連続しない(非負の)整数番の実数初期値が与えられている場合も考えましょう。
#### 整数初期値
まずは整数初期値に絞って話を進めます。
##### 再帰法
再帰法での実装は一筋縄ではいきません。


任意の連続する整数番の整数初期値の項目で実装した関数が汎用的に使えるかどうかを確かめるために、試しに初期値$F_8=21,F_{11}=89$が与えられている際の$F_{15}$を求めてみましょう(この数列は普通のフィボナッチ数列です)。

```.py
def cFib(n,ini):  #後々のために名前を変えています
    if n in ini:
        return ini[n]
    elif n > max(ini.keys()):
        return cFib(n-1,ini)+cFib(n-2,ini)
    elif n < min(ini.keys()):
        return cFib(n+2,ini)-cFib(n+1,ini)

print(Fib(15, {8:21,11:89}))
# TypeError: unsupported operand type(s) for +: 'int' and 'NoneType'
```

エラーが出てしまいます。この実装だと、不連続な整数番初期値$F_i,F_j$が与えられているときの$F_n(i<n<j)$を求める際に答えを出すことができないことが原因です。

再帰法で解を出すときは、初期値が連続する整数番でないと都合が悪いのです。不連続な整数番の初期値に対応するには、我々が良く知る「方程式」の考えを導入します。
連続する整数番の初期値が与えられていないなら、自分で未知数として置いてしまえばいいのです。
プロセスの例を以下に示します。

| $F_{8}$ | $F_{9}$ | $F_{10}$ | $F_{11}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $21$ |  |   | $89$  | 

| $F_{8}$ | $F_{9}$ | $F_{10}$ | $F_{11}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $21$ | $x$ |   | $89$  | 

| $F_{8}$ | $F_{9}$ | $F_{10}$ | $F_{11}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $21$ | $x$ | $x+21$  | $89$  | 

| $F_{8}$ | $F_{9}$ | $F_{10}$ | $F_{11}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $21$ | $x$ | $x+21$  | $2x+21=89$  | 

$2x+21=89 ~\therefore~ x=34$

これで、不連続な整数番の初期値から、連続する整数番の初期値を得ることができました。この初期値を使えば再帰法で解を出せます。

さて、このプロセスを実装していきましょう。

初期値`ini`で、$F_a,F_b(a<b)$の値が与えられているとします。

| $F_{a}$ | $\cdots$ | $F_{b}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $ini[a]$ | $\cdots$ | $ini[b]$  | 

$F_{a+1}=x$とおき、さきほどの連続整数番初期値用の関数`cFib`を用いて$F_{b}$を求めます。

| $F_{a}$ |$F_{a+1}$ |$\cdots$ | $F_{b}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $ini[a]$ |$x$| $\cdots$ | $cFib(b,\\{a:ini[a]~,~a+1:x\\})$  | 

一次方程式$cFib(b,\\{a:ini[a]~,~a+1:x\\})-ini[b]=0$を解きます。解を$x=k$とすると、任意の$n$に対して、$F_n=cFib(n,\\{a:ini[a]~,~a+1:k\\})$が成り立ちます。

コード例を以下に示します。

```.py
import sympy as sym

def Fib(n,ini):    
    
    #関数cFibの定義
    def cFib(n,ini):  
        if n in ini:
            return ini[n]
        elif n > max(ini.keys()):
            return cFib(n-1,ini)+cFib(n-2,ini)
        elif n < min(ini.keys()):
            return cFib(n+2,ini)-cFib(n+1,ini)
    
    x = sym.symbols('x')
    a = min(ini.keys())
    b = max(ini.keys())
    
    #単にcFibで解ける場合はcFibで解く
    if n in ini or abs(a-b) == 1:
        return cFib(n, ini)
    
    continuousIni = {a:ini[a], a+1:x} #未知数を用いて連続する整数番の初期値を作成
    equation = cFib(b, continuousIni) - ini[b] #方程式の定義(sympyでは(数式)=0となるように数式を定義すると方程式が解ける)
    
    continuousIni[a+1] = sym.solve(equation)[0] #xを解kに置き換え
    return cFib(n, continuousIni)
   
 
print([Fib(n, {8:21,11:89}) for n in range(15)])
# [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377] 
```

##### 一般項
一般項による方法での実装では、初期値を辞書で自由に与えられるので、不連続な整数番の整数初期値についても同様の実装でOKです。

```.py
import sympy as sym

def Fib(n,ini):
    x=sym.symbols('x')
    f = sym.simplify("f(x)") 
    s = sym.simplify("f(x)-f(x-1)-f(x-2)")
    ini = sym.simplify(ini)
    Fib_term=sym.rsolve(s,f,ini) #一般項の式
    result=Fib_term.subs(x,n) #xにnを代入
    result=sym.simplify(result) #式の整理
    return result

print([Fib(n, {8:21,11:89}) for n in range(15)])
# [0, 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144, 233, 377] 
```

##### 有理数の扱い
連続する整数番の整数初期値を持つフィボナッチ数列(の類似数列)においては、数列のどの数をとっても整数です。これは、整数同士の和が必ず整数であることから明らかです。
しかし、不連続整数番の整数初期値を持つフィボナッチ数列(の類似数列)においては、必ず整数であるわけではありません(非整数の有理数が解になることがあり得ます)。
実際に、例えば$F_8=1,F_{11}=8$のときの$F_9$は分数となります。

| $F_{8}$ | $F_{9}$ | $F_{10}$ | $F_{11}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $1$ |  |   | $8$  | 

| $F_{8}$ | $F_{9}$ | $F_{10}$ | $F_{11}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $1$ | $x$ |   | $8$  | 

| $F_{8}$ | $F_{9}$ | $F_{10}$ | $F_{11}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $1$ | $x$ | $x+1$  | $8$  | 

| $F_{8}$ | $F_{9}$ | $F_{10}$ | $F_{11}$ | 
|:-:|:-:|:-:|:-:|:-:|
|  $1$ | $x$ | $x+1$  | $2x+1=8$  | 
$2x+1=8 ~\therefore~ x=\frac{7}{2}$

このとき、上で示した実装法では両方、以下のように分数のまま表示されます。

```.py
print(Fib(9,{8:1,11:8}))
# 7/2
```
関数内で`float(result)`のようにすれば、`3.5`というような出力が得られますが、ここでは有理数のままにしておきます。

#### 実数に拡張
初期値に実数を用いることを考えます。
しかし、実装を変えなくても実数初期値に対応しています。
実際に、実数初期値についても以下のように出力を得ることができます。

```.py
print(Fib(9,{8:5.8,11:sym.pi})) #sym.piは円周率π(無理数)
# -1.32920367320510 (再帰法での出力)
# -3826.19572436128*sqrt(5) + 8554.30453121767 (一般項による方法での出力)
```
これでもいいですが、一般項による方法での出力が多項式となってしまっているので、変えておきましょう。
ついでに精度も任意で与えられるようにします。`Fib(n,ini,accu)`というような`accu`で精度(有効桁数)を指定する関数をつくります。

```.py
import sympy as sym

def Fib(n,ini,accu):
    x=sym.symbols('x')
    f = sym.simplify("f(x)") 
    s = sym.simplify("f(x)-f(x-1)-f(x-2)")
    ini = sym.simplify(ini)
    Fib_term=sym.rsolve(s,f,ini) #一般項の式
    result=Fib_term.subs(x,n) #xにnを代入
    result=result.evalf(accu) #有効桁数の指定
    return result

print(Fib(9,{8:5.8,11:sym.pi},20))
# -1.3292036732051035798
```


#### 複素数初期値について
考えることもできますし、上記の実装で出力もできますが、あまり面白くはありません。

```.py
print(Fib(9,{8:5+sym.I,11:3})) #sym.Iは虚数単位i
# -1 - I/2 (再帰法および精度指定なしの一般項による方法での出力)
```


# おわりに
今回はフィボナッチ数の初期値を任意に与えられるようにしました。

誤字脱字や誤った記述、更に良いアルゴリズム、気づいたこと、分かりにくい箇所、感想などがあればぜひコメントをください。


# 参考文献
[フィボナッチ数 - Wikipedia](https://ja.wikipedia.org/wiki/フィボナッチ数)