Linear Equation (線形方程式)
(math/linear_equation.hpp)

• View this file on GitHub
• Last update: 2026-04-11 01:01:21+09:00
• Include: #include "math/linear_equation.hpp"

使い方

// 線形方程式 Ax = b を解き, {rank, x} を返す
auto [rank, x] = linear_equation(A, b);

参考文献

• けんちょんの競プロ精進記録, Gauss-Jordan の掃き出し法と、連立一次方程式の解き方
• 競プロ典型 90 問 057 - Flip Flap（★6）

Depends on

• Gauss-Jordan elimination (ガウスジョルダンの掃き出し法) (math/gauss_jordan_elimination.hpp)

Code

#pragma once

#include "math/gauss_jordan_elimination.hpp"

#include <cassert>
#include <utility>
#include <vector>

// 線形方程式 Ax = b を解く
// {rank, x} を返す
template <class T> std::pair<int, std::vector<T>> linear_equation(const std::vector<std::vector<T>>& A, const std::vector<T>& b) {
    assert(A.size() > 0);
    const int m = (int)(A.size()), n = (int)(A[0].size());
    std::vector M(m, std::vector<T>(n + 1));
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            M[i][j] = A[i][j];
        }
        M[i][n] = b[i];
    }
    int rank = gauss_jordan_elimination(M, n);
    show(M);

    // 解が存在しない場合, rank = -1 を返す
    for (int i = rank; i < m; i++) {
        if (M[i][n] != T(0)) {
            return {-1, {}};
        }
    }
    std::vector<T> x(n);
    for (int i = 0; i < rank; i++) x[i] = M[i][n];
    return {rank, x};
}

#line 2 "math/linear_equation.hpp"

#line 2 "math/gauss_jordan_elimination.hpp"

#include <cassert>
#include <vector>

// 行列 a を掃き出し, ランクを返す
template <class T> int gauss_jordan_elimination(std::vector<std::vector<T>>& a, int pivot_end = -1) {
    assert(a.size() > 0);
    const int h = (int)(a.size()), w = (int)(a[0].size());
    int rank = 0;
    if (pivot_end == -1) pivot_end = w;

    for (int j = 0; j < pivot_end; j++) {
        // ピボットを探す
        int pivot = -1;
        for (int i = rank; i < h; i++) {
            if (a[i][j] != T(0)) {
                pivot = i;
                break;
            }
        }
        if (pivot == -1) continue;

        std::swap(a[pivot], a[rank]);

        // a[rank][j] = 1 になるように a[rank] を定数倍
        const T inv = T(1) / a[rank][j];
        for (int j2 = 0; j2 < w; j2++) a[rank][j2] *= inv;

        // a[rank] を使って他の行を掃き出す
        for (int i = 0; i < h; i++) {
            if (i != rank and a[i][j] != T(0)) {
                const T fac = a[i][j];  // a[rank] を fac 倍して a[i] から引く
                for (int j2 = 0; j2 < w; j2++) {
                    a[i][j2] -= a[rank][j2] * fac;
                }
            }
        }
        rank++;
    }
    return rank;
}
#line 4 "math/linear_equation.hpp"

#line 6 "math/linear_equation.hpp"
#include <utility>
#line 8 "math/linear_equation.hpp"

// 線形方程式 Ax = b を解く
// {rank, x} を返す
template <class T> std::pair<int, std::vector<T>> linear_equation(const std::vector<std::vector<T>>& A, const std::vector<T>& b) {
    assert(A.size() > 0);
    const int m = (int)(A.size()), n = (int)(A[0].size());
    std::vector M(m, std::vector<T>(n + 1));
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            M[i][j] = A[i][j];
        }
        M[i][n] = b[i];
    }
    int rank = gauss_jordan_elimination(M, n);
    show(M);

    // 解が存在しない場合, rank = -1 を返す
    for (int i = rank; i < m; i++) {
        if (M[i][n] != T(0)) {
            return {-1, {}};
        }
    }
    std::vector<T> x(n);
    for (int i = 0; i < rank; i++) x[i] = M[i][n];
    return {rank, x};
}

Back to top page