問題
$d$ 次の線形漸化式 $a_n = c_1 a_{n-1} + c_2 a_{n-2} + \cdots + c_d a_{n-d}$ と初期値 $a_0, \ldots, a_{d-1}$、整数 $N$（$0 \le N \le 10^{18}$）が与えられる。$a_N \bmod (10^9 + 7)$ を求めよ。
制約
| パラメータ | 範囲 |
|---|---|
| $d$ | $1 \le d \le 500$ |
| $N$ | $0 \le N \le 10^{18}$ |
| 係数 | $0 \le c_i, a_i < 10^9 + 7$ |
| 制限 | 時間 3 sec / メモリ 256MB |
入出力例
入力例 1 (Fibonacci)
2 10
1 1
0 1出力例 1
55入力例 2 (d=3)
3 1000000000000000000
1 1 1
0 1 1出力例 2
（mod 10^9+7 の値）概念図: Cayley-Hamilton 定理の応用
コンパニオン行列 $M$ の特性多項式 $p(\lambda)$ について $p(M) = 0$（零行列）。よって $M^N$ は $p$ で割った余りの多項式で表せる。
ヒント（段階的開示）
ヒント1: 方向性
行列累乗では $O(d^3 \log N)$ となり $d=500$ では TLE。Cayley-Hamilton 定理: 特性多項式 $p(\lambda)$ でコンパニオン行列 $M$ を割ると $p(M)=0$。よって $M^N \bmod p$ を多項式として計算すると $O(d^2 \log N)$ に改善できる（Kitamasa 法）。
ヒント2: アプローチ（Kitamasa 法）
1. 特性多項式 $p(\lambda) = \lambda^d - c_1\lambda^{d-1} - \cdots - c_d$ を構築（昇順係数）
2. $x^N \bmod p(x)$ を多項式の繰り返し二乗法で計算
3. 結果 $r(x) = \sum_{k=0}^{d-1} r_k x^k$ として $a_N = \sum_k r_k \cdot a_k$
2. $x^N \bmod p(x)$ を多項式の繰り返し二乗法で計算
3. 結果 $r(x) = \sum_{k=0}^{d-1} r_k x^k$ として $a_N = \sum_k r_k \cdot a_k$
ヒント3: poly_mod の実装
def poly_mod(f, p, d):
"""f mod p（p はモニック、次数 d）、昇順係数"""
f = f[:]
while len(f) > d:
coef = f[-1]
if coef == 0:
f.pop(); continue
deg = len(f) - 1
# x^deg = coef * (c₁x^{d-1} + ... + cₐ) に変換
for i in range(d):
f[deg - d + i] = (f[deg - d + i] - coef * p[i]) % MOD
f.pop()
while f and f[-1] == 0: f.pop()
return f or [0]模範解答 (Python)
import sys
input = sys.stdin.readline
MOD = 10**9 + 7
def poly_mod(f, g, dg):
"""f mod g（g はモニック、次数 dg）、昇順係数"""
f = f[:]
while len(f) > dg:
coef = f[-1]
if coef == 0:
f.pop()
continue
deg = len(f) - 1
for i in range(dg + 1):
f[deg - dg + i] = (f[deg - dg + i] - coef * g[i]) % MOD
f.pop()
while f and f[-1] == 0:
f.pop()
return f if f else [0]
def poly_mul(f, g):
if not f or not g:
return [0]
result = [0] * (len(f) + len(g) - 1)
for i, fi in enumerate(f):
if fi == 0:
continue
for j, gj in enumerate(g):
result[i+j] = (result[i+j] + fi * gj) % MOD
return result
def poly_mul_mod(f, g, p, dg):
return poly_mod(poly_mul(f, g), p, dg)
def poly_pow_mod(n, p, dg):
"""x^n mod p(x)、p はモニックで次数 dg"""
result = [1] # 定数多項式 1
base = [0, 1] # x
while n:
if n & 1:
result = poly_mul_mod(result, base, p, dg)
base = poly_mul_mod(base, base, p, dg)
n >>= 1
return result
def solve():
line1 = input().split()
d, N = int(line1[0]), int(line1[1])
c = list(map(int, input().split())) # c[0]=c1, ..., c[d-1]=cd
a = list(map(int, input().split())) # a[0]=a0, ..., a[d-1]=a(d-1)
if N < d:
print(a[N] % MOD)
return
# 特性多項式（昇順係数）:
# p(x) = x^d - c1*x^(d-1) - ... - cd
# p[i] = x^i の係数
# p[d] = 1, p[i] = -c[d-1-i] for i < d
char_poly = [0] * (d + 1)
char_poly[d] = 1
for i in range(d):
char_poly[i] = (-c[d-1-i]) % MOD
r = poly_pow_mod(N, char_poly, d)
ans = 0
for k in range(min(len(r), d)):
ans = (ans + r[k] * a[k]) % MOD
print(ans)
solve()
Step-by-Step 解説
1特性多項式の構築
漸化式 $a_n = \sum_{k=1}^{d} c_k a_{n-k}$ の特性多項式 $p(\lambda) = \lambda^d - c_1\lambda^{d-1} - \cdots - c_d$。昇順係数で
漸化式 $a_n = \sum_{k=1}^{d} c_k a_{n-k}$ の特性多項式 $p(\lambda) = \lambda^d - c_1\lambda^{d-1} - \cdots - c_d$。昇順係数で
char_poly[i] = $x^i$ の係数。
2poly_mod の実装
$f$ の最高次項から $p$ の倍数を引いて次数を下げる。$p$ はモニック（最高次係数 1）なので除算不要。
$f$ の最高次項から $p$ の倍数を引いて次数を下げる。$p$ はモニック（最高次係数 1）なので除算不要。
3$x^N \bmod p(x)$ の計算
繰り返し二乗法で多項式累乗。各ステップで `poly_mul_mod` を適用。$O(d^2 \log N)$。
繰り返し二乗法で多項式累乗。各ステップで `poly_mul_mod` を適用。$O(d^2 \log N)$。
4初期値との内積
結果 $r = \sum_{k=0}^{d-1} r_k x^k$ として $a_N = \sum_{k=0}^{d-1} r_k \cdot a_k$。これが Cayley-Hamilton 定理の帰結。$O(d)$。
結果 $r = \sum_{k=0}^{d-1} r_k x^k$ として $a_N = \sum_{k=0}^{d-1} r_k \cdot a_k$。これが Cayley-Hamilton 定理の帰結。$O(d)$。
5特殊ケース処理
$N < d$ のとき初期値をそのまま返す。
$N < d$ のとき初期値をそのまま返す。
計算量
poly_mul: $O(d^2)$
poly_mod: $O(d^2)$
poly_pow_mod: $O(d^2 \log N)$
合計: $O(d^2 \log N)$ — $d=500, N=10^{18}$ で約 $10^7 \times 60 = 6\times10^8$（定数小）
poly_mod: $O(d^2)$
poly_pow_mod: $O(d^2 \log N)$
合計: $O(d^2 \log N)$ — $d=500, N=10^{18}$ で約 $10^7 \times 60 = 6\times10^8$（定数小）
よくあるミス
| ミス | 原因 | 正しい書き方 |
|---|---|---|
| 昇順/降順の混同 | char_poly の添字ミス | char_poly[i] = $x^i$ の係数（昇順） |
| poly_mod で最高次係数を割る | 非モニックな場合の処理 | 特性多項式はモニック（最高次=1）なので除算不要 |
| $N < d$ の特殊ケース未処理 | 配列外参照 | if N < d: print(a[N]); return |
| 多項式乗算の結果サイズ | 2d-1 まで増える中間値 | result = [0] * (len(f) + len(g) - 1) |
次のステップ
- 発展問題: Berlekamp-Massey 法で漸化式の係数 $c_i$ 自体を数列から推定 + Kitamasa 組み合わせ
- 応用: グラフ上の経路数え上げ（隣接行列累乗の代替）
- 理論: Minimal Polynomial の理解（最小次数の消去多項式）