Day 058-Q5 — 特性多項式・Cayley-Hamilton・行列ランク

問題

$N \times N$ の整数行列 $A$（mod $p$）が与えられる。以下を計算せよ：

1. $A$ の特性多項式 $\chi_A(\lambda) = \det(\lambda I - A) \pmod{p}$ の係数列
2. Cayley-Hamilton 定理を用いて $A^K \pmod{p}$ を計算せよ（$K \le 10^{18}$）
3. $A$ の rank（行列のランク）を出力せよ

入出力例

3 5 1000000007
1 2 3
4 5 6
7 8 9

特性多項式: [1, ...]
A^5: (各要素 mod p)
rank: 2

$\chi_A(\lambda) = \lambda^3 - 15\lambda^2 - 18\lambda$（係数は mod p）。
rank = 2（行列式 = 0）。$A^K$ は Cayley-Hamilton により $\{I, A, A^2\}$ の線形結合。

概念図: Cayley-Hamilton を使った A^K の高速計算

ヒント（段階的開示）

def poly_mul_rem(a, b, charpoly, n, mod):
    """多項式 a * b を charpoly で剰余"""
    c = [0] * (2*n - 1)
    for i in range(len(a)):
        for j in range(len(b)):
            c[i+j] = (c[i+j] + a[i] * b[j]) % mod
    # λ^n = -(c_{n-1}λ^{n-1} + ... + c_0)
    for i in range(len(c)-1, n-1, -1):
        if c[i] == 0:
            continue
        for j in range(n):
            c[i-n+j] = (c[i-n+j] - c[i] * charpoly[j+1]) % mod
        c[i] = 0
    return [x % mod for x in c[:n]]

# λ^K mod charpoly を繰り返し2乗法で計算
result = [1] + [0]*(n-1)   # = 1
base   = [0, 1] + [0]*(n-2) # = λ

模範解答 (Python)

コピーimport sys
input = sys.stdin.readline

def mat_mul(A, B, mod):
    n = len(A)
    C = [[0]*n for _ in range(n)]
    for i in range(n):
        for k in range(n):
            if A[i][k] == 0: continue
            for j in range(n):
                C[i][j] = (C[i][j] + A[i][k] * B[k][j]) % mod
    return C

def gauss_mod(A, mod):
    n, m = len(A), len(A[0])
    rank = 0
    for col in range(m):
        pivot = next((row for row in range(rank, n) if A[row][col] != 0), -1)
        if pivot == -1: continue
        A[rank], A[pivot] = A[pivot], A[rank]
        inv = pow(A[rank][col], mod-2, mod)
        for j in range(m): A[rank][j] = A[rank][j] * inv % mod
        for row in range(n):
            if row != rank and A[row][col] != 0:
                f = A[row][col]
                for j in range(m): A[row][j] = (A[row][j] - f * A[rank][j]) % mod
        rank += 1
    return rank

def characteristic_poly(A, mod):
    n = len(A)
    H = [row[:] for row in A]
    for j in range(n-2):
        for i in range(j+2, n):
            if H[i][j] == 0: continue
            H[i], H[j+1] = H[j+1], H[i]
            for k in range(n): H[k][i], H[k][j+1] = H[k][j+1], H[k][i]
            break
        if H[j+1][j] == 0: continue
        inv = pow(H[j+1][j], mod-2, mod)
        for i in range(j+2, n):
            if H[i][j] == 0: continue
            f = H[i][j] * inv % mod
            for k in range(n): H[i][k] = (H[i][k] - f * H[j+1][k]) % mod
            for k in range(n): H[k][j+1] = (H[k][j+1] + f * H[k][i]) % mod
    p = [[1], [(-H[0][0]) % mod, 1]]
    for k in range(1, n):
        prev = p[k]
        new_p = [0] * (len(prev) + 1)
        for i, v in enumerate(prev):
            new_p[i+1] = (new_p[i+1] + v) % mod
            new_p[i] = (new_p[i] - H[k][k] * v) % mod
        coef = 1
        for m in range(k-1, -1, -1):
            coef = coef * H[m+1][m] % mod
            hm = H[k][m]
            for i, v in enumerate(p[m]):
                new_p[i] = (new_p[i] - coef * hm % mod * v) % mod
        p.append([x % mod for x in new_p])
    return [x % mod for x in p[n]]

def poly_mod_pow(exp, charpoly, mod):
    n = len(charpoly) - 1
    def mul(a, b):
        c = [0] * (2*n - 1)
        for i in range(len(a)):
            for j in range(len(b)):
                if i+j < len(c): c[i+j] = (c[i+j] + a[i] * b[j]) % mod
        for i in range(len(c)-1, n-1, -1):
            if c[i] == 0: continue
            ci = c[i]
            for j in range(n): c[i-n+j] = (c[i-n+j] - ci * charpoly[j+1]) % mod
            c[i] = 0
        return [x % mod for x in c[:n]]
    result = [0]*n; result[0] = 1
    base = [0]*n
    if n > 1: base[1] = 1
    else: base[0] = 1
    while exp > 0:
        if exp & 1: result = mul(result, base)
        base = mul(base, base)
        exp >>= 1
    return result

def poly_eval_matrix(coeffs, A, mod):
    n = len(A)
    result = [[0]*n for _ in range(n)]
    power = [[1 if i==j else 0 for j in range(n)] for i in range(n)]
    for c in coeffs:
        if c != 0:
            for i in range(n):
                for j in range(n):
                    result[i][j] = (result[i][j] + c * power[i][j]) % mod
        power = mat_mul(power, A, mod)
    return result

def main():
    N, K, p = map(int, input().split())
    A = [list(map(int, input().split())) for _ in range(N)]
    charpoly = characteristic_poly(A, p)
    print("特性多項式:", charpoly)
    if K == 0:
        AK = [[1 if i==j else 0 for j in range(N)] for i in range(N)]
    else:
        coeffs = poly_mod_pow(K, charpoly, p)
        AK = poly_eval_matrix(coeffs, A, p)
    print(f"A^{K}:")
    for row in AK: print(*row)
    A_copy = [row[:] for row in A]
    print("rank:", gauss_mod(A_copy, p))

main()

Step-by-Step 解説

Step 1: 特性多項式（Hessenberg 法）

行列を上ヘッセンベルク型に変形（$O(N^3)$）後、再帰式で特性多項式を $O(N^2)$ で計算。全体 $O(N^3)$。

Step 2: Cayley-Hamilton 定理の応用

$\chi_A(A) = 0$ なので $A^N$ は $\{I, A, \ldots, A^{N-1}\}$ の線形結合。$\lambda^K \bmod \chi_A(\lambda)$ を多項式の繰り返し2乗法で $O(N^2 \log K)$ 計算し、その係数で $A^K$ を評価。

Step 3: mod $p$ でのランク

素体 $\mathbb{F}_p$ 上のガウス消去で $O(N^3)$。ピボット消去時に mod $p$ の逆元（フェルマーの小定理）を使用。

よくあるミス

次のステップ

発展問題: 複数の $K_1, K_2, \ldots, K_Q$ に対して $A^{K_i}$ を高速に計算せよ（$Q \le 10^5$, $K_i \le 10^{18}$）。

自己評価