Day 090-Q5 — Berlekamp-Welch Algorithm(Reed-Solomon誤り訂正)

2026-07-13 赤色 Master / Phase 8+ ★★★★★★★★★ 有限体・線形代数・符号理論

問題

$\mathrm{GF}(p)$ 上で、次数 $K$ 未満の多項式 $f(x)=c_0+c_1x+\dots+c_{K-1}x^{K-1}$ を点 $x=0,\dots,N-1$ で評価した符号語のうち、最大 $E=\lfloor(N-K)/2\rfloor$ 個が改ざんされた値 $y_0,\dots,y_{N-1}$ が届いた。元の係数 $c_0,\dots,c_{K-1}$ を復元せよ。

制約

パラメータ範囲備考
$p$素数, $3\le p\le10^6$
$N,K$$1\le K\le N\le200$, $N-K\ge2$符号語長・メッセージ長
改ざん数$\le\lfloor(N-K)/2\rfloor$保証される

入出力例

入力例1

17 7 3
3 2 0 10 11 2 1

出力例1

3 5 1

入力例2

17 5 2
5 8 11 14 0

出力例2

5 3

例1: 正しい符号語$[3,9,0,10,5,2,1]$の位置1,4が改ざん。$f(x)=1x^2+5x+3\pmod{17}$を復元。

概念図: 誤り位置検出多項式 $E(x)$ で誤りを無効化

$y_i \cdot E(x_i) = Q(x_i)$ がすべての点で成立するよう線形連立方程式を解く x=0 y=3 x=1 y=2 (誤り) x=2 y=0 x=3 y=10 x=4 y=11 (誤り) x=5 y=2 $E(x)=(x-1)(x-4)$ なら誤り位置で $E(x_i)=0$ → 正常点・誤り点を区別せず統一的に $Q(x)=f(x)E(x)$ を解ける 最後に $f(x)=Q(x)/E(x)$(余り0)で元の係数を復元

ヒント

ヒント1(方向性)

誤り位置を全探索すると $\binom{N}{E}$ 通りになり非効率。誤り位置を明示的に特定せず、連立方程式を解くだけで元の多項式を求めるのがBerlekamp-Welchの発想。

ヒント2(アプローチ)

誤り位置を根に持つモニック多項式 $E(x)$(次数$e$)を導入すると、$y_i\cdot E(x_i)=Q(x_i)$($Q=fE$、次数$e+K-1$)が全ての点で成立する。これは$E,Q$の係数について線形なので、連立1次方程式として解ける。

ヒント3(ほぼ答え)
# 未知数: E(x)=x^e+sum_{j

模範解答

import sys

def modinv(a, p):
    return pow(a, p - 2, p)

def gauss_solve(A, b, p):
    n = len(A)
    M = [row[:] + [b[i]] for i, row in enumerate(A)]
    for col in range(n):
        piv = None
        for r in range(col, n):
            if M[r][col] % p != 0:
                piv = r; break
        M[col], M[piv] = M[piv], M[col]
        inv = modinv(M[col][col], p)
        M[col] = [(x * inv) % p for x in M[col]]
        for r in range(n):
            if r != col and M[r][col] % p != 0:
                factor = M[r][col]
                M[r] = [(M[r][k] - factor * M[col][k]) % p for k in range(n + 1)]
    return [M[i][n] % p for i in range(n)]

def poly_divmod(num, den, p):
    num = num[:]
    while len(den) > 0 and den[-1] == 0:
        den.pop()
    quotient = [0] * (max(0, len(num) - len(den)) + 1)
    while len(num) - 1 >= len(den) - 1 and any(c != 0 for c in num):
        deg_diff = (len(num) - 1) - (len(den) - 1)
        coef = (num[-1] * modinv(den[-1], p)) % p
        quotient[deg_diff] = coef
        for i, c in enumerate(den):
            num[deg_diff + i] = (num[deg_diff + i] - coef * c) % p
        while len(num) > 1 and num[-1] == 0:
            num.pop()
    return quotient, num

def berlekamp_welch(xs, ys, k, e, p):
    n = len(xs)
    num_unknown = e + (e + k)
    A, b = [], []
    for i in range(n):
        x, y = xs[i], ys[i]
        row = [0] * num_unknown
        for j in range(e):
            row[j] = (y * pow(x, j, p)) % p
        for j in range(e + k):
            row[e + j] = (-pow(x, j, p)) % p
        rhs = (-y * pow(x, e, p)) % p
        A.append(row); b.append(rhs)
    sol = gauss_solve(A, b, p)
    e_coeffs = sol[:e] + [1]
    q_coeffs = sol[e:]
    return e_coeffs, q_coeffs

def main():
    data = sys.stdin.read().split()
    idx = 0
    p = int(data[idx]); idx += 1
    n = int(data[idx]); idx += 1
    k = int(data[idx]); idx += 1
    ys = [int(data[idx + i]) for i in range(n)]
    xs = list(range(n))

    e = (n - k) // 2
    e_coeffs, q_coeffs = berlekamp_welch(xs, ys, k, e, p)
    quotient, remainder = poly_divmod(q_coeffs[:], e_coeffs[:], p)
    coeffs = (quotient + [0] * k)[:k]
    print(' '.join(map(str, coeffs)))

main()

計算量: ガウスの消去法が $O(N^3)$(未知数 $2e+K=N$)。多項式除算は $O(N^2)$。全体 $O(N^3)$、$N\le200$ に十分高速。

Step-by-Step 解説

Step 1: 誤り位置検出多項式という発想

誤り位置集合 $T$ の根を持つ $E(x)=\prod_{t\in T}(x-t)$ を考えれば、$i\in T$でも$i\notin T$でも統一的に $y_iE(x_i)=Q(x_i)$ が成立する。

Step 2: 線形方程式として解く

$E$の次数$e=\lfloor(N-K)/2\rfloor$、$Q$の次数$e+K-1$とすると未知数はちょうど$N$個になり、$N$点の観測から一意に解ける。

Step 3: $f=Q/E$ が割り切れる理由

条件結果
実際の誤り数 $\le e$$Q/E$の余りが厳密に0、$f$が復元できる

Step 4: 出力の整形

商の次数が$K-1$未満の場合はゼロ埋めして長さ$K$に揃える。

よくあるミス

ミス原因正しい書き方
$E(x)$の最高次係数を未知数に含めるモニック制約を忘れる$e$次は定数1として式に組み込む
mod pの逆元を取らず通常除算有限体上の計算を忘れるpow(x, p-2, p)で逆元を掛ける
$e$を実際の誤り数ぴったりに設定しようとする誤り数は未知常に上界$\lfloor(N-K)/2\rfloor$を使う

次のステップ

  • 発展問題: 誤り数が未知のまま$e$を0から増やし最初に余り0になる$e$を探索する
  • 発展問題: Berlekamp-Massey + Chien探索 + Forneyアルゴリズムとの計算量比較

自己評価

理解度: / /

自分の回答:

気づき・メモ: