Matrix Linear Algebra
(math/matrix/linear_algebra.hpp)

View this file on GitHub
Last update: 2026-07-07 14:26:59+09:00
Include: #include "math/matrix/linear_algebra.hpp"

Overview

This header provides Gaussian-elimination algorithms for Matrix<T>:

reduced row echelon form;
rank;
determinant;
inverse;
solving linear systems, including inconsistent and underdetermined systems.

These algorithms require T to behave like a field. Good choices are a floating-point type or ModInt with a prime modulus. Plain integer types are not suitable because division truncates.

Zero Tests and Precision

For non-floating types, a value is zero exactly when it equals T(). default_epsilon<T>() therefore returns T() for these types.

For float, double, and long double, values with absolute value at most eps are treated as zero, and partial pivoting chooses the largest available absolute value in each column. default_epsilon<T>() returns 1e-10, which is also the default argument. Pass an explicit tolerance when the scale or accuracy requirements of the problem differ:

auto rank = m1une::matrix::matrix_rank(matrix, 1e-12);

Floating-point answers remain approximate. For exact arithmetic modulo a prime, use ModInt.

Row Reduction and Rank

reduced_row_echelon_form(matrix, eps) returns RowReduction<T>.

Member	Meaning
`matrix`	The reduced row echelon form.
`pivot_columns`	Pivot column indices from left to right.
`rank()`	The number of pivots.

matrix_rank(matrix, eps) computes only the rank and avoids eliminating above each pivot, making it faster than constructing the full reduced form.

For an H x W matrix, both operations take $O(HW\min(H, W))$ time. The returned objects use $O(HW)$ memory because the input is copied and preserved.

Determinant and Inverse

Function	Result	Complexity
`determinant(matrix, eps)`	The determinant of a square matrix. The determinant of the `0 x 0` matrix is `1`.	$O(N^3)$ time, $O(N^2)$ memory
`inverse(matrix, eps)`	`std::optional<Matrix<T>>`; empty when the square matrix is singular.	$O(N^3)$ time, $O(N^2)$ memory

Linear Systems

solve_linear_system(coefficients, constants, eps) solves

\[A x = b.\]

The number of rows in coefficients must equal constants.size(). Rectangular coefficient matrices are supported.

The returned LinearSystemResult<T> contains:

Member / Method	Meaning
`consistent`	Whether at least one solution exists.
`particular_solution`	One solution when the system is consistent. Free variables are set to zero.
`nullspace_basis`	Basis vectors for all homogeneous solutions of `A x = 0`.
`pivot_columns`	Pivot variable indices.
`rank()`	Rank of `A`.
`nullity()`	Number of free variables when the system is consistent.
`has_unique_solution()`	Whether exactly one solution exists.

When the system is consistent, every solution can be written as the particular solution plus a linear combination of nullspace_basis. When it is inconsistent, particular_solution and nullspace_basis are empty.

The complexity is $O(HW\min(H, W))$ time and $O(HW)$ memory.

Example

#include "math/modint.hpp"
#include "math/matrix/linear_algebra.hpp"

#include <iostream>
#include <vector>

int main() {
    using mint = m1une::math::modint998244353;
    m1une::matrix::Matrix<mint> coefficients(2, 2);
    coefficients[0][0] = 2;
    coefficients[0][1] = 1;
    coefficients[1][0] = 1;
    coefficients[1][1] = 3;

    std::vector<mint> constants = {5, 7};
    auto result =
        m1une::matrix::solve_linear_system(coefficients, constants);

    if (result.has_unique_solution()) {
        std::cout << result.particular_solution[0] << ' '
                  << result.particular_solution[1] << '\n';
    }
}

Depends on

Dense Matrix (math/matrix/matrix.hpp)

Required by

Verified with

Code

#ifndef M1UNE_MATRIX_LINEAR_ALGEBRA_HPP
#define M1UNE_MATRIX_LINEAR_ALGEBRA_HPP 1

#include <optional>
#include <type_traits>
#include <vector>

#include "matrix.hpp"

namespace m1une {
namespace matrix {

template <class T>
constexpr T default_epsilon() {
    if constexpr (std::is_floating_point_v<T>) {
        return T(1e-10);
    } else {
        return T();
    }
}

namespace detail {

template <class T>
T matrix_abs(T value) {
    return value < T() ? T() - value : value;
}

template <class T>
bool is_zero(const T& value, const T& eps) {
    if constexpr (std::is_floating_point_v<T>) {
        return matrix_abs(value) <= eps;
    } else {
        (void)eps;
        return value == T();
    }
}

template <class T>
int choose_pivot(const Matrix<T>& matrix, int first_row, int col, const T& eps) {
    int pivot = -1;
    if constexpr (std::is_floating_point_v<T>) {
        for (int row = first_row; row < matrix.rows(); row++) {
            if (is_zero(matrix[row][col], eps)) continue;
            if (pivot == -1 || matrix_abs(matrix[pivot][col]) < matrix_abs(matrix[row][col])) {
                pivot = row;
            }
        }
    } else {
        for (int row = first_row; row < matrix.rows(); row++) {
            if (!is_zero(matrix[row][col], eps)) {
                pivot = row;
                break;
            }
        }
    }
    return pivot;
}

template <class T>
std::vector<int> row_reduce(Matrix<T>& matrix, int pivot_col_limit, const T& eps,
                            bool reduced) {
    std::vector<int> pivot_columns;
    int pivot_row = 0;
    for (int col = 0; col < pivot_col_limit && pivot_row < matrix.rows(); col++) {
        int pivot = choose_pivot(matrix, pivot_row, col, eps);
        if (pivot == -1) continue;
        matrix.swap_rows(pivot_row, pivot);

        const T pivot_value = matrix[pivot_row][col];
        if (reduced) {
            for (int j = col; j < matrix.cols(); j++) matrix[pivot_row][j] /= pivot_value;
        }

        const int first_row = reduced ? 0 : pivot_row + 1;
        for (int row = first_row; row < matrix.rows(); row++) {
            if (row == pivot_row || is_zero(matrix[row][col], eps)) continue;
            T factor = matrix[row][col];
            if (!reduced) factor /= pivot_value;
            matrix[row][col] = T();
            for (int j = col + 1; j < matrix.cols(); j++) {
                matrix[row][j] -= factor * matrix[pivot_row][j];
            }
        }

        pivot_columns.push_back(col);
        pivot_row++;
    }

    if constexpr (std::is_floating_point_v<T>) {
        for (T& value : matrix.data()) {
            if (is_zero(value, eps)) value = T();
        }
    }
    return pivot_columns;
}

}  // namespace detail

template <class T>
struct RowReduction {
    Matrix<T> matrix;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }
};

template <class T>
RowReduction<T> reduced_row_echelon_form(Matrix<T> matrix,
                                         T eps = default_epsilon<T>()) {
    RowReduction<T> result;
    result.pivot_columns = detail::row_reduce(matrix, matrix.cols(), eps, true);
    result.matrix = std::move(matrix);
    return result;
}

template <class T>
int matrix_rank(Matrix<T> matrix, T eps = default_epsilon<T>()) {
    return int(detail::row_reduce(matrix, matrix.cols(), eps, false).size());
}

template <class T>
T determinant(Matrix<T> matrix, T eps = default_epsilon<T>()) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    T result = T(1);
    bool negate = false;

    for (int col = 0; col < size; col++) {
        int pivot = detail::choose_pivot(matrix, col, col, eps);
        if (pivot == -1) return T();
        if (pivot != col) {
            matrix.swap_rows(pivot, col);
            negate = !negate;
        }

        const T pivot_value = matrix[col][col];
        result *= pivot_value;
        for (int row = col + 1; row < size; row++) {
            if (detail::is_zero(matrix[row][col], eps)) continue;
            const T factor = matrix[row][col] / pivot_value;
            matrix[row][col] = T();
            for (int j = col + 1; j < size; j++) {
                matrix[row][j] -= factor * matrix[col][j];
            }
        }
    }
    return negate ? T() - result : result;
}

template <class T>
std::optional<Matrix<T>> inverse(const Matrix<T>& matrix,
                                 T eps = default_epsilon<T>()) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    Matrix<T> augmented(size, size * 2);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            augmented[row][col] = matrix[row][col];
        }
        augmented[row][size + row] = T(1);
    }

    const std::vector<int> pivots = detail::row_reduce(augmented, size, eps, true);
    if (int(pivots.size()) != size) return std::nullopt;

    Matrix<T> result(size, size);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            result[row][col] = augmented[row][size + col];
        }
    }
    return result;
}

template <class T>
struct LinearSystemResult {
    bool consistent = false;
    std::vector<T> particular_solution;
    std::vector<std::vector<T>> nullspace_basis;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }

    int nullity() const {
        return consistent ? int(nullspace_basis.size()) : 0;
    }

    bool has_unique_solution() const {
        return consistent && nullspace_basis.empty();
    }
};

template <class T>
LinearSystemResult<T> solve_linear_system(const Matrix<T>& coefficients,
                                          const std::vector<T>& constants,
                                          T eps = default_epsilon<T>()) {
    assert(coefficients.rows() == int(constants.size()));
    const int equation_count = coefficients.rows();
    const int variable_count = coefficients.cols();
    Matrix<T> augmented(equation_count, variable_count + 1);
    for (int row = 0; row < equation_count; row++) {
        for (int col = 0; col < variable_count; col++) {
            augmented[row][col] = coefficients[row][col];
        }
        augmented[row][variable_count] = constants[std::size_t(row)];
    }

    LinearSystemResult<T> result;
    result.pivot_columns =
        detail::row_reduce(augmented, variable_count, eps, true);

    for (int row = result.rank(); row < equation_count; row++) {
        bool zero_left = true;
        for (int col = 0; col < variable_count; col++) {
            if (!detail::is_zero(augmented[row][col], eps)) {
                zero_left = false;
                break;
            }
        }
        if (zero_left && !detail::is_zero(augmented[row][variable_count], eps)) {
            return result;
        }
    }

    result.consistent = true;
    result.particular_solution.assign(std::size_t(variable_count), T());
    std::vector<bool> is_pivot(std::size_t(variable_count), false);
    for (int row = 0; row < result.rank(); row++) {
        const int col = result.pivot_columns[std::size_t(row)];
        is_pivot[std::size_t(col)] = true;
        result.particular_solution[std::size_t(col)] = augmented[row][variable_count];
    }

    for (int free_col = 0; free_col < variable_count; free_col++) {
        if (is_pivot[std::size_t(free_col)]) continue;
        std::vector<T> direction(static_cast<std::size_t>(variable_count));
        direction[std::size_t(free_col)] = T(1);
        for (int row = 0; row < result.rank(); row++) {
            const int pivot_col = result.pivot_columns[std::size_t(row)];
            direction[std::size_t(pivot_col)] = T() - augmented[row][free_col];
        }
        result.nullspace_basis.push_back(std::move(direction));
    }
    return result;
}

}  // namespace matrix
}  // namespace m1une

#endif  // M1UNE_MATRIX_LINEAR_ALGEBRA_HPP

#line 1 "math/matrix/linear_algebra.hpp"



#include <optional>
#include <type_traits>
#include <vector>

#line 1 "math/matrix/matrix.hpp"



#include <cassert>
#include <cstddef>
#include <cstdint>
#include <utility>
#line 9 "math/matrix/matrix.hpp"

namespace m1une {
namespace matrix {

template <class T>
class Matrix {
   private:
    int _rows;
    int _cols;
    std::vector<T> _data;

    static std::size_t storage_size(int rows, int cols) {
        assert(rows >= 0);
        assert(cols >= 0);
        return std::size_t(rows) * std::size_t(cols);
    }

   public:
    using value_type = T;

    Matrix() : _rows(0), _cols(0) {}

    Matrix(int rows, int cols, const T& value = T())
        : _rows(rows), _cols(cols), _data(storage_size(rows, cols), value) {}

    Matrix(int rows, int cols, std::vector<T> values)
        : _rows(rows), _cols(cols), _data(std::move(values)) {
        assert(rows >= 0);
        assert(cols >= 0);
        assert(_data.size() == std::size_t(rows) * std::size_t(cols));
    }

    explicit Matrix(const std::vector<std::vector<T>>& values)
        : _rows(int(values.size())), _cols(values.empty() ? 0 : int(values[0].size())),
          _data(storage_size(_rows, _cols)) {
        for (int row = 0; row < _rows; row++) {
            assert(int(values[std::size_t(row)].size()) == _cols);
            for (int col = 0; col < _cols; col++) {
                (*this)[row][col] = values[std::size_t(row)][std::size_t(col)];
            }
        }
    }

    int rows() const {
        return _rows;
    }

    int cols() const {
        return _cols;
    }

    bool empty() const {
        return _rows == 0 || _cols == 0;
    }

    std::vector<T>& data() {
        return _data;
    }

    const std::vector<T>& data() const {
        return _data;
    }

    T* operator[](int row) {
        assert(0 <= row && row < _rows);
        return _data.data() + std::size_t(row) * std::size_t(_cols);
    }

    const T* operator[](int row) const {
        assert(0 <= row && row < _rows);
        return _data.data() + std::size_t(row) * std::size_t(_cols);
    }

    T& operator()(int row, int col) {
        assert(0 <= col && col < _cols);
        return (*this)[row][col];
    }

    const T& operator()(int row, int col) const {
        assert(0 <= col && col < _cols);
        return (*this)[row][col];
    }

    static Matrix identity(int size) {
        assert(size >= 0);
        Matrix result(size, size);
        for (int i = 0; i < size; i++) result[i][i] = T(1);
        return result;
    }

    Matrix transposed() const {
        Matrix result(_cols, _rows);
        for (int row = 0; row < _rows; row++) {
            for (int col = 0; col < _cols; col++) {
                result[col][row] = (*this)[row][col];
            }
        }
        return result;
    }

    void swap_rows(int first, int second) {
        assert(0 <= first && first < _rows);
        assert(0 <= second && second < _rows);
        if (first == second) return;
        for (int col = 0; col < _cols; col++) {
            std::swap((*this)[first][col], (*this)[second][col]);
        }
    }

    Matrix& operator+=(const Matrix& rhs) {
        assert(_rows == rhs._rows && _cols == rhs._cols);
        for (std::size_t i = 0; i < _data.size(); i++) _data[i] += rhs._data[i];
        return *this;
    }

    Matrix& operator-=(const Matrix& rhs) {
        assert(_rows == rhs._rows && _cols == rhs._cols);
        for (std::size_t i = 0; i < _data.size(); i++) _data[i] -= rhs._data[i];
        return *this;
    }

    Matrix& operator*=(const T& scalar) {
        for (T& value : _data) value *= scalar;
        return *this;
    }

    Matrix& operator/=(const T& scalar) {
        for (T& value : _data) value /= scalar;
        return *this;
    }

    Matrix& operator*=(const Matrix& rhs) {
        return *this = *this * rhs;
    }

    Matrix operator+() const {
        return *this;
    }

    Matrix operator-() const {
        Matrix result = *this;
        for (T& value : result._data) value = T() - value;
        return result;
    }

    friend Matrix operator+(Matrix lhs, const Matrix& rhs) {
        return lhs += rhs;
    }

    friend Matrix operator-(Matrix lhs, const Matrix& rhs) {
        return lhs -= rhs;
    }

    friend Matrix operator*(Matrix lhs, const T& rhs) {
        return lhs *= rhs;
    }

    friend Matrix operator*(const T& lhs, Matrix rhs) {
        return rhs *= lhs;
    }

    friend Matrix operator/(Matrix lhs, const T& rhs) {
        return lhs /= rhs;
    }

    friend Matrix operator*(const Matrix& lhs, const Matrix& rhs) {
        assert(lhs._cols == rhs._rows);
        Matrix result(lhs._rows, rhs._cols);
        for (int row = 0; row < lhs._rows; row++) {
            T* output = result[row];
            for (int middle = 0; middle < lhs._cols; middle++) {
                const T coefficient = lhs[row][middle];
                if (coefficient == T()) continue;
                const T* input = rhs[middle];
                for (int col = 0; col < rhs._cols; col++) {
                    output[col] += coefficient * input[col];
                }
            }
        }
        return result;
    }

    friend std::vector<T> operator*(const Matrix& lhs, const std::vector<T>& rhs) {
        assert(lhs._cols == int(rhs.size()));
        std::vector<T> result(std::size_t(lhs._rows));
        for (int row = 0; row < lhs._rows; row++) {
            T value = T();
            for (int col = 0; col < lhs._cols; col++) {
                value += lhs[row][col] * rhs[std::size_t(col)];
            }
            result[std::size_t(row)] = value;
        }
        return result;
    }

    friend std::vector<T> operator*(const std::vector<T>& lhs, const Matrix& rhs) {
        assert(int(lhs.size()) == rhs._rows);
        std::vector<T> result(std::size_t(rhs._cols));
        for (int row = 0; row < rhs._rows; row++) {
            if (lhs[std::size_t(row)] == T()) continue;
            for (int col = 0; col < rhs._cols; col++) {
                result[std::size_t(col)] += lhs[std::size_t(row)] * rhs[row][col];
            }
        }
        return result;
    }

    bool operator==(const Matrix& rhs) const {
        return _rows == rhs._rows && _cols == rhs._cols && _data == rhs._data;
    }

    bool operator!=(const Matrix& rhs) const {
        return !(*this == rhs);
    }

    Matrix pow(std::uint64_t exponent) const {
        assert(_rows == _cols);
        Matrix result = identity(_rows);
        Matrix base = *this;
        while (exponent > 0) {
            if (exponent & 1) result *= base;
            exponent >>= 1;
            if (exponent > 0) base *= base;
        }
        return result;
    }
};

}  // namespace matrix
}  // namespace m1une


#line 9 "math/matrix/linear_algebra.hpp"

namespace m1une {
namespace matrix {

template <class T>
constexpr T default_epsilon() {
    if constexpr (std::is_floating_point_v<T>) {
        return T(1e-10);
    } else {
        return T();
    }
}

namespace detail {

template <class T>
T matrix_abs(T value) {
    return value < T() ? T() - value : value;
}

template <class T>
bool is_zero(const T& value, const T& eps) {
    if constexpr (std::is_floating_point_v<T>) {
        return matrix_abs(value) <= eps;
    } else {
        (void)eps;
        return value == T();
    }
}

template <class T>
int choose_pivot(const Matrix<T>& matrix, int first_row, int col, const T& eps) {
    int pivot = -1;
    if constexpr (std::is_floating_point_v<T>) {
        for (int row = first_row; row < matrix.rows(); row++) {
            if (is_zero(matrix[row][col], eps)) continue;
            if (pivot == -1 || matrix_abs(matrix[pivot][col]) < matrix_abs(matrix[row][col])) {
                pivot = row;
            }
        }
    } else {
        for (int row = first_row; row < matrix.rows(); row++) {
            if (!is_zero(matrix[row][col], eps)) {
                pivot = row;
                break;
            }
        }
    }
    return pivot;
}

template <class T>
std::vector<int> row_reduce(Matrix<T>& matrix, int pivot_col_limit, const T& eps,
                            bool reduced) {
    std::vector<int> pivot_columns;
    int pivot_row = 0;
    for (int col = 0; col < pivot_col_limit && pivot_row < matrix.rows(); col++) {
        int pivot = choose_pivot(matrix, pivot_row, col, eps);
        if (pivot == -1) continue;
        matrix.swap_rows(pivot_row, pivot);

        const T pivot_value = matrix[pivot_row][col];
        if (reduced) {
            for (int j = col; j < matrix.cols(); j++) matrix[pivot_row][j] /= pivot_value;
        }

        const int first_row = reduced ? 0 : pivot_row + 1;
        for (int row = first_row; row < matrix.rows(); row++) {
            if (row == pivot_row || is_zero(matrix[row][col], eps)) continue;
            T factor = matrix[row][col];
            if (!reduced) factor /= pivot_value;
            matrix[row][col] = T();
            for (int j = col + 1; j < matrix.cols(); j++) {
                matrix[row][j] -= factor * matrix[pivot_row][j];
            }
        }

        pivot_columns.push_back(col);
        pivot_row++;
    }

    if constexpr (std::is_floating_point_v<T>) {
        for (T& value : matrix.data()) {
            if (is_zero(value, eps)) value = T();
        }
    }
    return pivot_columns;
}

}  // namespace detail

template <class T>
struct RowReduction {
    Matrix<T> matrix;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }
};

template <class T>
RowReduction<T> reduced_row_echelon_form(Matrix<T> matrix,
                                         T eps = default_epsilon<T>()) {
    RowReduction<T> result;
    result.pivot_columns = detail::row_reduce(matrix, matrix.cols(), eps, true);
    result.matrix = std::move(matrix);
    return result;
}

template <class T>
int matrix_rank(Matrix<T> matrix, T eps = default_epsilon<T>()) {
    return int(detail::row_reduce(matrix, matrix.cols(), eps, false).size());
}

template <class T>
T determinant(Matrix<T> matrix, T eps = default_epsilon<T>()) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    T result = T(1);
    bool negate = false;

    for (int col = 0; col < size; col++) {
        int pivot = detail::choose_pivot(matrix, col, col, eps);
        if (pivot == -1) return T();
        if (pivot != col) {
            matrix.swap_rows(pivot, col);
            negate = !negate;
        }

        const T pivot_value = matrix[col][col];
        result *= pivot_value;
        for (int row = col + 1; row < size; row++) {
            if (detail::is_zero(matrix[row][col], eps)) continue;
            const T factor = matrix[row][col] / pivot_value;
            matrix[row][col] = T();
            for (int j = col + 1; j < size; j++) {
                matrix[row][j] -= factor * matrix[col][j];
            }
        }
    }
    return negate ? T() - result : result;
}

template <class T>
std::optional<Matrix<T>> inverse(const Matrix<T>& matrix,
                                 T eps = default_epsilon<T>()) {
    assert(matrix.rows() == matrix.cols());
    const int size = matrix.rows();
    Matrix<T> augmented(size, size * 2);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            augmented[row][col] = matrix[row][col];
        }
        augmented[row][size + row] = T(1);
    }

    const std::vector<int> pivots = detail::row_reduce(augmented, size, eps, true);
    if (int(pivots.size()) != size) return std::nullopt;

    Matrix<T> result(size, size);
    for (int row = 0; row < size; row++) {
        for (int col = 0; col < size; col++) {
            result[row][col] = augmented[row][size + col];
        }
    }
    return result;
}

template <class T>
struct LinearSystemResult {
    bool consistent = false;
    std::vector<T> particular_solution;
    std::vector<std::vector<T>> nullspace_basis;
    std::vector<int> pivot_columns;

    int rank() const {
        return int(pivot_columns.size());
    }

    int nullity() const {
        return consistent ? int(nullspace_basis.size()) : 0;
    }

    bool has_unique_solution() const {
        return consistent && nullspace_basis.empty();
    }
};

template <class T>
LinearSystemResult<T> solve_linear_system(const Matrix<T>& coefficients,
                                          const std::vector<T>& constants,
                                          T eps = default_epsilon<T>()) {
    assert(coefficients.rows() == int(constants.size()));
    const int equation_count = coefficients.rows();
    const int variable_count = coefficients.cols();
    Matrix<T> augmented(equation_count, variable_count + 1);
    for (int row = 0; row < equation_count; row++) {
        for (int col = 0; col < variable_count; col++) {
            augmented[row][col] = coefficients[row][col];
        }
        augmented[row][variable_count] = constants[std::size_t(row)];
    }

    LinearSystemResult<T> result;
    result.pivot_columns =
        detail::row_reduce(augmented, variable_count, eps, true);

    for (int row = result.rank(); row < equation_count; row++) {
        bool zero_left = true;
        for (int col = 0; col < variable_count; col++) {
            if (!detail::is_zero(augmented[row][col], eps)) {
                zero_left = false;
                break;
            }
        }
        if (zero_left && !detail::is_zero(augmented[row][variable_count], eps)) {
            return result;
        }
    }

    result.consistent = true;
    result.particular_solution.assign(std::size_t(variable_count), T());
    std::vector<bool> is_pivot(std::size_t(variable_count), false);
    for (int row = 0; row < result.rank(); row++) {
        const int col = result.pivot_columns[std::size_t(row)];
        is_pivot[std::size_t(col)] = true;
        result.particular_solution[std::size_t(col)] = augmented[row][variable_count];
    }

    for (int free_col = 0; free_col < variable_count; free_col++) {
        if (is_pivot[std::size_t(free_col)]) continue;
        std::vector<T> direction(static_cast<std::size_t>(variable_count));
        direction[std::size_t(free_col)] = T(1);
        for (int row = 0; row < result.rank(); row++) {
            const int pivot_col = result.pivot_columns[std::size_t(row)];
            direction[std::size_t(pivot_col)] = T() - augmented[row][free_col];
        }
        result.nullspace_basis.push_back(std::move(direction));
    }
    return result;
}

}  // namespace matrix
}  // namespace m1une

m1une's library

This documentation is automatically generated by online-judge-tools/verification-helper

Matrix Linear Algebra
(math/matrix/linear_algebra.hpp)

Overview

Zero Tests and Precision

Row Reduction and Rank

Determinant and Inverse

Linear Systems

Example

Depends on

Required by

Verified with

Code

Matrix Linear Algebra (math/matrix/linear_algebra.hpp)

Overview

Zero Tests and Precision

Row Reduction and Rank

Determinant and Inverse

Linear Systems

Example

Depends on

Required by

Verified with

Code

Matrix Linear Algebra
(math/matrix/linear_algebra.hpp)