Day 124-Q5 — Berlekamp's Q-Matrix Algorithm（有限体上の多項式因数分解）

問題

素数 $p$ と、有限体 $\mathrm{GF}(p)$ 上の $n$ 次多項式 $f(x)$（重複因子を持たないことが保証される、すなわち $f$ とその形式微分 $f'$ が互いに素）が与えられる。$f(x)$ を $\mathrm{GF}(p)$ 上で既約な多項式の積に因数分解し、既約因子をすべて出力せよ（順序は次数昇順、同次数内は係数の辞書順で構わない）。

入力形式

p n
c_0 c_1 ... c_n

制約

入出力例

5 4
1 0 0 0 1

2
2 2 0 1
2 3 0 1

概念図

ヒント（段階的開示）

factors = [f]
for basis_vec in null_space_basis:  # 定数ベクトル以外
    if len(factors) == k:
        break
    new_factors = []
    for g in factors:
        if degree(g) == 1:
            new_factors.append(g)
            continue
        # 1つの s で打ち切らず、s=0..p-1 全部を元の g に対して試す。
        # 見つかった gcd 達は互いに素で、掛け合わせると g に戻る性質を利用し、
        # 1回の基底ベクトルで可能な限り深く分解する。
        pieces, deg_sum = [], 0
        for s in range(p):
            h = poly_gcd(g, poly_sub(basis_vec, [s], p), p)
            if degree(h) > 0:
                pieces.append(h)
                deg_sum += degree(h)
                if deg_sum == degree(g):
                    break
        new_factors.extend(pieces if len(pieces) >= 2 else [g])
    factors = new_factors

模範解答 (Python)

📋 コピーimport sys

def solve():
    data = sys.stdin.read().split()
    idx = 0
    p = int(data[idx]); idx += 1
    n = int(data[idx]); idx += 1
    f = [int(data[idx + i]) for i in range(n + 1)]
    idx += n + 1

    def trim(a):
        a = a[:]
        while len(a) > 1 and a[-1] == 0:
            a.pop()
        return a

    def padd(a, b):
        m = max(len(a), len(b))
        r = [0] * m
        for i, x in enumerate(a): r[i] = (r[i] + x) % p
        for i, x in enumerate(b): r[i] = (r[i] + x) % p
        return trim(r)

    def psub(a, b):
        m = max(len(a), len(b))
        r = [0] * m
        for i, x in enumerate(a): r[i] = (r[i] + x) % p
        for i, x in enumerate(b): r[i] = (r[i] - x) % p
        return trim(r)

    def pmul(a, b):
        r = [0] * (len(a) + len(b) - 1)
        for i, ai in enumerate(a):
            if ai == 0: continue
            for j, bj in enumerate(b):
                r[i + j] = (r[i + j] + ai * bj) % p
        return trim(r)

    def pdivmod(a, b):
        a = trim(a); b = trim(b)
        if len(a) < len(b):
            return [0], a
        q = [0] * (len(a) - len(b) + 1)
        r = a[:]
        inv_lead = pow(b[-1], p - 2, p)
        while len(trim(r)) >= len(b) and any(r):
            r = trim(r)
            if len(r) < len(b):
                break
            coef = (r[-1] * inv_lead) % p
            shift = len(r) - len(b)
            q[shift] = coef
            for i, bi in enumerate(b):
                r[shift + i] = (r[shift + i] - coef * bi) % p
            r = trim(r)
        return trim(q), trim(r)

    def pmod(a, mod):
        return pdivmod(a, mod)[1]

    def pgcd(a, b):
        a, b = trim(a), trim(b)
        while not (len(b) == 1 and b[0] == 0):
            _, r = pdivmod(a, b)
            a, b = b, (r if r else [0])
        inv = pow(a[-1], p - 2, p)
        return trim([c * inv % p for c in a])

    def degree(a):
        a = trim(a)
        return len(a) - 1 if not (len(a) == 1 and a[0] == 0) else -1

    # --- x^p mod f を高速べき乗で求める ---
    def ppow_mod(e, mod):
        result = [1]
        base = [0, 1]
        while e:
            if e & 1:
                result = pmod(pmul(result, base), mod)
            base = pmod(pmul(base, base), mod)
            e >>= 1
        return result

    xp = ppow_mod(p, f)

    # --- Q行列構築: 行i = x^(i*p) mod f の係数（長さnにパディング） ---
    Q = [[0] * n for _ in range(n)]
    cur = [1]
    for i in range(n):
        for j, c in enumerate(cur):
            if j < n:
                Q[i][j] = c
        cur = pmod(pmul(cur, xp), f)

    # Q - I
    for i in range(n):
        Q[i][i] = (Q[i][i] - 1) % p

    # --- ガウス消去で Q-I の核空間を求める ---
    mat = [row[:] for row in Q]
    pivot_cols = []
    row_i = 0
    for col in range(n):
        piv = -1
        for r in range(row_i, n):
            if mat[r][col] != 0:
                piv = r
                break
        if piv == -1:
            continue
        mat[row_i], mat[piv] = mat[piv], mat[row_i]
        inv = pow(mat[row_i][col], p - 2, p)
        mat[row_i] = [(x * inv) % p for x in mat[row_i]]
        for r in range(n):
            if r != row_i and mat[r][col] != 0:
                factor = mat[r][col]
                mat[r] = [(mat[r][c2] - factor * mat[row_i][c2]) % p for c2 in range(n)]
        pivot_cols.append(col)
        row_i += 1

    free_cols = [c for c in range(n) if c not in pivot_cols]
    k = len(free_cols)  # 既約因子の個数

    basis = []
    for fc in free_cols:
        vec = [0] * n
        vec[fc] = 1
        r = 0
        for pc in pivot_cols:
            vec[pc] = (-mat[r][fc]) % p
            r += 1
        basis.append(trim(vec))

    factors = [trim(f)]
    if k > 1:
        for bv in basis:
            if degree(bv) <= 0:
                continue  # 定数ベクトルは分解に使えない
            if len(factors) >= k:
                break
            new_factors = []
            for g in factors:
                if degree(g) == 1:
                    new_factors.append(g)
                    continue
                # s=0..p-1 それぞれについて gcd(g, v-s) を元の g に対して計算する。
                # v は各既約因子上で異なる定数値を取るため、こうして集めた gcd 達は
                # 互いに素であり、掛け合わせると g に戻る（中国剰余定理の帰結）。
                # 1つの s で打ち切らず全 s を試すことで、1回の基底ベクトルで
                # 最大限分解できる。
                pieces = []
                deg_sum = 0
                for s in range(p):
                    h = pgcd(g, psub(bv, [s]))
                    if degree(h) > 0:
                        pieces.append(h)
                        deg_sum += degree(h)
                        if deg_sum == degree(g):
                            break
                if len(pieces) >= 2:
                    new_factors.extend(pieces)
                else:
                    new_factors.append(g)
            factors = new_factors

    factors = [trim(fa) for fa in factors]
    factors.sort(key=lambda a: (degree(a), a))

    print(len(factors))
    for fa in factors:
        print(degree(fa), *fa)

solve()

Step-by-Step 解説

よくあるミス

次のステップ

• 発展: $p$ が大きい場合（$p > n$ 程度)は $s$ の全探索が非効率になるため、Cantor-Zassenhaus法（ランダム多項式による確率的分割）を組み合わせるのが実用的。重複因子がある一般の入力には square-free分解（$f/\gcd(f,f')$ の繰り返し）を前処理として組み合わせる。
• 次回予告: Master Levelローテーション継続（次サイクルの5テーマ）

自己評価

自分の回答

気づき・メモ