Day 016-Q1 | Berlekamp-Massey法

問題

数列 $a_0, a_1, \ldots, a_{N-1}$（$N \le 2000$）が与えられる。この数列は、ある長さ $L$ の線形漸化式

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

を満たすという（$L$ は最小のものとする）。

Berlekamp-Massey法を実装し、最小の線形漸化式の係数 $[c_1, c_2, \ldots, c_L]$ を求めよ。

また、求めた漸化式を用いて $a_M$（$M \le 10^{18}$）を $998244353$ で割った余りを出力せよ。

入力形式

N M
a_0 a_1 ... a_{N-1}

制約

$1 \le N \le 2000$

$0 \le M \le 10^{18}$

$0 \le a_i < 998244353$

長さ $N$ の数列が少なくとも1つの線形漸化式を持つことが保証される

入出力例

入力例 1

8 100
1 1 2 3 5 8 13 21

出力例 1

927197596

（$F_{100} \mod 998244353$）

入力例 2

6 10
1 2 4 8 16 32

出力例 2

ヒント (段階的開示)

ヒント1: 方向性

Berlekamp-Massey は、数列を説明する最短の LFSR（線形フィードバックシフトレジスタ）を $O(N^2)$ で求めるアルゴリズム。

ヒント2: アプローチ

BM 法で係数列 $c$ を求めた後、Cayley-Hamilton定理 + Kitamasa法または行列累乗で $a_M$ を $O(L^2 \log M)$ で求める。

ヒント3: 誘導

# Berlekamp-Massey の核心
# C: 現在の最短 LFSR 係数（1-indexed、C[0]=1 は形式的）
# b: 最後にCを更新したときのLFSR
# L: 現在のLFSR長
# x: 最後の更新からのステップ数

模範解答 (Python)

import sys
from typing import List

MOD = 998244353

def power(a: int, b: int, mod: int) -> int:
    res = 1
    a %= mod
    while b:
        if b & 1:
            res = res * a % mod
        a = a * a % mod
        b >>= 1
    return res

def berlekamp_massey(s: List[int]) -> List[int]:
    """最短LFSRの係数を返す（1-indexed, c[i]がa[n-i]にかかる係数）"""
    n = len(s)
    C = [1]  # C[0]=1, C[i] = -c_i の形式
    B = [1]
    L = 0
    x = 1
    b = 1

    for i in range(n):
        # 差分 d = a[i] - Σ c_j * a[i-j]
        d = s[i]
        for j in range(1, L + 1):
            d = (d + C[j] * s[i - j]) % MOD
        if d == 0:
            x += 1
            continue
        T = C[:]
        coef = d * power(b, MOD - 2, MOD) % MOD
        # C = C - coef * x^shift * B
        if len(C) < len(B) + x:
            C += [0] * (len(B) + x - len(C))
        for j in range(len(B)):
            C[j + x] = (C[j + x] - coef * B[j]) % MOD
        if 2 * L <= i:
            L = i + 1 - L
            B = T
            b = d
            x = 1
        else:
            x += 1

    # C = [1, -c1, -c2, ...] -> 返すのは [c1, c2, ..., cL]
    return [(-C[i]) % MOD for i in range(1, L + 1)]


def nth_term(coefs: List[int], init: List[int], n: int) -> int:
    """漸化式 a[i] = Σ coefs[j]*a[i-1-j] で a[n] を O(L^2 log n) で求める"""
    L = len(coefs)
    if n < len(init):
        return init[n] % MOD

    # 初期数列が L 未満なら拡張
    a = list(init[:L])
    while len(a) < L:
        nxt = 0
        for j in range(len(a)):
            nxt = (nxt + coefs[j] * a[-(j+1)]) % MOD
        a.append(nxt)

    # Kitamasa法: 多項式の mod 演算で x^n を求める
    def poly_mul(p, q, mod_poly):
        """p * q mod mod_poly"""
        res = [0] * (len(p) + len(q) - 1)
        for i, pi in enumerate(p):
            for j, qj in enumerate(q):
                res[i + j] = (res[i + j] + pi * qj) % MOD
        # mod mod_poly: x^L = Σ coefs[j] * x^(L-1-j)
        while len(res) > len(mod_poly):
            lead = res.pop()
            if lead == 0:
                continue
            deg = len(res)  # res の次数
            for k in range(len(mod_poly)):
                res[deg - len(mod_poly) + k] = (
                    res[deg - len(mod_poly) + k] + lead * mod_poly[k]
                ) % MOD
        while len(res) < L:
            res.append(0)
        return res[:L]

    # mod_poly は漸化式に対応: x^L - c1*x^(L-1) - ... - cL
    # poly_mul で mod をとるとき係数として使う
    mod_poly = coefs[:]  # coefs[0]=c1, ..., coefs[L-1]=cL

    # x^n を求める (binary lifting)
    result = [0] * L
    result[0] = 1  # = 1 (= x^0 を表す単項式)
    base = [0] * L
    if L > 1:
        base[1] = 1  # = x
    else:
        base[0] = coefs[0]  # x = c1 (mod x - c1)

    m = n
    while m:
        if m & 1:
            result = poly_mul(result, base, mod_poly)
        base = poly_mul(base, base, mod_poly)
        m >>= 1

    # result は x^n の係数列: x^n ≡ Σ result[i] * x^i
    # a[n] = Σ result[i] * a[i]
    ans = 0
    for i in range(L):
        ans = (ans + result[i] * a[i]) % MOD
    return ans


def main():
    input_data = sys.stdin.read().split()
    N, M = int(input_data[0]), int(input_data[1])
    a = [int(x) for x in input_data[2:2+N]]

    coefs = berlekamp_massey(a)
    print(nth_term(coefs, a, M))

if __name__ == "__main__":
    main()

Step-by-Step 解説

1LFSR（線形フィードバックシフトレジスタ）の理解
線形漸化式 $a_i = c_1 a_{i-1} + \cdots + c_L a_{i-L}$ は LFSR そのもの。BM 法は観測した数列から最短の LFSR を同定する。

2Berlekamp-Massey の反復

$C$ を「現在の候補多項式」として管理
各ステップで「差分 $d$」を計算。$d \neq 0$ なら $C$ を更新
更新量は $d \cdot b^{-1}$ にシフトした旧多項式 $B$

3Kitamasa法で $a_M$ を高速計算

$a_M = $ 係数ベクトル $r$ と初期値 $[a_0, \ldots, a_{L-1}]$ の内積
$r$ は多項式冪 $x^M \bmod (x^L - c_1 x^{L-1} - \cdots - c_L)$ で得られる
多項式冪を繰り返し二乗法で $O(L^2 \log M)$

4計算量

BM 法: $O(N^2)$
Kitamasa: $O(L^2 \log M)$（$L \le N/2$）
全体: $O(N^2 + L^2 \log M)$

よくあるミス

ミス	原因	正しい書き方
`b` の逆元計算忘れ	$d/b$ のつもりが $d \cdot b$ になる	`power(b, MOD-2, MOD)`
mod_poly の次数ずれ	`x^L` の削減で添字オフセットミス	`res[deg - L + k]` で正確に
M < L のケース	BM 後の初期値 `a[0..L-1]` が足りない	`init[n]` を直接返す

次のステップ

発展問題: Berlekamp-Massey + FPS の組み合わせで「$N$ 項の和を $O(L \log L \log N)$」で求める

自己評価

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

Day 016-Q1 — Berlekamp-Massey法（線形漸化式の係数同定）

問題