linear_equation
(src/matrix/linear_equation.hpp)

• View this file on GitHub
• Last update: 2024-11-09 02:30:41+09:00
• Include: #include "src/matrix/linear_equation.hpp"

linear_equation

vector<vector<T>> linear_equation(Matrix<T> A, Matrix<T> b)

線形方程式 $Ax = b$ を解きます．

返り値は ${v, w_1, w_2, …, w_k}$ という vector<vector<T>> です．
$v$ は方程式の解のうちの $1$ つである vector<T> です．
$w_i$ は解空間の（それぞれ独立な）基底ベクトルである vector<T> です．

制約

• $A.H() = b.H()$
• $b.W() = 1$

計算量

行列 $A$ のサイズを $H \times W$ として

• $O(HW^2)$

Depends on

• gauss_elimination (src/matrix/gauss_elimination.hpp)
• Matrix (src/matrix/matrix.hpp)
• template (src/template/template.hpp)

Verified with

• verify/library_checker/linear_algebra/system_of_linear_equations.test.cpp

Code

#pragma once
#include "../template/template.hpp"
#include "./matrix.hpp"
#include "./gauss_elimination.hpp"
template <typename T>
vector<vector<T>> linear_equation(const Matrix<T>& a, const Matrix<T>& b) {
    assert(a.H() == b.H() and b.W() == 1);
    const int h = a.H(), w = a.W();
    Matrix<T> A(h, w + 1);
    for(int i = 0; i < h; ++i) {
        for(int j = 0; j < w; ++j) {
            A[i][j] = a[i][j];
        }
        A[i][w] = b[i][0];
    }
    const int rank = gauss_elimination(A, w).first;
    for(int i = rank; i < h; ++i) {
        if(A[i][w] != 0) return vector<vector<T>>{};
    }
    vector<vector<T>> res(1, vector<T>(w));
    vector<int> pivot(w, -1);
    for(int i = 0, j = 0; i < rank; ++i) {
        while(A[i][j] == 0) ++j;
        res[0][j] = A[i][w], pivot[j] = i;
    }
    for(int j = 0; j < w; ++j) {
        if(pivot[j] == -1) {
            vector<T> x(w);
            x[j] = 1;
            for(int k = 0; k < j; ++k) {
                if(pivot[k] != -1) x[k] = -A[pivot[k]][j];
            }
            res.emplace_back(x);
        }
    }
    return res;
}

#line 2 "src/template/template.hpp"
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
using P = pair<long long, long long>;
#define rep(i, a, b) for(long long i = (a); i < (b); ++i)
#define rrep(i, a, b) for(long long i = (a); i >= (b); --i)
constexpr long long inf = 4e18;
struct SetupIO {
    SetupIO() {
        ios::sync_with_stdio(0);
        cin.tie(0);
        cout << fixed << setprecision(30);
    }
} setup_io;
#line 3 "src/matrix/matrix.hpp"
template <typename T>
struct Matrix {
    Matrix(const int h, const int w, const T& val = 0)
        : h(h), w(w), A(h, vector<T>(w, val)) {}
    int H() const {
        return h;
    }
    int W() const {
        return w;
    }
    const vector<T>& operator[](const int i) const {
        assert(0 <= i and i < h);
        return A[i];
    }
    vector<T>& operator[](const int i) {
        assert(0 <= i and i < h);
        return A[i];
    }
    static Matrix I(const int n) {
        Matrix mat(n, n);
        for(int i = 0; i < n; ++i) mat[i][i] = 1;
        return mat;
    }
    Matrix& operator+=(const Matrix& B) {
        assert(h == B.h and w == B.w);
        for(int i = 0; i < h; ++i) {
            for(int j = 0; j < w; ++j) {
                (*this)[i][j] += B[i][j];
            }
        }
        return (*this);
    }
    Matrix& operator-=(const Matrix& B) {
        assert(h == B.h and w == B.w);
        for(int i = 0; i < h; ++i) {
            for(int j = 0; j < w; ++j) {
                (*this)[i][j] -= B[i][j];
            }
        }
        return (*this);
    }
    Matrix& operator*=(const Matrix& B) {
        assert(w == B.h);
        vector<vector<T>> C(h, vector<T>(B.w, 0));
        for(int i = 0; i < h; ++i) {
            for(int k = 0; k < w; ++k) {
                for(int j = 0; j < B.w; ++j) {
                    C[i][j] += (*this)[i][k] * B[k][j];
                }
            }
        }
        A.swap(C);
        return (*this);
    }
    Matrix& pow(long long t) {
        assert(h == w);
        assert(t >= 0);
        Matrix B = Matrix::I(h);
        while(t > 0) {
            if(t & 1ll) B *= (*this);
            (*this) *= (*this);
            t >>= 1ll;
        }
        A.swap(B.A);
        return (*this);
    }
    Matrix operator+(const Matrix& B) const {
        return (Matrix(*this) += B);
    }
    Matrix operator-(const Matrix& B) const {
        return (Matrix(*this) -= B);
    }
    Matrix operator*(const Matrix& B) const {
        return (Matrix(*this) *= B);
    }
    bool operator==(const Matrix& B) const {
        assert(h == B.H() and w == B.W());
        for(int i = 0; i < h; ++i) {
            for(int j = 0; j < w; ++j) {
                if(A[i][j] != B[i][j]) return false;
            }
        }
        return true;
    }
    bool operator!=(const Matrix& B) const {
        assert(h == B.H() and w == B.W());
        for(int i = 0; i < h; ++i) {
            for(int j = 0; j < w; ++j) {
                if(A[i][j] != B[i][j]) return true;
            }
        }
        return false;
    }

   private:
    int h, w;
    vector<vector<T>> A;
};
#line 4 "src/matrix/gauss_elimination.hpp"
template <typename T>
pair<int, T> gauss_elimination(Matrix<T>& a, int pivot_end = -1) {
    const int h = a.H(), w = a.W();
    int rank = 0;
    assert(-1 <= pivot_end and pivot_end <= w);
    if(pivot_end == -1) pivot_end = w;
    T det = 1;
    for(int j = 0; j < pivot_end; ++j) {
        int idx = -1;
        for(int i = rank; i < h; ++i) {
            if(a[i][j] != T(0)) {
                idx = i;
                break;
            }
        }
        if(idx == -1) {
            det = 0;
            continue;
        }
        if(rank != idx) det = -det, swap(a[rank], a[idx]);
        det *= a[rank][j];
        if(a[rank][j] != T(1)) {
            const T coeff = T(1) / a[rank][j];
            for(int k = j; k < w; ++k) a[rank][k] *= coeff;
        }
        for(int i = 0; i < h; ++i) {
            if(i == rank) continue;
            if(a[i][j] != T(0)) {
                const T coeff = a[i][j] / a[rank][j];
                for(int k = j; k < w; ++k) a[i][k] -= a[rank][k] * coeff;
            }
        }
        ++rank;
    }
    return {rank, det};
}
#line 5 "src/matrix/linear_equation.hpp"
template <typename T>
vector<vector<T>> linear_equation(const Matrix<T>& a, const Matrix<T>& b) {
    assert(a.H() == b.H() and b.W() == 1);
    const int h = a.H(), w = a.W();
    Matrix<T> A(h, w + 1);
    for(int i = 0; i < h; ++i) {
        for(int j = 0; j < w; ++j) {
            A[i][j] = a[i][j];
        }
        A[i][w] = b[i][0];
    }
    const int rank = gauss_elimination(A, w).first;
    for(int i = rank; i < h; ++i) {
        if(A[i][w] != 0) return vector<vector<T>>{};
    }
    vector<vector<T>> res(1, vector<T>(w));
    vector<int> pivot(w, -1);
    for(int i = 0, j = 0; i < rank; ++i) {
        while(A[i][j] == 0) ++j;
        res[0][j] = A[i][w], pivot[j] = i;
    }
    for(int j = 0; j < w; ++j) {
        if(pivot[j] == -1) {
            vector<T> x(w);
            x[j] = 1;
            for(int k = 0; k < j; ++k) {
                if(pivot[k] != -1) x[k] = -A[pivot[k]][j];
            }
            res.emplace_back(x);
        }
    }
    return res;
}

Back to top page