【例題】数列$t_n$が以下の漸化式で表されるとき$t_{10^{12}}$ mod $10^4$を求めよ
\begin{align}
&t_0=0,\space t_1=t_2=1,\space t_3=4 \\
&\large t_n = 2t_{n-1}+2t_{n-2}-2t_{n-3}+t_{n-4}
\end{align}
この問題は$n$の値が$10^{12}$と非常に大きいのとmod計算がキーとなります。
再帰プログラム
まずは安易に再帰プログラム。lru_cacheを付けても$10^3$で撃沈。
from functools import lru_cache
@lru_cache(maxsize=None)
def t(n):
if n < 4: return [0,1,1,4][n]
return (2*t(n-1)+2*t(n-2)-2*t(n-3)+t(n-4)) % 10**4
print(t(10**2))
# 8344
ループ実装
次にサイクルバッファでループ実装。 なんとか$10^8$までが精いっぱい。
T = [0,1,1,4]
for n in range(4,N+1):
T[n%4] = (2*T[(n-1)%4]+2*T[(n-2)%4]-2*T[(n-3)%4]+T[(n-4)%4]) % 10**4
print(T[n%4])
# 688
行列で実装
漸化式は以下の行列の計算に変換できます。
\begin{pmatrix}
t_{n} \\t_{n-1} \\ t_{n-2} \\ t_{n-3}
\end{pmatrix} =
\begin{pmatrix}
2 & 2 & -2 & 1\\
1 & 0 & 0 & 0\\
0 & 1 & 0 & 0\\
0 & 0 & 1 & 0
\end{pmatrix}
\begin{pmatrix}
t_{n-1} \\t_{n-2} \\ t_{n-3} \\ t_{n-4}
\end{pmatrix}
プログラムで書きやすいように$t(n)$を横長の行列に変換します。
\begin{align}
&\begin{pmatrix}
t_{n} & t_{n-1} & t_{n-2} & t_{n-3}
\end{pmatrix} = \\
&\begin{pmatrix}
t_{n-1} & t_{n-2} & t_{n-3} & t_{n-4}
\end{pmatrix} \times
\begin{pmatrix}
2 & 1 & 0 & 0\\
2 & 0 & 1 & 0\\
-2 & 0 & 0 & 1\\
1 & 0 & 0 & 0
\end{pmatrix}
\end{align}
行列のべき乗
元の問題を解くためには以下の行列のべき乗を計算すれば良いことになります。
\begin{pmatrix}
2 & 1 & 0 & 0\\
2 & 0 & 1 & 0\\
-2 & 0 & 0 & 1\\
1 & 0 & 0 & 0
\end{pmatrix} ^{10^{12}}\space mod \space 10^4\\
これは「繰り返し二乗法で行列のべき乗の高速化」で作ったpowmtxを使うと高速に計算できるので、これを使ってプログラムにすると以下のようになり。1秒以内で答えがでました。
M= 10**4
N = 10**12
mtx = np.matrix([[2,2,-2,1],[1,0,0,0],[0,1,0,0],[0,0,1,0]], dtype='object').T
tv = np.matrix([4,1,1,0])
tv = tv*powmtx(mtx,N, M) % M
print(f"tn={tv[0,3]}")
# 6928
（開発環境：Google Colab）
この考え方はEuler Project: Problem 237を解くのに役に立ちます。