掃き出し法 (row reduction), ガウスの消去法 (Gaussian Elimination)

概要

連立方程式 \(A \boldsymbol{x} = \boldsymbol{b}\) を解くためのアルゴリズム。

拡大係数行列を行基本変形で上三角の形に変形し、唯一解、解なし、いくつも解ありのいずれかかを判別する。

計算量

\(O(N^3)\)

実装

EPS = 1e-9
def gaussian_elimination(N, A, bs):
    M = [[0]*(N+1) for i in range(N)]
    for i in range(N):
        M[i][:N] = A[i]
        M[i][N] = bs[i]
    cur = 0
    for i in range(N):
        k = cur; ma = abs(M[cur][i])
        for j in range(cur+1, N):
            if ma < abs(M[j][i]):
                k = j
                ma = abs(M[j][i])
        if ma < EPS:
            continue
        if k != cur:
            M[cur], M[k] = M[k], M[cur]
        Mc = M[cur]
        for Mj in M[cur+1:]:
            if abs(Mj[i]) > EPS:
                e = Mj[i] / Mc[i]
                Mj[i] = 0
                for k in range(i+1, N+1):
                    Mj[k] -= e*Mc[k]
        cur += 1
    if cur == N:
        xs = [Mi[N] for Mi in M]
        for i in range(N-1, -1, -1):
            Mi = M[i]
            xs[i] -= sum(Mi[j] * xs[j] for j in range(i+1, N))
            xs[i] /= Mi[i]
        return 1, xs
    if any(Mi[N] > EPS for Mi in M[cur:]):
        return 0, None
    return 2, None

Verified

  • AOJ: "1232: Super Star": source (Python3, 7.62sec)

実装 (logical matrix)

from operator import xor
# calculate a vector x s.t. Ax = b
# - A is a N×N logical matrix
# - x and b are N-dimensional logical vectors
def gaussian_elimination(N, A, bs):
    M = [[0]*(N+1) for i in range(N)]
    for i in range(N):
        M[i][:N] = A[i]
        M[i][N] = bs[i]
    cur = 0
    for i in range(N):
        if M[cur][i] == 0:
            k = cur+1
            while k < N and M[k][i] == 0:
                k += 1
            if k == N:
                continue
            M[cur], M[k] = M[k], M[cur]
        Mc = M[cur]
        for Mj in (Mj for Mj in M[cur+1:] if Mj[i]):
            Mj[i:] = map(xor, Mj[i:], Mc[i:])
        cur += 1
    if cur == N:
        # an unique solution
        xs = [Mi[N] for Mi in M]
        for i in range(N-1, -1, -1):
            Mi = M[i]
            xs[i] ^= sum(Mi[j] & xs[j] for j in range(i+1, N)) & 1
        return 1, xs
    if any(Mi[N] for Mi in M[cur:]):
        # no solution
        return 0, None
    # infinitely many solutions
    return 2, None

N = 3
A = [[1, 0, 1], [0, 1, 1], [1, 0, 0]]
bs = [1, 0, 1]
print(gaussian_elimination(N, A, bs))
# => "(1, [1, 0, 0])"

Verified

  • AOJ: "2624: Graph Automata Player": source (Python3, 36.71sec)

参考

戻る