(スマホからだと数式が表示できない場合があるようなので、その場合はPCなどからご覧ください。申し訳ありません。)
拡張フィボナッチ数列とは？
フィボナッチ数列は、直前の2つの項を足して次の項を作る数列です。式で表すと、
\begin{cases}
\begin{align}
a_{n+2} &= a_{n+1} + a_{n}\\
a_0 &= 0\\
a_1 &= 1
\end{align}
\end{cases}
となります。これを拡張して直前の$k$個の項を足した数を次の項にしたものが拡張フィボナッチ数列です。名前を探しましたが、見つからなかったのでここでは $k$ -ボナッチ数列と呼ぶことにします。式で表すと、
\begin{cases}
\begin{align}
a_{n+k} &= a_{n+k-1}+a_{n+k-2}+\dots+a_{n}\\
&=\displaystyle \sum_{i=n}^{n+k-1}a_i
\\ a_0 &= a_1 = \dots = a_{k-2} = 0\\a_{k-1} &= 1
\end{align}
\end{cases}
となります。前回の記事(トリボナッチ数列の一般項を線形代数とPythonで求める)で $k=3$ である、トリボナッチ数列を扱いました。今回はそれを拡張していきます。
前準備
(トリボナッチ数列の一般項を線形代数とPythonで求めるとほぼ同じですので、飛ばしていただいて構いません。)
漸化式を満たす実数の数列のなす空間を$V$とします。数列${ a_n }$の最初の$k$項 $a_0,a_1,\dots,a_{k-1}$ が与えられているので、$n\geq k$において数列 ${ a_n }$ は一意に定まります。
\begin{align}
\boldsymbol{v_0}&=\{\overbrace{1,0,0,\dots,0}^{k個},1,\dots \}\\
\boldsymbol{v_1}&=\{0,1,0,\dots,0,1,\dots \}\\
& \vdots \\
\boldsymbol{v_{k-1}}&=\{0,0,0,\dots,1,1,\dots\}
\end{align}
とします。$c_0,c_1,\dots,c_{k-1}$に対して、$c_0\boldsymbol{v_0}+c_1\boldsymbol{v_1}+\dots+c_{k-1}\boldsymbol{v_{k-1}}=\boldsymbol{0}$が成り立ったとすると、
c_0\{1,0,0,\dots \}+c_1\{0,1,0,\dots \}+\dots+c_{k-1}\{0,\dots,1,\dots \}\\=\{c_0,c_1,\dots,c_{k-1},\dots \}=\{0,0,\dots,0,\dots \}
となるので、$c_0=c_1=\dots =c_{k-1}=0$ です。よって、$\boldsymbol{v_0},\boldsymbol{v_1},\dots,\boldsymbol{v_{k-1}}$は一次独立とわかります。
次に、
\boldsymbol{a}= \{ a_0,a_1,\dots,a_{k-1},\dots \}
を$V$の勝手な元とすると、
\begin{align}
\boldsymbol{a}&=\{ a_0,0,0,\dots \}+\{0,a_1,0,\dots \} +\dots+\{0,\dots,0,a_{k-1},\dots \} \\
&=a_0\{1,0,0,\dots \}+a_1\{0,1,0,\dots \}+\dots+a_{k-1}\{0,\dots,0,1,\dots \} \\
&=a_0\boldsymbol{v_0}+a_1\boldsymbol{v_1}+\dots+a_{k-1}\boldsymbol{v_{k-1}}
\end{align}
となり、$\boldsymbol{v_0},\boldsymbol{v_1},\dots,\boldsymbol{v_{k-1}}$ の一次結合で表せます。よって、$\boldsymbol{v_0},\boldsymbol{v_1},\dots,\boldsymbol{v_{k-1}}$ は$V$を生成します。
以上より、$\boldsymbol{v_0},\boldsymbol{v_1},\dots,\boldsymbol{v_{k-1}}$ は一次独立かつ$V$を生成するので、$V$の基底となります。
ここで、写像 $f:V\rightarrow V$を
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$です。
\boldsymbol{a}=\{ a_n\}_{n=0}^{\infty}\in V \\
\boldsymbol{b}=\{ b_n\}_{n=0}^{\infty}\in V \\
とすると、定数 $c$ に対して、
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{v_0},\boldsymbol{v_1},\dots,\boldsymbol{v_{k-1}}$ に関して、
\begin{align}
f(\boldsymbol{v_0})&=\{0,0,\dots,0,1,\dots \}=\boldsymbol{v_{k-1}}\\
f(\boldsymbol{v_1})&=\{1,0,\dots,0,1,\dots \}=\boldsymbol{v_0}+\boldsymbol{v_{k-1}}\\
f(\boldsymbol{v_2})&=\{0,1,\dots,0,1,\dots \}=\boldsymbol{v_1}+\boldsymbol{v_{k-1}}\\
\vdots \\
f(\boldsymbol{v_{k-1}})&=\{0,0,\dots,1,1,\dots \}=\boldsymbol{v_{k-2}}+\boldsymbol{v_{k-1}}
\end{align}
なので、
[f(\boldsymbol{v_0})\quad f(\boldsymbol{v_1})\quad f(\boldsymbol{v_2}) \quad \cdots \quad f(\boldsymbol{v_{k-1}})]=
[\boldsymbol{v_0}\ \boldsymbol{v_1}\ \boldsymbol{v_2} \cdots \boldsymbol{v_{k-1}}]
\begin{bmatrix}
0 & 1 & 0 & 0 & \cdots & 0\\
0 & 0 & 1 & 0 & \cdots & 0\\
0 & 0 & 0 & 1 & \cdots & 0\\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & 0 & \cdots & 1\\
1 & 1 & 1 & 1 & \cdots & 1
\end{bmatrix}
と表せます。よって、$f$の基底 $\boldsymbol{v_0},\boldsymbol{v_1},\dots,\boldsymbol{v_{k-1}}$ に関する表現行列は、
\overbrace{
\begin{bmatrix}
0 & 1 & 0 & 0 & \cdots & 0\\
0 & 0 & 1 & 0 & \cdots & 0\\
0 & 0 & 0 & 1 & \cdots & 0\\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & 0 & \cdots & 1\\
1 & 1 & 1 & 1 & \cdots & 1
\end{bmatrix}
}^{k行k列}
です。この表現行列を$A$とします。
固有値と公比の関係
(トリボナッチ数列の一般項を線形代数とPythonで求めると同じですので、飛ばしていただいて構いません。)
\boldsymbol{p}=\{ r^{n} \}_{n=0}^{\infty}=\{1,r,r^2,\dots \}
が$V$に属したとすると、
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$を単位行列として、
\begin{align}
\varphi_k(\lambda)&=|A-\lambda I|\\
&=\begin{vmatrix}
-\lambda & 1 & 0 & 0 & \cdots & 0\\
0 & -\lambda & 1 & 0 & \cdots & 0\\
0 & 0 & -\lambda & 1 & \cdots & 0\\
\vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & 0 & \cdots & 1\\
1 & 1 & 1 & 1 & \cdots & 1-\lambda
\end{vmatrix}
\end{align}
1列目に沿って余因子展開すると、
\begin{align}
\varphi_k&=(-1)^{1+1}(-\lambda) \begin{vmatrix}
-\lambda & 1 & 0 & \cdots & 0\\
0 & -\lambda & 1 & \cdots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & 1\\
1 & 1 & 1 & \cdots & 1-\lambda
\end{vmatrix} + (-1)^{k+1}\cdot 1 \begin{vmatrix}
1 & 0 & 0 & \cdots & 0\\
-\lambda & 1 & 0 & \cdots & 0\\
0 & -\lambda & 1 & \cdots & 0\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
0 & 0 & 0 & \cdots & 1\\
\end{vmatrix} \\
&=-\lambda \varphi_{k-1} - (-1)^{k}\cdot 1\cdot 1\\
&=-\lambda \varphi_{k-1} - (-1)^{k}
\end{align}
と漸化式のように表せます。両辺を$(-1)^k$で割ると、
\begin{align}
\frac{\varphi_k}{(-1)^k} &= -\lambda \frac{\varphi_{k-1}}{(-1)^k} - 1\\
&=\lambda \frac{\varphi_{k-1}}{(-1)^{k-1}} - 1
\end{align}
$b_k = \displaystyle \frac{\varphi_k}{(-1)^k}$とおくと、
\begin{align}
b_k &= \lambda b_{k-1} - 1\\
b_k - \frac{1}{\lambda - 1} &= \lambda \left(b_{k-1} - \frac{1}{\lambda - 1}\right)\\
&=\cdots \\
&= \lambda^{k-2}\left(b_{2} - \frac{1}{\lambda - 1}\right)
\end{align}
であり、初項は、
\begin{align}
b_2- \frac{1}{\lambda - 1} &= \displaystyle \frac{\varphi_2}{(-1)^2} - \frac{1}{\lambda - 1}\\
&=\varphi_2 - \frac{1}{\lambda - 1}\\
&=\begin{vmatrix}
-\lambda & 1 \\
1 & 1-\lambda
\end{vmatrix}- \frac{1}{\lambda - 1}\\
&=\lambda^2 - \lambda - 1- \frac{1}{\lambda - 1}\\
&=\frac{\lambda^3-2\lambda^2}{\lambda-1}
\end{align}
です。よって、
\begin{align}
b_k &- \frac{1}{\lambda - 1}= \lambda^{k-2} \frac{\lambda^3-2\lambda^2}{\lambda-1}\\
b_k &= \frac{\lambda^{k+1}-2\lambda^k}{\lambda-1} + \frac{1}{\lambda - 1} \\
&= \frac{\lambda^{k+1}-2\lambda^k+1}{\lambda-1}\\
&=\lambda^k - \lambda^{k-1} - \lambda^{k-2} - \dots - \lambda^2 - \lambda - 1\quad (\because 組み立て除法)
\end{align}
が導出されました。方程式 $\varphi_k(\lambda) = (-1)^k b_k = 0$ は、$b_k = 0$ と同値であることがわかります。
しかし、$b_k=0$は簡単には解くことはできません。あとの章で計算機を用いて数値的に求めることにします。
$k$ 個の解を $\lambda_0,\lambda_1,\dots,\lambda_{k-1}$とすると、それぞれの解に対応する固有ベクトルである等比数列の
\boldsymbol{u_0}=\{ \lambda_0^{n} \}_{n=0}^{\infty}\\
\boldsymbol{u_1}=\{ \lambda_1^{n} \}_{n=0}^{\infty}\\
\vdots \\
\boldsymbol{u_{k-1}}=\{ \lambda_{k-1}^{n} \}_{n=0}^{\infty}
は$V$に属します。これらは $\lambda_0,\lambda_1,\dots,\lambda_{k-1}$ が全て異なるとすると、$f$の相異なる固有値に対応する固有ベクトルなので一次独立です。また、$\dim V=k$ です。よって、$\boldsymbol{u_0},\boldsymbol{u_1},\dots,\boldsymbol{u_{k-1}}$ は$V$の基底となります。
係数の決定
以上より、ある $c_0,c_1,\dots,c_{k-1}$が存在して、
\boldsymbol{a_n}=c_0\boldsymbol{u_0}+c_1\boldsymbol{u_1}+\dots+c_{k-1}\boldsymbol{u_{k-1}}
と表せるので、
\begin{align}
\boldsymbol{a_n}&=c_0\boldsymbol{u_0}+c_1\boldsymbol{u_1}+\dots+c_{k-1}\boldsymbol{u_{k-1}}\\
\Leftrightarrow \{a_0,a_1,\dots,a_{k-1},\dots \}&=c_0\{ \lambda_0^{n} \}_{n=0}^{\infty}+c_1\{ \lambda_1^{n} \}_{n=0}^{\infty}+\dots+c_{k-1}\{ \lambda_{k-1}^{n} \}_{n=0}^{\infty}\\
\Leftrightarrow \{0,0,\dots,1,\dots \}&=c_0\{ 1,\lambda_0,\dots,\lambda_0^{k-1},\dots \}+c_1\{1,\lambda_1,\dots,\lambda_1^{k-1},\dots \}+\dots+c_{k-1}\{ 1,\lambda_{k-1},\dots,\lambda_{k-1}^{k-1},\dots\} \\
\Leftrightarrow \{0,0,\dots,1,\dots \}&=\{ c_0+c_1+\dots+c_{k-1},\ c_0\lambda_0+c_1\lambda_1+\dots+c_{k-1}\lambda_{k-1},\dots, c_0\lambda_0^{k-1}+c_1\lambda^{k-1}+\dots+c_{k-1}\lambda^{k-1}, \dots\}
\end{align}
これを行列で表現すると、
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1\\
\lambda_0 & \lambda_1 & \lambda_2 & \cdots & \lambda_{k-1} \\
\lambda_0^2 & \lambda_1^2 & \lambda_2^2 & \cdots & \lambda_{k-1}^2\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\lambda_0^{k-1} & \lambda_1^{k-1} & \lambda_2^{k-1} & \cdots & \lambda_{k-1}^{k-1}
\end{bmatrix}
\begin{bmatrix}
c_1 \\ c_2 \\ c_3 \\ \vdots \\ c_{k-1}
\end{bmatrix}
=
\begin{bmatrix}
0 \\ 0 \\ 0 \\ \vdots \\ 1
\end{bmatrix}
です。左の行列の行列式は、ヴァンデルモンド行列式 $V_k$ であり、
V_k = \begin{vmatrix}
1 & 1 & 1 & \cdots & 1\\
\lambda_0 & \lambda_1 & \lambda_2 & \cdots & \lambda_{k-1} \\
\lambda_0^2 & \lambda_1^2 & \lambda_2^2 & \cdots & \lambda_{k-1}^2\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
\lambda_0^{k-1} & \lambda_1^{k-1} & \lambda_2^{k-1} & \cdots & \lambda_{k-1}^{k-1}
\end{vmatrix}
=\prod_{0\leq i < j \leq k-1}\big(\lambda_j - \lambda_i\big)
と知られています。クラメルの公式を用いると、$0\leq m\leq k-1$ に対して
\begin{align}
c_m &= \frac{\begin{vmatrix}
1 & 1 & \cdots & 1 & 0 & 1 & \cdots & 1\\
\lambda_0 & \lambda_1 & \cdots & \lambda_{m-1} & 0 & \lambda_{m+1} & \cdots & \lambda_{k-1} \\
\lambda_0^2 & \lambda_1^2 & \cdots & \lambda_{m-1}^2 & 0 & \lambda_{m+1}^2 & \cdots & \lambda_{k-1}^2\\
\vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
\lambda_0^{k-1} & \lambda_1^{k-1} & \cdots & \lambda_{m-1}^{k-1} & 1 & \lambda_{m+1}^{k-1} & \cdots & \lambda_{k-1}^{k-1}
\end{vmatrix}}{V_k}\\
&=\frac{(-1)^{m}\begin{vmatrix}
0 & 1 & 1 & \cdots & 1 & 1 & \cdots & 1\\
0 & \lambda_0 & \lambda_1 & \cdots & \lambda_{m-1} & \lambda_{m+1} & \cdots & \lambda_{k-1} \\
0 & \lambda_0^2 & \lambda_1^2 & \cdots & \lambda_{m-1}^2 & \lambda_{m+1}^2 & \cdots & \lambda_{k-1}^2\\
0 & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
1 & \lambda_0^{k-1} & \lambda_1^{k-1} & \cdots & \lambda_{m-1}^{k-1} & \lambda_{m+1}^{k-1} & \cdots & \lambda_{k-1}^{k-1}
\end{vmatrix}}{V_k}\quad (\because 列をm回交換)\\
&=\frac{(-1)^{m+k-1}\begin{vmatrix}
1 & \lambda_0^{k-1} & \lambda_1^{k-1} & \cdots & \lambda_{m-1}^{k-1} & \lambda_{m+1}^{k-1} & \cdots & \lambda_{k-1}^{k-1}\\
0 & 1 & 1 & \cdots & 1 & 1 & \cdots & 1\\
0 & \lambda_0 & \lambda_1 & \cdots & \lambda_{m-1} & \lambda_{m+1} & \cdots & \lambda_{k-1} \\
0 & \lambda_0^2 & \lambda_1^2 & \cdots & \lambda_{m-1}^2 & \lambda_{m+1}^2 & \cdots & \lambda_{k-1}^2\\
0 & \vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
0 & \lambda_0^{k-2} & \lambda_1^{k-2} & \cdots & \lambda_{m-1}^{k-2} & \lambda_{m+1}^{k-2} & \cdots & \lambda_{k-1}^{k-2}
\end{vmatrix}}{V_k}\quad (\because 行をk-1回交換)\\
&=\frac{(-1)^{m+k-1}\begin{vmatrix}
1 & 1 & \cdots & 1 & 1 & \cdots & 1\\
\lambda_0 & \lambda_1 & \cdots & \lambda_{m-1} & \lambda_{m+1} & \cdots & \lambda_{k-1} \\
\lambda_0^2 & \lambda_1^2 & \cdots & \lambda_{m-1}^2 & \lambda_{m+1}^2 & \cdots & \lambda_{k-1}^2\\
\vdots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\
\lambda_0^{k-2} & \lambda_1^{k-2} & \cdots & \lambda_{m-1}^{k-2} & \lambda_{m+1}^{k-2} & \cdots & \lambda_{k-1}^{k-2}
\end{vmatrix}}{V_k}\\
&=\frac{(-1)^{m+k-1}\displaystyle \prod_{0\leq i < j \leq k-1かつi\neq mかつj \neq m}\big(\lambda_j - \lambda_i\big)}{\displaystyle \prod_{0\leq i < j \leq k-1}\big(\lambda_j - \lambda_i\big)}\quad (\because ヴァンデルモンド行列式)\\
&=\frac{(-1)^{m+k-1}}{\displaystyle \prod_{0\leq i < j \leq k-1かつ(i=mまたはj=m)}\big(\lambda_j - \lambda_i\big)}\\
&=\frac{(-1)^{m+k-1}}{\displaystyle \prod_{0\leq i < m}\big(\lambda_m - \lambda_i\big)\prod_{m < j \leq k-1}\big(\lambda_j - \lambda_m\big)}\\
&=\frac{(-1)^{m+k-1}}{\displaystyle \prod_{0\leq i < m}\big(\lambda_m - \lambda_i\big)\prod_{m < j \leq k-1}(-1)\big(\lambda_m - \lambda_j\big)}\\
&=\frac{(-1)^{m+k-1}}{\displaystyle \prod_{0\leq i < m}\big(\lambda_m - \lambda_i\big)(-1)^{k-m-1}\prod_{m < j \leq k-1}\big(\lambda_m - \lambda_j\big)}\\
&=\frac{(-1)^{(m+k-1)-(k-m-1)}}{\displaystyle \prod_{0\leq i < m}\big(\lambda_m - \lambda_i\big)\prod_{m < j \leq k-1}\big(\lambda_m - \lambda_j\big)}\\
&=\frac{(-1)^{2m}}{\displaystyle \prod_{0\leq i < m}\big(\lambda_m - \lambda_i\big)\prod_{m < j \leq k-1}\big(\lambda_m - \lambda_j\big)}\\
&=\frac{1}{\displaystyle \prod_{0\leq i <m,m<i\leq k-1}\big(\lambda_m - \lambda_i\big)}
\end{align}
が求められます。
以上より、
\begin{align}
\boldsymbol{a_n} &= c_0 \boldsymbol{u_0} + c_1 \boldsymbol{u_1} +\dots+ c_{k-1} \boldsymbol{u_{k-1}}\\
&=\left\{ \frac{\lambda_0^n}{\displaystyle \prod_{1\leq i \leq k-1}\big(\lambda_0 - \lambda_i\big)} + \frac{\lambda_1^n}{\displaystyle \prod_{0\leq i <1,1 <j\leq k-1}\big(\lambda_1 - \lambda_i\big)} +\dots+ \frac{\lambda_{k-1}^n}{\displaystyle \prod_{0\leq i \leq k-2}\big(\lambda_{k-1} - \lambda_i\big)}\right\}_{n=0}^{\infty}\\
&=\left\{\sum_{i=0}^{k-1}\frac{\lambda_i^n}{\displaystyle \prod_{0\leq j <i, i< j\leq k-1}\big(\lambda_i - \lambda_j\big)}\right\}_{n=0}^{\infty}
\end{align}
なので、$\lambda^k - \lambda^{k-1} - \lambda^{k-2} - \dots - \lambda^2 - \lambda - 1 = 0$の$k$個の解を用いて、
a_n = \sum_{i=0}^{k-1}\frac{\lambda_i^n}{\displaystyle \prod_{0\leq j <i, i< j\leq k-1}\big(\lambda_i - \lambda_j\big)}
が導けました。 上で省略した1が固有値ではないことや固有値が全て異なることを含めた説明はPDFファイルをご覧ください。Qiitaでは埋め込みができなかったため他の場所にアップロードしました。
https://magolors.hatenablog.com/entry/2020/09/03/170050
実装例と計算量比較
以上までで、一般項の形はわかりましたが、具体的な $\lambda_0,\lambda_1,\dots,\lambda_{k-1}$ の値はわかっていません。5次以上の方程式には代数的解法は必ずしも存在しないので数値解を求めます。今回もPythonのnumpyライブラリを利用して求めます。以下の問題を考えます。
自然数 $K,N$ が与えられます。$K$ -ボナッチ数列の第 $N$ 項を出力してください。
実装例1 : 定義通りに計算
teigi.py
k, n = map(int, input().split())
a = [0] * max(n+1,k)
a[k-1] = 1
for i in range(k,n+1):
for j in range(1,k+1):
a[i] += a[i-j]
print(a[n])
定義通りに計算していくと、各 $i (K \leq i \leq N)$ に対して $K$ 回の加算が必要なので、$O(K(N-K)) = O(KN)$ となります。
実装例2 : 一般項を用いる
ippankou.py
import numpy as np
import warnings
warnings.filterwarnings('ignore') # 複素数の計算時にWarningが発生するので消しておく
k, n = map(int, input().split())
equation = [1]
for i in range(k):
equation.append(-1) # equation = [1, -1, -1, ... , -1] = lambda^k - lambda^{k-1} - ... - lambda -1 = 0
lam = np.roots(equation) # lambda_i
a = 0
for i in range(k):
prod = 1
for j in range(k):
if j != i:
prod *= lam[i] - lam[j]
a += np.power(lam[i],n) / prod
print(int(round(a)))
一般項を用いて計算すると、各 $i (0 \leq i \leq K)$ に対して、各 $j (0 \leq j < i, i < j \leq K)$ の乗算と $N$ 乗計算で $O(\log N)$ が必要なので、$O(K(K+\log N)) = O(K^2+K\log N)$ となります。この方法では、前回の $k=3$ のときと同じように、浮動小数点数による数値の$10^{14}$ 分の1程度の誤差が発生していると思われます。
実装例3 : 以下の漸化式を用いる(k 項の和を保存して最初と最後だけ入れ替え)
\begin{align}
a_{n+k} &= \sum_{i=n}^{n+k-1}a_i\\
a_{n+k+1} &= \sum_{i=n+1}^{n+k}a_i \\
&= a_{n+k}+\sum_{i=n}^{n+k-1}a_i-a_n\\
&= a_{n+k}+a_{n+k}-a_n\\
&= 2a_{n+k} - a_n
\end{align}
zenkashiki.py
k, n = map(int, input().split())
a = [0] * max(n+1,k)
a[k-1] = 1
a[k] = 1
for i in range(0,n-k):
a[i+k+1] = 2 * a[i+k] - a[i]
print(a[n])
このように、一般項を求めなくても $k$ 項の和を保存しておく累積和のようなアルゴリズムを考えれば計算量 $O(N)$ で求めることができます。
実行例
2 5
5
3 50
3122171529233
4 25
1055026
10 50
1082392024832
50 30
0
100 150
1125899906842618
1000 1010
1024
計算量の比較
| 実装1 | 実装2 | 実装3 |
|---|---|---|
| 定義通り | 一般項 | 漸化式 |
| $O(NK)$ | $O(K^2+K\log N)$ | $O(N)$ |
実装2と実装3の計算量の比較により、$N$ と $K^2+K\log N$ を比較してどちらの方が高速か考えてから計算すると良いと思います。
実行時間の比較
measure_time.py
import time
import numpy as np
import warnings
warnings.filterwarnings('ignore') # 複素数の計算時にWarningが発生するので消しておく
k, n = map(int, input().split())
m = 100 # loop count
"""
# 実装1
start = time.time()
for loop in range(m):
a = [0] * max(n+1,k)
a[k-1] = 1
for i in range(k,n+1):
for j in range(1,k+1):
a[i] += a[i-j]
elapsed_time = time.time() - start
elapsed_time /= m
print ("実装1 : {0}".format(elapsed_time) + "[sec]")
"""
# 実装2
start = time.time()
for loop in range(m):
equation = [1]
for i in range(k):
equation.append(-1) # equation = [1, -1, -1, ... , -1] = lambda^k - lambda^{k-1} - ... - lambda -1 = 0
lam = np.roots(equation) # lambda_i
a = 0
for i in range(k):
prod = 1
for j in range(k):
if j != i:
prod *= lam[i] - lam[j]
a += np.power(lam[i],n) / prod
elapsed_time = time.time() - start
elapsed_time /= m
print ("実装2 : {0}".format(elapsed_time) + "[sec]")
# 実装3
start = time.time()
for loop in range(m):
a = [0] * max(n+1,k)
a[k-1] = 1
a[k] = 1
for i in range(0,n-k):
a[i+k+1] = 2 * a[i+k] - a[i]
elapsed_time = time.time() - start
elapsed_time /= m
print ("実装3 : {0}".format(elapsed_time) + "[sec]")
計測結果例
# N ~ K^2+Klog N
10 500
実装2 : 0.00041536092758178713[sec]
実装3 : 0.0003914594650268555[sec]
# N > K^2+Klog N
3 10000
実装2 : 0.0002478790283203125[sec]
実装3 : 0.012493629455566407[sec]
# N < K^2+Klog N
100 500
実装2 : 0.01716961145401001[sec]
実装3 : 0.0002404618263244629[sec]
計算量比較でも考えていたように、$N > K^2+K\log N$ のときは実装2、$N < K^2+K\log N$ のときは実装3が高速と確かめられます。
参考文献
こちらのpdfを参考にさせていただきました。ありがとうございました。