cp-library-cpp

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

View the Project on GitHub suisen-cp/cp-library-cpp

:heavy_check_mark: test/src/linear_algebra/gaussian_elimination/system_of_linear_equations.test.cpp

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/system_of_linear_equations"

#include <iostream>
#include <atcoder/modint>

#include "library/linear_algebra/gaussian_elimination.hpp"

using mint = atcoder::modint998244353;

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, m;
    std::cin >> n >> m;
    std::vector A(n, std::vector(m, mint(0)));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            int val;
            std::cin >> val;
            A[i][j] = val;
        }
    }
    std::vector<mint> b(n);
    for (int i = 0; i < n; ++i) {
        int val;
        std::cin >> val;
        b[i] = val;
    }
    suisen::GaussianEliminationArithmetic<mint> solution(A, b);
    if (solution.has_solution()) {
        int r = solution.dimension();
        const auto c = *solution.get_solution();
        const auto &basis = solution.get_basis();
        std::cout << r << '\n';
        for (int i = 0; i < m; ++i) {
            std::cout << c[i].val() << " \n"[i == m - 1];
        }
        for (const auto &x : basis) {
            for (int i = 0; i < m; ++i) {
                std::cout << x[i].val() << " \n"[i == m - 1];
            }
        }
    } else {
        std::cout << -1 << '\n';
    }
    return 0;
}
#line 1 "test/src/linear_algebra/gaussian_elimination/system_of_linear_equations.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/system_of_linear_equations"

#include <iostream>
#include <atcoder/modint>

#line 1 "library/linear_algebra/gaussian_elimination.hpp"



#include <algorithm>
#include <cmath>
#include <optional>
#include <vector>

namespace suisen {
    namespace internal::gauss_jordan {
        template <typename T, typename EqZero, std::enable_if_t<std::negation_v<std::is_floating_point<T>>, std::nullptr_t> = nullptr>
        std::pair<size_t, size_t> pivoting(const std::vector<std::vector<T>>& Ab, const size_t i, EqZero equals_zero) {
            const size_t n = Ab.size(), m = Ab[0].size() - 1;
            size_t mse = m, pivot = n;
            for (size_t row = i; row < n; ++row) {
                const auto &v = Ab[row];
                size_t col = std::find_if_not(v.begin(), v.begin() + mse, equals_zero) - v.begin();
                if (col < mse) mse = col, pivot = row;
            }
            return { mse, pivot };
        }
        // Gauss pivoting
        template <typename T, typename EqZero, std::enable_if_t<std::is_floating_point_v<T>, std::nullptr_t> = nullptr>
        std::pair<size_t, size_t> pivoting(const std::vector<std::vector<T>>& Ab, const size_t i, EqZero equals_zero) {
            const size_t n = Ab.size(), m = Ab[0].size() - 1;
            size_t mse = m, pivot = n;
            T max_val = 0;
            for (size_t row = i; row < n; ++row) {
                const auto &v = Ab[row];
                if (mse < m and std::abs(v[mse]) > max_val) pivot = row, max_val = std::abs(v[mse]);
                size_t col = std::find_if_not(v.begin(), v.begin() + mse, equals_zero) - v.begin();
                if (col < mse) mse = col, pivot = row, max_val = std::abs(Ab[row][col]);
            }
            return { mse, pivot };
        }

        template <typename T> constexpr T add_fp_f2(T x, T y) { return x ^ y; }
        template <typename T> constexpr T add_inv_fp_f2(T x) { return x; }
        template <typename T> constexpr T mul_fp_f2(T x, T y) { return x & y; }
        template <typename T> constexpr T mul_inv_fp_f2(T x) { return x; }

        template <typename T> constexpr T add_fp_arithmetic(T x, T y) { return x + y; }
        template <typename T> constexpr T add_inv_fp_arithmetic(T x) { return 0 - x; }
        template <typename T> constexpr T mul_fp_arithmetic(T x, T y) { return x * y; }
        template <typename T> constexpr T mul_inv_fp_arithmetic(T x) { return 1 / x; }
    }

    template <typename T, T(*add_fp)(T, T), T(*add_inv_fp)(T), T(*mul_fp)(T, T), T(*mul_inv_fp)(T)>
    struct GaussianElimination {
        GaussianElimination(std::vector<std::vector<T>> A, const std::vector<T>& b, const T &zero = T(0), const T &one = T(1)) {
            size_t n = A.size();
            for (size_t i = 0; i < n; ++i) A[i].push_back(b[i]);
            solve(zero, one, A);
        }
        bool has_solution() const { return not _empty; }
        bool has_unique_solution() const { return not _empty and _basis.size() == 0; }
        bool has_multiple_solutions() const { return _basis.size() > 0; }
        const std::optional<std::vector<T>> get_solution() const { return _empty ? std::nullopt : std::make_optional(_x0); }
        const std::vector<std::vector<T>>& get_basis() const { return _basis; }
        int dimension() const { return _empty ? -1 : _basis.size(); }
    private:
        std::vector<T> _x0;
        std::vector<std::vector<T>> _basis;
        bool _empty = false;

        void solve(const T &zero, const T &one, std::vector<std::vector<T>>& Ab) {
            using namespace internal::gauss_jordan;
            auto equals_zero = [&](const T& v) {
                if constexpr (std::is_floating_point_v<T>) return std::abs(v) < 1e-9;
                else return v == zero;
            };
            const size_t n = Ab.size(), m = Ab[0].size() - 1;
            for (size_t i = 0; i < n; ++i) {
                auto [mse, pivot] = pivoting(Ab, i, equals_zero);
                if (pivot == n) break;
                Ab[i].swap(Ab[pivot]);
                T mse_val_inv = mul_inv_fp(Ab[i][mse]);
                for (size_t row = i + 1; row < n; ++row) if (not equals_zero(Ab[row][mse])) {
                    T coef = add_inv_fp(mul_fp(Ab[row][mse], mse_val_inv));
                    for (size_t col = mse; col <= m; ++col) Ab[row][col] = add_fp(Ab[row][col], mul_fp(coef, Ab[i][col]));
                }
            }
            size_t basis_num = m;
            std::vector<char> down(m, false);
            _x0.assign(m, zero);
            for (size_t i = n; i-- > 0;) {
                size_t mse = m + 1;
                for (size_t col = 0; col <= m; ++col) if (not equals_zero(Ab[i][col])) {
                    mse = col;
                    break;
                }
                if (mse < m) {
                    T mse_val_inv = mul_inv_fp(Ab[i][mse]);
                    for (size_t row = 0; row < i; ++row) if (not equals_zero(Ab[row][mse])) {
                        T coef = add_inv_fp(mul_fp(Ab[row][mse], mse_val_inv));
                        for (size_t col = mse; col <= m; ++col) Ab[row][col] = add_fp(Ab[row][col], mul_fp(coef, Ab[i][col]));
                    }
                    for (size_t col = mse; col <= m; ++col) Ab[i][col] = mul_fp(Ab[i][col], mse_val_inv);
                    _x0[mse] = Ab[i][m];
                    down[mse] = true;
                    --basis_num;
                } else if (mse == m) {
                    _empty = true;
                    return;
                }
            }
            _basis.assign(basis_num, std::vector<T>(m));
            int basis_id = 0;
            for (size_t j = 0; j < m; ++j) if (not down[j]) {
                for (size_t j2 = 0, i = 0; j2 < m; ++j2) _basis[basis_id][j2] = down[j2] ? Ab[i++][j] : zero;
                _basis[basis_id++][j] = add_inv_fp(one);
            }
        }
    };

    template <typename T>
    using GaussianEliminationF2 = GaussianElimination<
        T,
        internal::gauss_jordan::add_fp_f2, internal::gauss_jordan::add_inv_fp_f2,
        internal::gauss_jordan::mul_fp_f2, internal::gauss_jordan::mul_inv_fp_f2>;
    template <typename T>
    using GaussianEliminationArithmetic = GaussianElimination<
        T,
        internal::gauss_jordan::add_fp_arithmetic, internal::gauss_jordan::add_inv_fp_arithmetic,
        internal::gauss_jordan::mul_fp_arithmetic, internal::gauss_jordan::mul_inv_fp_arithmetic>;
} // namespace suisen


#line 7 "test/src/linear_algebra/gaussian_elimination/system_of_linear_equations.test.cpp"

using mint = atcoder::modint998244353;

int main() {
    std::ios::sync_with_stdio(false);
    std::cin.tie(nullptr);

    int n, m;
    std::cin >> n >> m;
    std::vector A(n, std::vector(m, mint(0)));
    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < m; ++j) {
            int val;
            std::cin >> val;
            A[i][j] = val;
        }
    }
    std::vector<mint> b(n);
    for (int i = 0; i < n; ++i) {
        int val;
        std::cin >> val;
        b[i] = val;
    }
    suisen::GaussianEliminationArithmetic<mint> solution(A, b);
    if (solution.has_solution()) {
        int r = solution.dimension();
        const auto c = *solution.get_solution();
        const auto &basis = solution.get_basis();
        std::cout << r << '\n';
        for (int i = 0; i < m; ++i) {
            std::cout << c[i].val() << " \n"[i == m - 1];
        }
        for (const auto &x : basis) {
            for (int i = 0; i < m; ++i) {
                std::cout << x[i].val() << " \n"[i == m - 1];
            }
        }
    } else {
        std::cout << -1 << '\n';
    }
    return 0;
}
Back to top page