トリボナッチ数列の一般項を線形代数とPythonで求める

(スマホからだと数式が表示できない場合があるようなので、その場合はPCなどからご覧ください。申し訳ありません。)

# トリボナッチ数列とは？

フィボナッチ数列は、直前の2つの項を足して次の項を作る数列です。式で表すと、

```math
\begin{cases}
a_{n+2} = a_{n+1} + a_{n}\\ a_0 = 0\\a_1 = 1
\end{cases}
```
となります。これを拡張して直前の3つの項を足した数を次の項にしたものがトリボナッチ数列です。式で表すと、

```math
\begin{cases}
a_{n+3} = a_{n+2} + a_{n+1} + a_{n}\\ a_0 = 0\\a_1 = 0\\a_2 = 1
\end{cases}
```
となります。少し計算してみると、

```math
a_3 = 1+0+0=1\\
a_4 = 1+1+0=2\\
a_5 = 2+1+1=4\\
a_6 = 4+2+1=7\\
\vdots
```
となります。

# 前準備

漸化式を満たす実数の数列のなす空間を$V$とします。数列$\{ a_n \}$の最初の3項 $a_0,a_1,a_2$ が与えられているので、$n\geq3$において数列 $\{ a_n \}$ は一意に定まります。

```math
\boldsymbol{u}=\{1,0,0,1,\dots \},\ \boldsymbol{v}=\{0,1,0,1,\dots \},\ \boldsymbol{w}=\{0,0,1,1,\dots\}
```

とします。$c_1,c_2,c_3$に対して、$c_1\boldsymbol{u}+c_2\boldsymbol{v}+c_3\boldsymbol{w}=\boldsymbol{o}$が成り立ったとすると、

```math
c_1\{1,0,0,\dots \}+c_2\{0,1,0,\dots \}+c_3\{0,0,1,\dots \}=\{c_1,c_2,c_3,\dots \}=\{0,0,0,\dots \}
```
となるので、$c_1=c_2=c_3=0$です。よって、$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$は一次独立とわかります。


次に、

```math
\boldsymbol{a}= \{ a_0,a_1,a_2,\dots \}
```

を$V$の勝手な元とすると、


```math
\begin{align}
\boldsymbol{a}&=\{ a_0,0,0,\dots \}+\{0,a_1,0,\dots \} +\{0,0,a_2,\dots \}  \\
 &=a_0\{1,0,0,\dots \}+a_1\{0,1,0,\dots \}+a_2\{0,0,1,\dots \} \\
&=a_0\boldsymbol{u}+a_1\boldsymbol{v}+a_2\boldsymbol{w}
\end{align}
```

となり、$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$ の一次結合で表せます。よって、$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$ は$V$を生成します。


以上より、$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$ は一次独立かつ$V$を生成するので、$V$の基底となります。

ここで、写像 $f:V\rightarrow V$を
 
```math
f(\{ a_n\}_{n=0}^{\infty})=\{ a_n\}_{n=1}^{\infty}
```

とします。$\boldsymbol{a}=\{ a_0,a_1,a_2,\dots \} \in V$のとき、$f(\boldsymbol{a})=\{ a_1,a_2,a_3,\dots \}$も漸化式を満たすので、$f(\boldsymbol{a})\in V$です。

```math
\boldsymbol{a}=\{ a_n\}_{n=0}^{\infty}\in V \\
\boldsymbol{b}=\{ b_n\}_{n=0}^{\infty}\in V \\
```
とすると、定数 $c$ に対して、

```math
f(\boldsymbol{a}+\boldsymbol{b})=f(\{ a_n+b_n\}_{n=0}^{\infty})=\{ a_n+b_n\}_{n=1}^{\infty}=\{ a_n\}_{n=1}^{\infty}+\{ b_n\}_{n=1}^{\infty}=f(\boldsymbol{a})+f(\boldsymbol{b})\\
f(c\boldsymbol{a})=f(c\{ a_n\}_{n=0}^{\infty})=c\{ a_n\}_{n=1}^{\infty}=cf(\boldsymbol{a})
```

が成り立つので、$f$は$V$の線形変換です。

# 漸化式を行列で表現

$\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$に関して、

```math
\begin{align}
f(\boldsymbol{u})&=\{0,0,1,\dots \}=\boldsymbol{w}\\
f(\boldsymbol{v})&=\{1,0,1,\dots \}=\boldsymbol{u}+\boldsymbol{w}\\
f(\boldsymbol{w})&=\{0,1,1,\dots \}=\boldsymbol{v}+\boldsymbol{w}
\end{align}
```

なので、

```math
[f(\boldsymbol{u})\quad f(\boldsymbol{v})\quad f(\boldsymbol{w})]=
[\boldsymbol{u}\ \boldsymbol{v}\ \boldsymbol{w}]
\begin{bmatrix}
0 & 1 & 0 \\
0 & 0 & 1 \\
1 & 1 & 1
\end{bmatrix}
```
と表せます。よって、$f$の基底 $\boldsymbol{u},\boldsymbol{v},\boldsymbol{w}$に関する表現行列は、

```math
\begin{bmatrix}
0 & 1 & 0 \\
0 & 0 & 1 \\
1 & 1 & 1
\end{bmatrix}
```

です。この表現行列を$A$とします。


# 固有値と公比の関係

```math
\boldsymbol{p}=\{ r^{n} \}_{n=0}^{\infty}=\{1,r,r^2,\dots \}
```
が$V$に属したとすると、

```math
f(\boldsymbol{p})=f(\{ r^{n} \}_{n=0}^{\infty})=\{ r^{n} \}_{n=1}^{\infty}=\{ r^{n+1} \}_{n=0}^{\infty}=r\{ r^{n} \}_{n=0}^{\infty}=r\boldsymbol{p}
```

より、$\boldsymbol{p}$は固有値$r$の固有ベクトルになります。逆に、上の式から$f$の固有値 $r$の固有ベクトルは、公比$r$の等比数列になることも分かります。よって、公比と固有値は等しいとわかります。

# 固有値を求める

$f$の固有値を求めます。$A$の固有多項式は、$I$を単位行列として、

```math
\begin{align}
\varphi_A(\lambda)&=|A-\lambda I|\\
&=\begin{vmatrix}
0-\lambda & 1 & 0 \\
0 & 0-\lambda & 1 \\
1 & 1 & 1-\lambda
\end{vmatrix}\\
&=-\lambda^3+\lambda^2+\lambda+1
\end{align}
```
となります。しかし、$\varphi_A(\lambda)=0$は簡単には解くことはできません。3つの解をそれぞれ $\alpha,\beta,\gamma$とします。$\alpha,\beta,\gamma$は、あとの章で求めることにします。解 $\alpha,\beta,\gamma$に対応する固有ベクトルである等比数列の

```math
\boldsymbol{a}=\{ \alpha^{n} \}_{n=0}^{\infty}\\
 \boldsymbol{b}=\{ \beta^{n} \}_{n=0}^{\infty}\\
 \boldsymbol{c}=\{ \gamma^{n} \}_{n=0}^{\infty}
```
は$V$に属します。これらは $\alpha,\beta,\gamma$が全て異なるとすると、$f$の相異なる固有値に対応する固有ベクトルなので一次独立です。また、$\dim V=3$です。よって、$\boldsymbol{a},\boldsymbol{b},\boldsymbol{c}$は$V$の基底となります。


# 係数の決定
以上より、ある $c_1,c_2,c_3$が存在して、

```math
\boldsymbol{a_n}=c_1\boldsymbol{a}+c_2\boldsymbol{b}+c_3\boldsymbol{c}
```
と表せるので、

```math
\begin{align}
\boldsymbol{a_n}&=c_1\boldsymbol{a}+c_2\boldsymbol{b}+c_3\boldsymbol{c}\\
\Leftrightarrow \{a_0,a_1,a_2,\dots \}&=c_1\{ \alpha^{n} \}_{n=0}^{\infty}+c_2\{ \beta^{n} \}_{n=0}^{\infty}+c_3\{ \gamma^{n} \}_{n=0}^{\infty}\\
\Leftrightarrow \{0,0,1,\dots \}&=c_1\{ 1,\alpha,\alpha^2,\dots \}+c_2\{1,\beta,\beta^2,\dots \}+c_3\{ 1,\gamma,\gamma^2,\dots\} \\
\Leftrightarrow \{0,0,1,\dots \}&=\{ c_1+c_2+c_3,c_1\alpha+c_2\beta+c_3\gamma,c_1\alpha^2+c_2\beta^2+c_3\gamma^2, \dots\} 
\end{align}
```
これを行列で表現すると、

```math
\begin{bmatrix}
1 & 1 & 1 \\
\alpha & \beta & \gamma \\
\alpha^2 & \beta^2 & \gamma^2
\end{bmatrix}
\begin{bmatrix}
c_1 \\ c_2 \\ c_3
\end{bmatrix}
=
\begin{bmatrix}
0 \\ 0 \\ 1
\end{bmatrix}
```

です。左の行列の行列式は、

```math
\begin{align}
\begin{vmatrix}
1 & 1 & 1 \\
\alpha & \beta & \gamma \\
\alpha^2 & \beta^2 & \gamma^2
\end{vmatrix}
&=
\begin{vmatrix}
1 & 0 & 0 \\
\alpha & \beta-\alpha & \gamma-\alpha \\
\alpha^2 & \beta^2-\alpha^2 & \gamma^2-\alpha^2
\end{vmatrix}\\
&=
\begin{vmatrix}
\beta-\alpha & \gamma-\alpha \\
\beta^2-\alpha^2 & \gamma^2-\alpha^2
\end{vmatrix}\\
&=\begin{vmatrix}
\beta-\alpha & \gamma-\alpha \\
\beta^2-\alpha^2 - (\beta-\alpha)(\beta + \alpha) & \gamma^2-\alpha^2 - (\gamma-\alpha)(\beta + \alpha)
\end{vmatrix}\\
&=\begin{vmatrix}
\beta-\alpha & \gamma-\alpha \\
0 & (\gamma-\alpha)(\gamma+\alpha) - (\gamma-\alpha)(\beta + \alpha)
\end{vmatrix}\\
&=(\beta-\alpha)(\gamma-\alpha)(\gamma - \beta)\\
&=(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
\end{align}
```



となります。これはヴァンデルモンド行列式としても知られています。クラメルの公式を用いると、

```math
\begin{align}
c_1 &= \frac{
\begin{vmatrix}
0 & 1 & 1 \\
0 & \beta & \gamma \\
1 & \beta^2 & \gamma^2
\end{vmatrix}
}{
\begin{vmatrix}
1 & 1 & 1 \\
\alpha & \beta & \gamma \\
\alpha^2 & \beta^2 & \gamma^2
\end{vmatrix}
}\\
&=
\frac{
-\begin{vmatrix}
1 & \beta^2 & \gamma^2 \\
0 & \beta & \gamma \\
0 & 1 & 1 
\end{vmatrix}
}{
(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}\\
&=
\frac{
-\begin{vmatrix}
\beta & \gamma\\
1 & 1 
\end{vmatrix}
}{(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}\\
&=\frac{
\gamma-\beta
}{(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}\\
&=\frac{1}{(\alpha-\beta)(\alpha-\gamma)
}
\end{align}
```


と求められます。同様に、

```math
c_2 = \frac{
\alpha-\gamma
}{(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}=\frac{1
}{(\beta-\alpha)(\beta-\gamma)
}\\
c_3 = \frac{
\beta-\alpha 
}{(\alpha-\beta)(\beta-\gamma)(\gamma-\alpha)
}
=\frac{1
}{(\gamma-\alpha)(\gamma-\beta)
}
```
と求められます。

以上より、

```math
\begin{align}
\boldsymbol{a_n} &= c_1 \boldsymbol{a} + c_2 \boldsymbol{b} + c_3 \boldsymbol{c}\\
&=\left\{ \frac{\alpha^n}{(\alpha-\beta)(\alpha-\gamma)} + \frac{\beta^n}{(\beta-\alpha)(\beta-\gamma)} + \frac{\gamma^n}{(\gamma-\alpha)(\gamma-\beta)}\right\}_{n=0}^{\infty}
\end{align}
```
なので、

```math
a_n = \frac{\alpha^n}{(\alpha-\beta)(\alpha-\gamma)} + \frac{\beta^n}{(\beta-\alpha)(\beta-\gamma)} + \frac{\gamma^n}{(\gamma-\alpha)(\gamma-\beta)}
```

が導けました。

# 具体的な数値の計算

以上までで、一般項の形はわかりましたが、具体的な $\alpha,\beta,\gamma$の値はわかっていません。3次方程式の解の公式を使って求める方法もありますが、今回はPythonで求めてみましょう。$a[i]$をトリボナッチ数列の定義通りに計算したもの、$b[i]$を求めた一般項で計算したものとします。sympyを使って代数的に解く方法もありますが、今回はnumpyを使って数値的に求めます。

#### ソースコード

```python
import numpy as np
import warnings 
warnings.filterwarnings('ignore')       # 複素数の計算時にWarningが発生するので消しておく

alpha, beta, gamma = np.roots([1,-1,-1,-1])     # solve 1*x^3 + (-1)*x^2+ (-1)*x + (-1) = 0
print("alpha =", alpha)
print("beta =", beta)
print("gamma =", gamma)

a = np.zeros(101)
b = np.zeros(101)

a[2] = 1

print(" i ","           a[i]           ","          b[i]           ", "  a[i]-b[i]    ")

for i in range(101):
    if i > 2:
        a[i] = a[i-1] + a[i-2] + a[i-3]
    b[i] = np.power(alpha,i)/((alpha-beta)*(alpha-gamma)) + np.power(beta,i)/((beta-alpha)*(beta-gamma)) + np.power(gamma,i)/((gamma-alpha)*(gamma-beta))
    print('{:3}'.format(i), '{:26}'.format(int(a[i])), '{:26}'.format(int(round(abs(b[i])))), '{:12}'.format(int(a[i]) - int(round(abs(b[i])))))
```


#### 出力結果

```
alpha = (1.839286755214161+0j)
beta = (-0.41964337760708065+0.6062907292071997j)
gamma = (-0.41964337760708065-0.6062907292071997j)
 i             a[i]                      b[i]              a[i]-b[i]    
  0                          0                          0            0
  1                          0                          0            0
  2                          1                          1            0
  3                          1                          1            0
  4                          2                          2            0
  5                          4                          4            0
  6                          7                          7            0
  7                         13                         13            0
  8                         24                         24            0
  9                         44                         44            0
 10                         81                         81            0
 11                        149                        149            0
 12                        274                        274            0
 13                        504                        504            0
 14                        927                        927            0
 15                       1705                       1705            0
 16                       3136                       3136            0
 17                       5768                       5768            0
 18                      10609                      10609            0
 19                      19513                      19513            0
 20                      35890                      35890            0
 21                      66012                      66012            0
 22                     121415                     121415            0
 23                     223317                     223317            0
 24                     410744                     410744            0
 25                     755476                     755476            0
 26                    1389537                    1389537            0
 27                    2555757                    2555757            0
 28                    4700770                    4700770            0
 29                    8646064                    8646064            0
 30                   15902591                   15902591            0
 31                   29249425                   29249425            0
 32                   53798080                   53798080            0
 33                   98950096                   98950096            0
 34                  181997601                  181997601            0
 35                  334745777                  334745777            0
 36                  615693474                  615693474            0
 37                 1132436852                 1132436852            0
 38                 2082876103                 2082876103            0
 39                 3831006429                 3831006429            0
 40                 7046319384                 7046319384            0
 41                12960201916                12960201916            0
 42                23837527729                23837527729            0
 43                43844049029                43844049029            0
 44                80641778674                80641778674            0
 45               148323355432               148323355432            0
 46               272809183135               272809183135            0
 47               501774317241               501774317241            0
 48               922906855808               922906855808            0
 49              1697490356184              1697490356184            0
 50              3122171529233              3122171529233            0
 51              5742568741225              5742568741225            0
 52             10562230626642             10562230626642            0
 53             19426970897100             19426970897100            0
 54             35731770264967             35731770264967            0
 55             65720971788709             65720971788709            0
 56            120879712950776            120879712950775            1
 57            222332455004452            222332455004451            1
 58            408933139743937            408933139743935            2
 59            752145307699165            752145307699160            5
 60           1383410902447554           1383410902447546            8
 61           2544489349890656           2544489349890641           15
 62           4680045560037375           4680045560037345           30
 63           8607945812375585           8607945812375530           55
 64          15832480722303616          15832480722303512          104
 65          29120472094716576          29120472094716384          192
 66          53560898629395776          53560898629395416          360
 67          98513851446415968          98513851446415296          672
 68         181195222170528320         181195222170527072         1248
 69         333269972246340096         333269972246337792         2304
 70         612979045863284352         612979045863280000         4352
 71        1127444240280152832        1127444240280144512         8320
 72        2073693258389777408        2073693258389762048        15360
 73        3814116544533214720        3814116544533186560        28160
 74        7015254043203144704        7015254043203092480        52224
 75       12903063846126137344       12903063846126039040        98304
 76       23732434433862496256       23732434433862311936       184320
 77       43650752323191775232       43650752323191439360       335872
 78       80286250603180408832       80286250603179769856       638976
 79      147669437360234692608      147669437360233480192      1212416
 80      271606440286606852096      271606440286604623872      2228224
 81      499562128250021937152      499562128250017808384      4128768
 82      918838005896863416320      918838005896855945216      7471104
 83     1690006574433492271104     1690006574433477853184     14417920
 84     3108406708580377427968     3108406708580351213568     26214400
 85     5717251288910732984320     5717251288910684749824     48234496
 86    10515664571924601634816    10515664571924511457280     90177536
 87    19341322569415712047104    19341322569415540080640    171966464
 88    35574238430251047714816    35574238430250737336320    310378496
 89    65431225571591361396736    65431225571590774194176    587202560
 90   120346786571258121158656   120346786571257030639616   1090519040
 91   221352250573100547047424   221352250573098500227072   2046820352
 92   407130262715950037991424   407130262715946279895040   3758096384
 93   748829299860308689420288   748829299860301844316160   6845104128
 94  1377311813149359375122432  1377311813149346221785088  13153337344
 95  2533271375725618236751872  2533271375725593540689920  24696061952
 96  4659412488735286167076864  4659412488735240533049344  45634027520
 97  8569995677610263510515712  8569995677610179758653440  83751862272
 98 15762679542071167914344448 15762679542071011148038144 156766306304
 99 28992087708416715444453376 28992087708416427681644544 287762808832
100 53324762928098150090539008 53324762928097265327276032 884763262976
```
このように、$i$ が小さいときは一致しています。$i$ が大きいときは浮動小数点数による誤差が出てきていますが、数値に比べて誤差が小さいので一般項は正しいと考えられます。

#### 計算量

Pythonでは累乗に用いるpow関数は繰り返し二乗法で実装されているため、ある数の$N$乗は$O(\log N)$ で高速に求められます。そのため、$a[N]$を 求めるのに必要な計算量は $O(N)$ ですが、$b[N]$ は $O(\log N)$ で求めることができます。計算機上では、$a[i]$では値が正確に求められるが計算量は大きい、$b[i]$では高速に計算できるが誤差が発生する、というそれぞれのメリットとデメリットがあります。

# トリボナッチ数列をさらに拡張
[拡張フィボナッチ数列(k-ボナッチ数列 : 直前のk項の和)の一般項を線形代数とPythonで求める](https://qiita.com/magolors/items/359b1da94830a71fec08)という記事を書きました。こちらもよければご覧ください。

# 参考文献
[こちらのpdf](http://www.math.twcu.ac.jp/~yamauchi/text/2016/linear_iii/2016_linear_iii_ans.pdf)を参考にさせていただきました。ありがとうございました。