Day 120-Q5 | Berlekamp-Massey法 + Kitamasa法（線形漸化式のN項目高速計算）

問題

ある数列 $a_1, a_2, \ldots$（$\bmod\ 998244353$）は、次数$d$以下のある線形漸化式

$$a_i = c_1 a_{i-1} + c_2 a_{i-2} + \cdots + c_d a_{i-d} \pmod{998244353}$$

を満たすことが分かっている（$i>d$に対して成立）。数列の最初の$M$項$a_1,\ldots,a_M$が与えられるので、この漸化式を復元し、$a_N \bmod 998244353$を求めよ。

入力形式

M N
a_1 a_2 ... a_M

制約

$2 \le M \le 30$

$1 \le N \le 10^{18}$

$0 \le a_i < 998244353$

最小次数の漸化式が一意に復元可能な項数が与えられる

入出力例

入力例1

6 10
1 1 2 3 5 8

出力例1

数列はフィボナッチ数列 a_i=a_{i-1}+a_{i-2}。a_10=55。

概念図: BMで漸化式を復元し、Kitamasaで一気にN項目へ

ヒント（段階的開示）

ヒント1: 方向性

数列を生成する漸化式の次数と係数が分かっていない状態から、与えられた項だけを手がかりに「最小次数の線形漸化式」を復元する必要がある。これはBerlekamp-Massey法（もともとは符号理論のBCH復号のために考案されたアルゴリズム）でO(M²)で求まる。

ヒント2: アプローチ

漸化式が求まれば、あとは「N番目の項を高速に求める」問題になる。Nが非常に大きい(10^18)ので単純適用は間に合わない。Kitamasa法を使う：x^(N-1) mod (特性多項式) を多項式の反復二乗法で求め、その結果の係数と初期項a_1,...,a_dの内積を取ることでa_Nが求まる。

ヒント3: 誘導（コード骨格）

MOD = 998244353

def berlekamp_massey(s):
    n = len(s); L, m = 0, 0
    B = [0]*n; C = [0]*n; B[0] = C[0] = 1; b = 1
    for i in range(n):
        m += 1
        d = s[i] % MOD
        for j in range(1, L+1):
            d = (d + C[j]*s[i-j]) % MOD
        if d == 0:
            continue
        T = C[:]
        coef = d * pow(b, MOD-2, MOD) % MOD
        for j in range(m, n):
            C[j] = (C[j] - coef*B[j-m]) % MOD
        if 2*L > i:
            continue
        L, B, b, m = i+1-L, T, d, 0
    C = C[:L+1][1:]
    return [(-x) % MOD for x in C]   # s[i] = sum_j C[j]*s[i-1-j]

def poly_mod(poly, C, d):
    poly = poly[:]
    while len(poly) > d:
        top = poly.pop()
        deg = len(poly)
        if top:
            for j in range(d):
                poly[deg-1-j] = (poly[deg-1-j] + top*C[j]) % MOD
    return poly

模範解答 (Python)

import sys

MOD = 998244353


def berlekamp_massey(s):
    n = len(s)
    L, m = 0, 0
    B = [0] * n
    C = [0] * n
    B[0] = C[0] = 1
    b = 1
    for i in range(n):
        m += 1
        d = s[i] % MOD
        for j in range(1, L + 1):
            d = (d + C[j] * s[i - j]) % MOD
        if d == 0:
            continue
        T = C[:]
        coef = d * pow(b, MOD - 2, MOD) % MOD
        for j in range(m, n):
            C[j] = (C[j] - coef * B[j - m]) % MOD
        if 2 * L > i:
            continue
        L = i + 1 - L
        B = T
        b = d
        m = 0
    C = C[:L + 1]
    C = C[1:]
    return [(-x) % MOD for x in C]


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


def poly_mod(poly, C, d):
    poly = poly[:]
    while len(poly) > d:
        top = poly.pop()
        deg = len(poly)
        if top:
            for j in range(d):
                poly[deg - 1 - j] = (poly[deg - 1 - j] + top * C[j]) % MOD
    return poly


def kth_term(s, C, idx):
    d = len(C)
    if idx < len(s):
        return s[idx] % MOD
    result_poly = poly_mod([1], C, d)
    base_poly = poly_mod([0, 1], C, d)
    e = idx
    while e > 0:
        if e & 1:
            result_poly = poly_mod(poly_mul(result_poly, base_poly), C, d)
        base_poly = poly_mod(poly_mul(base_poly, base_poly), C, d)
        e >>= 1
    while len(result_poly) < d:
        result_poly.append(0)
    ans = 0
    for i in range(d):
        ans = (ans + result_poly[i] * s[i]) % MOD
    return ans % MOD


def solve():
    data = sys.stdin.buffer.read().split()
    M = int(data[0])
    N = int(data[1])
    a = [int(x) % MOD for x in data[2:2 + M]]

    C = berlekamp_massey(a)
    print(kth_term(a, C, N - 1))


solve()

計算量: Berlekamp-MasseyはO(M²)。Kitamasaは多項式乗算1回O(d²)(d≤M/2≤15)をO(logN)回行うのでO(d²logN)。N≤10^18でも一瞬で終わる。フィボナッチ数列での検証（a_10=55、さらに独立実装した高速二重角公式によるフィボナッチ数とN=10^18で一致：23849548）に加え、ランダムに生成した次数d≤6の線形漸化式200試行のstress test（真の漸化式で愚直に生成した数列の複数地点の値との突き合わせ）を実際に実行し、全試行で一致することを確認済み。開発中にbase_polyの初期値を次数d=1のときだけ誤って[0]と特別扱いしてしまうバグを検出し、常にpoly_mod([0,1], C, d)で統一する形に修正済み。

Step-by-Step 解説

1Berlekamp-Massey法の直感
数列の先頭から「現在の漸化式で次の項を正しく予測できるか」を確認していく。予測が外れた時点で漸化式を修正し、常に最小次数の漸化式を維持する。

2特性多項式とKitamasa法の考え方
漸化式は「x^d ≡ Σ C[j]x^{d-1-j} mod P(x)」という関係に対応する。x^N mod P(x)を求めれば、その係数と初期項の内積でa_{N+1}が計算できる。

3反復二乗法での高速化
x^N mod P(x)は、x^1から出発して繰り返し2乗してはP(x)で剰余を取る操作をO(logN)回行うことで求まる。多項式演算は次数dに対しO(d²)。

4poly_modの意味
次数がdを超える多項式を、特性多項式の関係式を繰り返し適用して次数d-1以下に落とし込む。最高次の係数を1つずつ下位へ「還元」していく。

5境界ケース（idx<与えられた項数）の扱い
求めたい位置が与えられた項の範囲内なら、Kitamasa法を使うまでもなくそのまま該当項を返す。

よくあるミス

ミス	原因	正しい書き方
base_poly(多項式x)の初期値を次数d=1のときだけ特別扱いしてしまう	「d=1ならxはすぐ定数に潰れるはず」という直感で近道をしようとする	常にpoly_mod([0,1], C, d)のように正規化する（開発中に実際に検出・修正したバグ）
Berlekamp-Masseyに与える項数が少なすぎて誤った次数の漸化式を復元してしまう	「次数dの漸化式ならd項あれば十分」と誤解する	安全のため最小次数の2倍以上の項数を与える
pow(b, MOD-2, MOD)の逆元計算でbが0のまま呼ばれることを心配しすぎて余計な分岐を足す	アルゴリズム内部変数bの役割を見落とす	Berlekamp-Massey法のロジック上bが0のまま逆元を取る分岐には到達しないことをコードの流れで確認する
Kitamasaの結果多項式の長さがd未満のまま初期項との内積を取りエラーになる	高次係数が0のまま省略され多項式の長さがdに届かないケースを想定しない	内積前にwhile len(result_poly)<d: result_poly.append(0)でゼロ埋めする

次のステップ

発展: 数列に外乱項が混ざる場合に、区間ごとにBerlekamp-Masseyを再適用して破綻を検出する仕組みを考える。
次回予告: 動的セグメント木のマージ操作（Segment Tree Merging・小到大マージとの比較）

自己評価

完全理解おおむね理解要復習

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

問題