Day 140-Q4 — Berlekamp-Welch Algorithm（Reed-Solomon誤り訂正復号）

問題

素数 $p$、メッセージ長 $k$、訂正可能誤り数 $t$ が与えられる。

送信側は次数$k-1$以下の多項式 $P(x)=c_0+c_1x+\dots+c_{k-1}x^{k-1}\pmod p$ を、評価点 $x_i=i$（$i=0,\dots,n-1$、$n=k+2t$）における値として符号化し送信する（Reed-Solomon符号）。

受信側では最大$t$個までの値が改ざんされた符号語 $y_0,\dots,y_{n-1}$ が届く。Berlekamp-Welch アルゴリズムを用いて元の係数 $c_0,\dots,c_{k-1}$ を復元せよ。

手順:

1. $Q(x)$（次数$k+t-1$以下）と、モニックな誤り位置多項式$E(x)$（次数$t$）を未知数として用意
2. 各$i$について $Q(x_i)=y_i\cdot E(x_i)\pmod p$ の連立方程式（未知数$n$個、方程式$n$本）を解く
3. $Q(x)$を$E(x)$で多項式除算し、割り切れれば商が$P(x)$。割り切れない/連立方程式が解けない場合は訂正不能

訂正不能な場合は-1を出力せよ。

入力形式

p k t
y_0 y_1 ... y_{n-1}

制約

入出力例

17 3 2
2 6 16 3 13 15 5

2 3 1

17 3 2
2 0 1 3 13 2 5

-1

概念図

ヒント（段階的開示）

# 未知数: q_0..q_{k+t-1} (Q), e_0..e_{t-1} (E, 最高次はモニックで1固定)
# 方程式: sum q_j x_i^j - y_i*sum e_j x_i^j = y_i * x_i^t  (i=0..n-1)
A = []
b = []
for i in range(n):
    row = [pow(x[i], j, p) for j in range(k+t)] + [(-y[i]*pow(x[i], j, p)) % p for j in range(t)]
    A.append(row)
    b.append((y[i] * pow(x[i], t, p)) % p)

sol = gauss_solve(A, b, p)
if sol is None:
    print(-1)
else:
    Q = sol[:k+t]
    E = sol[k+t:] + [1]
    quot, rem = poly_divmod(Q, E, p)
    if 余りが全て0でない:
        print(-1)
    else:
        print(quot[:k])

模範解答 (Python)

📋 コピーimport sys

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

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

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

    def poly_divmod(num, den):
        num = num[:]
        dn = len(den) - 1
        while dn > 0 and den[dn] == 0:
            dn -= 1
        quot = [0] * max(1, len(num) - dn)
        while True:
            deg_num = len(num) - 1
            while deg_num >= 0 and num[deg_num] % p == 0:
                deg_num -= 1
            if deg_num < dn:
                break
            coef = (num[deg_num] * modinv(den[dn])) % p
            shift = deg_num - dn
            quot[shift] = coef
            for i in range(dn + 1):
                num[shift + i] = (num[shift + i] - coef * den[i]) % p
        return quot, num

    A = []
    b = []
    for i in range(n):
        xi = xs[i] % p
        row = []
        for j in range(k + t):
            row.append(pow(xi, j, p))
        for j in range(t):
            row.append((-ys[i] * pow(xi, j, p)) % p)
        A.append(row)
        b.append((ys[i] * pow(xi, t, p)) % p)

    sol = gauss_solve(A, b)
    if sol is None:
        print(-1)
        return

    Qc = sol[:k + t]
    Ec = sol[k + t:] + [1]
    quot, rem = poly_divmod(Qc, Ec)
    if any(x % p for x in rem):
        print(-1)
        return

    while len(quot) < k:
        quot.append(0)
    print(' '.join(str(x % p) for x in quot[:k]))

main()

Step-by-Step 解説

よくあるミス

次のステップ

• 発展: 実際のReed-Solomon復号ではBerlekamp-Massey法（$O(n^2)$）を使う「シンドローム復号」がより高速。両者の関係と高速化の理由を調べてまとめよ。
• 発展2: 評価点を任意の相異なる元にした場合でも動くことを確認し、$n>p$になるケースを弾く入力検証を追加せよ。

自己評価

自分の回答

気づき・メモ