Day 047-Q3 — Berlekamp-Massey + Kitamasa法（線形漸化式N項目高速計算）

2026-05-31 赤色 Master / Phase 8+ ★★★★★★★★★ Berlekamp-Massey / Kitamasa / Cayley-Hamilton

問題

長さ $M$ の数列 $a_0, a_1, \ldots, a_{M-1}$ が与えられる。この数列は未知の線形漸化式を満たす：

$$a_n = c_1 a_{n-1} + c_2 a_{n-2} + \cdots + c_K a_{n-K} \quad (n \ge K)$$

まず Berlekamp-Massey 法で最小漸化式を求め、次に Kitamasa 法（多項式累乗 + Cayley-Hamilton）で $T$ 個のクエリ $n_i$ に対して $a_{n_i} \pmod{998244353}$ を求めよ。

制約

$2 \le M \le 1000$
$1 \le T \le 10^5$
$0 \le n_i \le 10^{18}$
$0 \le a_i < 998244353$
時間制限: 3秒

入出力例

入力例 1 (フィボナッチ)

6 3
0 1 1 2 3 5
10
20
50

出力例 1

55
6765
12586269025

概念図: Berlekamp-Massey → Kitamasa 2段階パイプライン

① Berlekamp-Massey 入力: 数列 a_0..a_{M-1} LFSR 最短次数を発見 出力: [c_1, ..., c_K] O(M²) ② Kitamasa 法 特性多項式 p(x) を構成 x^n mod p(x) を多項式累乗 出力: r_0,...,r_{K-1} O(K² log n) per query ③ 復元 $a_n = \sum r_i \cdot a_i$ Cayley-Hamilton の線形性 Cayley-Hamilton 定理の活用 特性多項式: $p(x) = x^K - c_1 x^{K-1} - \cdots - c_K$ $x^K \equiv c_1 x^{K-1} + \cdots + c_K \pmod{p(x)}$ なので $x^n \bmod p(x)$ が計算可能 $a_n = r_0 a_0 + r_1 a_1 + \cdots + r_{K-1} a_{K-1}$ （$x^n \bmod p(x) = r_0 + r_1 x + \cdots$）

ヒント（段階的開示）

ヒント1: 方向性
2段階のアルゴリズムを組み合わせます：
1. Berlekamp-Massey (BM) 法：与えられた数列から最短の線形漸化式を $O(M^2)$ で発見
2. Kitamasa 法 / 多項式累乗：発見した漸化式を使い、$a_n$ を $O(K^2 \log n)$ で計算
ヒント2: アプローチ
  • BM 法の出力: $c_1, \ldots, c_K$（最短漸化式の係数）
  • 特性多項式: $p(x) = x^K - c_1 x^{K-1} - \cdots - c_K$
  • Kitamasa: $x^n \pmod{p(x)}$ を多項式繰り返し二乗法で計算（各ステップで多項式 mod を行う）
  • 結果の係数 $r_i$ を使って $a_n = \sum_{i=0}^{K-1} r_i \cdot a_i$
ヒント3: BM法のコア実装 
def berlekamp_massey(s):
    C, B = [1], [1]
    L, x, b = 0, 1, 1
    for n in range(len(s)):
        d = sum(C[j] * s[n-j] for j in range(L+1)) % MOD
        if d == 0:
            x += 1
        elif 2 * L <= n:
            T = C[:]
            coef = d * pow(b, MOD-2, MOD) % MOD
            # C -= coef * x^x * B
            ...
            L, B, b, x = n + 1 - L, T, d, 1
        else:
            # C -= coef * x^x * B
            ...
            x += 1
    return [-c % MOD for c in C[1:]]

模範解答 (Python)

import sys
input = sys.stdin.readline

MOD = 998244353

def berlekamp_massey(s):
    C, B = [1], [1]
    L, x, b = 0, 1, 1
    for n in range(len(s)):
        d = 0
        for j in range(min(L+1, len(C))):
            d = (d + C[j] * s[n-j]) % MOD
        if d == 0:
            x += 1
        elif 2 * L <= n:
            T = C[:]
            coef = d * pow(b, MOD-2, MOD) % MOD
            if len(C) < x + len(B):
                C += [0] * (x + len(B) - len(C))
            for i, v in enumerate(B):
                C[i+x] = (C[i+x] - coef * v) % MOD
            L, B, b, x = n + 1 - L, T, d, 1
        else:
            coef = d * pow(b, MOD-2, MOD) % MOD
            if len(C) < x + len(B):
                C += [0] * (x + len(B) - len(C))
            for i, v in enumerate(B):
                C[i+x] = (C[i+x] - coef * v) % MOD
            x += 1
    return [-c % MOD for c in C[1:]]

def poly_mul(a, b, mod):
    res = [0] * (len(a) + len(b) - 1)
    for i, x in enumerate(a):
        for j, y in enumerate(b):
            res[i+j] = (res[i+j] + x * y) % mod
    return res

def poly_mod_p(f, g, mod):
    f = f[:]
    dg = len(g) - 1
    for i in range(len(f) - 1, dg - 1, -1):
        if f[i] == 0:
            continue
        c = f[i]
        for j in range(dg):
            f[i - dg + j] = (f[i - dg + j] - c * g[j]) % mod
    return f[:dg]

def poly_pow_mod(base_poly, exp, char_poly, mod):
    result = [0] * len(char_poly)
    result[0] = 1
    base = base_poly[:]
    while exp:
        if exp & 1:
            result = poly_mod_p(poly_mul(result, base, mod), char_poly, mod)
        base = poly_mod_p(poly_mul(base, base, mod), char_poly, mod)
        exp >>= 1
    return result

def kitamasa(n, rec, init):
    K = len(rec)
    if n < K:
        return init[n] % MOD
    char = [-rec[K-1-j] % MOD for j in range(K)] + [1]
    xpoly = [0] * (K + 1)
    xpoly[1] = 1
    xpoly = poly_mod_p(xpoly, char, MOD)
    rn = poly_pow_mod(xpoly, n, char, MOD)
    ans = 0
    for i in range(K):
        ans = (ans + rn[i] * init[i]) % MOD
    return ans

def main():
    M, T = map(int, input().split())
    a = list(map(int, input().split()))
    rec = berlekamp_massey(a)
    K = len(rec)
    init = a[:K]
    out = []
    for _ in range(T):
        n = int(input())
        out.append(kitamasa(n, rec, init))
    print('\n'.join(map(str, out)))

main()

Step-by-Step 解説

1Berlekamp-Massey 法
与えられた数列から最短線形漸化式を発見するアルゴリズム。LFSR（Linear Feedback Shift Register）の最小次数を求める問題と同値。計算量 $O(M^2)$。
2特性多項式の構成
漸化式 $a_n = c_1 a_{n-1} + \cdots + c_K a_{n-K}$ の特性多項式: $p(x) = x^K - c_1 x^{K-1} - \cdots - c_K$。Cayley-Hamilton 定理: $x^K \equiv c_1 x^{K-1} + \cdots + c_K \pmod{p(x)}$。
3$x^n \bmod p(x)$ の計算
多項式の繰り返し二乗法。$x^n \bmod p(x) = r_0 + r_1 x + \cdots + r_{K-1} x^{K-1}$ を求める。各ステップで多項式乗算 + 多項式 mod を行う。乗算 $O(K^2)$（愚直）、全体 $O(K^2 \log n)$。
4$a_n$ の復元
$a_n = \sum_{i=0}^{K-1} r_i \cdot a_i$（線形性から）。この式の導出は「$a_n$ が $x^n$ の評価値に対応し、$x^n \bmod p(x)$ が $a_n$ の線形結合を与える」という事実から。
5NTT による高速化
$K$ が大きい（$K \sim 500$）とき、多項式乗算を NTT で $O(K \log K)$ にすると全体 $O(K \log K \log n)$ で大幅高速化。

計算量

Berlekamp-Massey: $O(M^2)$
Kitamasa per query: $O(K^2 \log n)$（愚直多項式乗算）
Kitamasa per query: $O(K \log K \log n)$（NTT版）
全体: $O(M^2 + T \cdot K^2 \log n)$

よくあるミス

ミス原因正しい書き方
BM 法の係数の符号が逆$C$ の定義が $a_n + c_1 a_{n-1} + \ldots = 0$ 形式[-c % MOD for c in C[1:]] で正の $c_i$ に変換
特性多項式の表現ミス次数の低い方から並べるか高い方かchar[i] が $x^i$ の係数であることを一貫させる
poly_pow_mod の初期値$x^0 = 1$ を誤って初期化result[0] = 1
$n < K$ の境界条件を忘れるkitamasa が初期値より小さい $n$ を参照if n < K: return init[n]

次のステップ

  • 発展問題: $K = 10^3$ の場合に NTT ベースの多項式乗算で高速化
  • 関連: Day038 Q4（Cayley-Hamilton定理 + 線形漸化式高速計算）の復習
  • 応用: グラフの隣接行列の冪乗をキタマサで高速化

自己評価