Day 079-Q2 — 多項式行列式 mod p（Hessenberg 変換 + Bareiss アルゴリズム）

2026-07-01 赤色 Master / Phase 8+ ★★★★★★★★★ 行列式・Bareiss・Gaussian elimination・mod p

問題

$N \times N$ の整数行列 $A$ が与えられる。素数 $p$ に対して $\det(A) \bmod p$ を求めよ。

ここで $\det(A)$ は行列 $A$ の行列式（determinant）である。 通常の Gaussian elimination では各ステップで割り算が発生して分数が現れるが、 Bareiss アルゴリズムを使うことで整数演算のみで行列式を計算できる。 mod $p$ 上ではモジュラー逆元を使った Gaussian elimination が実装上シンプル。

入力の全要素は既に $0 \le A_{i,j} < p$ の範囲で与えられる。

入力形式

N p
A_{1,1} A_{1,2} ... A_{1,N}
A_{2,1} A_{2,2} ... A_{2,N}
...
A_{N,1} A_{N,2} ... A_{N,N}

制約

パラメータ範囲備考
$N$$1 \le N \le 500$行列のサイズ
$p$$2 \le p \le 10^9 + 7$（素数）modulus（素数保証）
$A_{i,j}$$0 \le A_{i,j} < p$入力は既に mod p 済み

入出力例

入力例1（3×3）

3 1000000007
1 2 3
4 5 6
7 8 10

出力例1

1000000004

$\det(A) = -3$。$-3 \bmod (10^9+7) = 10^9+7-3 = 1000000004$。 展開: $1(50-48) - 2(40-42) + 3(32-35) = 2+4-9 = -3$。

入力例2（4×4）

4 998244353
2 1 3 4
0 3 1 2
5 2 4 1
1 0 2 3

出力例2

25

Bareiss アルゴリズムで計算すると $\det = 25$（整数値）。$25 \bmod 998244353 = 25$。

入力例3（特異行列）

3 7
1 2 3
4 5 6
7 8 9

出力例3

0

行3 = 行1 + 2×行2 なので $\det = 0$（特異行列）。ピボットが見つからない列が存在 → 0 を返す。

概念図: Gaussian Elimination のピボット操作

Gaussian Elimination mod p — ステップ別ピボット消去 元の行列 A col1 col2 col3 r1 1 2 3 r2 4 5 6 r3 7 8 10 k=0 k=0 消去後 1 2 3 r1 r2←r2-4r1 0 -3 -6 r3←r3-7r1 0 -6 -11 k=1 k=1 消去後（上三角） 1 2 3 0 -3 -6 r3←r3-2r2 0 0 1 行列式の計算: 対角積 × 符号 上三角行列の対角要素: [1, -3, 1] det = sign × A[0][0] × A[1][1] × A[2][2] = (+1) × 1 × (-3) × 1 = -3 mod p: -3 mod (10^9+7) = 10^9+7-3 = 1000000004 modinv(a, p) = pow(a, p-2, p) ← Fermat の小定理 Bareiss 漸化式: A[i][j] = (A[k][k]×A[i][j] - A[i][k]×A[k][j]) / A_prev_pivot

Gaussian elimination: 各列 $k$ でピボット（赤枠）を選び、下の行を消去して上三角行列を作る。 対角積に符号（行交換の回数）を掛けたものが行列式。mod $p$ 上では除算をモジュラー逆元に置き換える。

ヒント

ヒント1（方向性）

$\det(A)$ を素朴な Gaussian elimination で求めると、各ステップで割り算が発生して分数になる。 $p$ が素数ならモジュラー逆元 pow(x, p-2, p) で割り算を $O(\log p)$ で代替できるが、 ピボットが 0 の場合に逆元が存在しないため列（または行）交換が必要。

Bareiss アルゴリズムは「前のステップのピボットで割る」ことで分母を管理し、 整数行列なら整数演算だけで行列式を求める。mod $p$ 上では割り算をモジュラー逆元に変えて使える。

ヒント2（アプローチ）

方法A: Gaussian elimination mod p（推奨）

  1. 各列 $k$ でピボット行を見つける（A[k][k] != 0 を探して行交換）
  2. 行交換のたびに sign *= -1
  3. 行 $i > k$ に対して: A[i] -= A[i][k] * modinv(A[k][k]) * A[k]（mod p）
  4. 対角積 $\prod_{k} A[k][k]$ に符号を掛けて答え

方法B: Bareiss アルゴリズム

漸化式: $$A^{(k)}_{i,j} = \frac{A^{(k-1)}_{k,k} \cdot A^{(k-1)}_{i,j} - A^{(k-1)}_{i,k} \cdot A^{(k-1)}_{k,j}}{A^{(k-2)}_{k-1,k-1}}$$ mod $p$ では分母をモジュラー逆元で処理。

ヒント3（コードスケルトン） 
def modinv(a, p):
    # Fermat の小定理: a^(p-1) ≡ 1 (mod p) → a^(-1) ≡ a^(p-2) (mod p)
    return pow(a, p - 2, p)

def det_mod(A, n, p):
    M = [row[:] for row in A]   # コピー（元を破壊しない）
    sign = 1
    for k in range(n):
        # ピボット選択: A[k][k] != 0 となる行を探す
        pivot = -1
        for i in range(k, n):
            if M[i][k] != 0:
                pivot = i
                break
        if pivot == -1:
            return 0   # 特異行列: det = 0
        if pivot != k:
            M[k], M[pivot] = M[pivot], M[k]
            sign = -sign   # 行交換で符号反転
        # 前進消去
        inv_piv = modinv(M[k][k], p)
        for i in range(k + 1, n):
            factor = M[i][k] * inv_piv % p
            for j in range(k, n):
                M[i][j] = (M[i][j] - factor * M[k][j]) % p
    # 対角積
    result = sign % p
    for k in range(n):
        result = result * M[k][k] % p
    return result % p

模範解答

方法A: Gaussian elimination mod p（推奨・シンプル）

import sys

def solve():
    data = sys.stdin.buffer.read().split()
    idx = 0
    N, MOD = int(data[idx]), int(data[idx + 1])
    idx += 2

    A = []
    for i in range(N):
        row = [int(data[idx + j]) % MOD for j in range(N)]
        A.append(row)
        idx += N

    print(det_mod_gaussian(A, N, MOD))


def modinv(a, m):
    """Fermat の小定理による逆元（m は素数であること）"""
    return pow(a, m - 2, m)


def det_mod_gaussian(A, n, mod):
    """
    Gaussian elimination mod p で行列式を計算。
    時間計算量: O(N^3 * log(mod))    ← modinv 呼び出しが O(N^2) 回
    ただし pow() の定数が小さいため実測は O(N^3) に近い
    空間計算量: O(N^2)
    """
    M = [row[:] for row in A]   # 元の行列を破壊しないためコピー
    sign = 1                     # 行交換の符号

    for k in range(n):
        # ─── ピボット選択: 列 k で最初の非ゼロ要素を持つ行を探す ─
        pivot_row = -1
        for i in range(k, n):
            if M[i][k] != 0:
                pivot_row = i
                break

        if pivot_row == -1:
            return 0    # 列 k が全て 0 → 特異行列 → det = 0

        # 行交換（必要な場合のみ）
        if pivot_row != k:
            M[k], M[pivot_row] = M[pivot_row], M[k]
            sign = -sign    # 行交換で行列式の符号が反転

        # ─── 前進消去: ピボット行より下の行を消去 ───────────────
        inv_pivot = modinv(M[k][k], mod)
        for i in range(k + 1, n):
            if M[i][k] == 0:
                continue    # 既に 0 → スキップで定数倍高速化
            factor = M[i][k] * inv_pivot % mod
            for j in range(k, n):
                M[i][j] = (M[i][j] - factor * M[k][j]) % mod

    # ─── 上三角行列の対角積が行列式（符号も掛ける） ────────────
    result = sign % mod    # sign は -1 または 1
    for k in range(n):
        result = result * M[k][k] % mod

    return result % mod    # Python の % は常に非負


solve()

方法B: Bareiss アルゴリズム mod p（分数フリー Gaussian elimination）

import sys

def solve_bareiss():
    data = sys.stdin.buffer.read().split()
    idx = 0
    N, MOD = int(data[idx]), int(data[idx + 1])
    idx += 2

    A = []
    for i in range(N):
        row = [int(data[idx + j]) % MOD for j in range(N)]
        A.append(row)
        idx += N

    print(det_mod_bareiss(A, N, MOD))


def det_mod_bareiss(A, n, mod):
    """
    Bareiss アルゴリズム（分数フリー Gaussian elimination）の mod p 版。

    更新式:
        M[i][j] = (M[k][k] * M[i][j] - M[i][k] * M[k][j]) / M_prev_pivot
    整数版では分母が必ず割り切れる（Sylvester の恒等式による被除性）。
    mod p 版では / をモジュラー逆元に変換。

    時間計算量: O(N^3)
    空間計算量: O(N^2)
    """
    M = [row[:] for row in A]
    sign = 1
    prev_pivot = 1    # A^{(k-1)}_{k-1,k-1}（初期値 = 1 = 単位元）

    for k in range(n):
        # ─── ピボット選択 ────────────────────────────────────────
        pivot_row = -1
        for i in range(k, n):
            if M[i][k] != 0:
                pivot_row = i
                break
        if pivot_row == -1:
            return 0
        if pivot_row != k:
            M[k], M[pivot_row] = M[pivot_row], M[k]
            sign = -sign

        cur_pivot = M[k][k]
        # 前ステップのピボットの逆元（prev_pivot = 1 のときは不要）
        inv_prev = pow(prev_pivot, mod - 2, mod) if prev_pivot != 1 else 1

        # ─── Bareiss の更新 ──────────────────────────────────────
        for i in range(k + 1, n):
            for j in range(k + 1, n):
                # (cur_pivot * M[i][j] - M[i][k] * M[k][j]) / prev_pivot
                val = (cur_pivot * M[i][j] - M[i][k] * M[k][j]) % mod
                M[i][j] = val * inv_prev % mod
            M[i][k] = 0    # 消去済みに設定（省略可能だが明示）

        prev_pivot = cur_pivot

    result = sign % mod
    for k in range(n):
        result = result * M[k][k] % mod
    return result % mod


solve_bareiss()

Hessenberg 変換を使ったアプローチ（参考）

def to_hessenberg_det(A, n, mod):
    """
    相似変換で上 Hessenberg 形に変換後、O(N^2) で行列式を計算。
    Hessenberg 形: H[i][j] = 0 for i > j + 1
    行列式は漸化式で計算: D[k] = H[k-1][k-1]*D[k-1] - H[k-1][k-2]*H[k-2][k-1]*D[k-2] + ...
    """
    M = [row[:] for row in A]

    # 列 k で k+2 行目以下を消去（相似変換: 行操作 + 対応する列操作）
    for k in range(n - 2):
        # サブ対角要素 M[k+1][k] を非ゼロにする
        pivot = -1
        for i in range(k + 1, n):
            if M[i][k] != 0:
                pivot = i
                break
        if pivot == -1:
            continue
        if pivot != k + 1:
            # 行交換
            M[k + 1], M[pivot] = M[pivot], M[k + 1]
            # 対応する列交換（相似変換: det 不変）
            for i in range(n):
                M[i][k + 1], M[i][pivot] = M[i][pivot], M[i][k + 1]

        inv_sub = pow(M[k + 1][k], mod - 2, mod)
        for i in range(k + 2, n):
            if M[i][k] == 0:
                continue
            factor = M[i][k] * inv_sub % mod
            # 行 i から factor * 行(k+1) を引く
            for j in range(n):
                M[i][j] = (M[i][j] - factor * M[k + 1][j]) % mod
            # 列 (k+1) に factor * 列 i を足す（行列式不変）
            for j in range(n):
                M[j][k + 1] = (M[j][k + 1] + factor * M[j][i]) % mod

    # 上 Hessenberg 行列の行列式（漸化式）
    D = [0] * (n + 1)
    D[0] = 1
    for k in range(1, n + 1):
        D[k] = M[k - 1][k - 1] * D[k - 1] % mod
        sub_diag_prod = 1
        for i in range(k - 2, -1, -1):
            sub_diag_prod = sub_diag_prod * M[i + 1][i] % mod
            D[k] = (D[k] - M[k - 1][i] * sub_diag_prod % mod * D[i]) % mod

    return D[n] % mod

Step-by-Step 解説

Step 1: 行列式の基本性質と Gaussian elimination の根拠

行列式には以下の基本変換が成り立つ（多重線形性・交代性）:

操作行列式への影響使用場面
行 $i$ を $c$ 倍$\det$ も $c$ 倍対角正規化
行 $i$ に $c \times$ 行 $j$ を加算（$i \ne j$）$\det$ は変わらない前進消去の核心
行 $i$ と行 $j$ を交換$\det$ の符号が反転ピボット選択

Gaussian elimination は性質2と3を使って上三角行列を作り、対角積を求める。

Step 2: mod p での除算とモジュラー逆元

mod $p$ 上での除算（$a/b$）は Fermat の小定理で実装する:

$$b^{-1} \equiv b^{p-2} \pmod{p}$$

これが成り立つのは $p$ が素数のとき（$b \not\equiv 0$）。Python では pow(b, p-2, p) が最速（内部で繰り返し二乗法）。

注意: ピボットが $0 \pmod{p}$ のとき逆元が存在しない。行交換で非ゼロ要素を持つ行を持ってくる必要がある。全ての行で 0 なら特異行列 → $\det = 0$。

Step 3: Bareiss アルゴリズムの仕組み

通常の Gaussian では各ステップで逆元計算（$O(\log p)$）が $O(N^2)$ 回発生する。Bareiss の改善点: 各ステップで「前のピボット $A^{(k-1)}_{k-1,k-1}$」で割ることで、整数行列なら整数のまま計算を進められる（Sylvester の恒等式による被除性の保証）。

更新式:

$$A^{(k)}_{i,j} = \frac{A^{(k-1)}_{k,k} \cdot A^{(k-1)}_{i,j} - A^{(k-1)}_{i,k} \cdot A^{(k-1)}_{k,j}}{A^{(k-1)}_{k-1,k-1}}$$

  • 整数版: 分母が必ず割り切れる → 純粋整数演算（多倍長整数注意）
  • mod p 版: 分母をモジュラー逆元で処理（ただし逆元計算コストは同じ）

Step 4: Hessenberg 変換の概念

上 Hessenberg 形とは $H[i][j] = 0$ for $i > j+1$ の行列（対角の 1 つ下のサブ対角には非ゼロが残る）:

上 Hessenberg 形:
  * * * *
  * * * *
  0 * * *
  0 0 * *

変換方法: 相似変換 $H = P^{-1}AP$ を使うことで $\det(H) = \det(A)$ を保ちつつ、列 $k$ の $k+2$ 行目以下を消去。消去するときに対応する列操作を同時に施す（相似変換なので $\det$ は不変）。Hessenberg 行列の行列式は展開式の漸化式で $O(N^2)$ だが、実装が複雑なため競技プログラミングでは方法A が主流。

Step 5: 計算量まとめ

手法時間計算量空間計算量実装難易度
Gaussian elimination mod p$O(N^3 \log p)$$O(N^2)$容易
Bareiss mod p$O(N^3 \log p)$$O(N^2)$やや難
Hessenberg + 行列式$O(N^3)$ 変換 + $O(N^2)$$O(N^2)$難

$N=500$ で $N^3 = 1.25\times10^8$。modinv は pow で高速化済み。実際は $O(N^3)$ 程度で動作。PyPy 推奨。

よくあるミス

ミス原因正しい書き方
modinv(0, p) を呼んでしまう ピボットが 0 のときに逆元を求めようとする ピボット選択で A[k][k] != 0 を確認してから modinv を呼ぶ
行交換を忘れる A[k][k] = 0 でそのまま進むと 0 除算エラー 必ずピボット選択ループで非ゼロ行を探す
符号 sign を mod しない sign = -1 のまま最終計算で負になる result * sign % mod で Python の % が常に非負を保証
Bareiss で prev_pivot = 0 特異行列のケースで 0 除算 Bareiss も行交換 + ピボットチェックが必要
Hessenberg の列操作を忘れる 行だけ消去すると行列式が変わる 行消去のたびに対応する列操作を必ず施す（相似変換）
$N=1$ の特殊ケース ループが 0 回でも result = A[0][0] で正しい 特別処理不要。N=1 テストは確認する

次のステップ

  • 発展問題1: 行列式ではなく永久式（Permanent）を求めよ。行列式と異なり $O(N!)$ か Ryser の公式で $O(2^N N)$ 必要。
  • 発展問題2: $\mathbb{F}_p[x]$ 上の多項式行列式（各要素が $x$ の多項式）を求めよ。結果は $x$ の多項式になり、係数ごとに mod p 計算が必要。
  • 発展問題3: Sylvester 行列の行列式（= 2つの多項式の終結式 resultant）を $O(N^3)$ で計算し、多項式 GCD の次数を求めよ。
  • AtCoder 関連問題: 行列木定理（Kirchhoff のマトリクス木定理）、行列累乗の行列式（det の $K$ 乗 = det(A^K)）。

自己評価