test/library-checker/Matrix/Inverse.test.cpp

View this file on GitHub
Last update: 2025-11-24 18:49:07+09:00
Problem: https://judge.yosupo.jp/problem/inverse_matrix

Depends on

Code

#define PROBLEM "https://judge.yosupo.jp/problem/inverse_matrix"
#include <bits/stdc++.h>

#include "library/linearalgebra/Matrix.hpp"
#include "library/mod/Modint.hpp"

using mint = Mint<long long, 998244353>;
using M = Matrix<mint>;

#define REP(i, n) for (int i = 0; i < (n); i++)

int main() {
    int n;
    std::cin >> n;
    M A(n, n);
    std::cin >> A;
    auto B = A.inv();
    if (B.has_value())
        REP (i, n)
            REP (j, n)
                std::cout << B.value()[i][j] << "\n "[j + 1 < n];
    else
        std::cout << -1 << "\n";
}

#line 1 "test/library-checker/Matrix/Inverse.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/inverse_matrix"
#include <bits/stdc++.h>

#line 2 "library/linearalgebra/Matrix.hpp"
#define REP_(i, n) for (int i = 0; i < (n); i++)
#define REP2_(i, s, n) for (int i = (s); i < (n); i++)
template <typename K> struct Matrix {
    using value_type = K;
    using vec = std::vector<K>;
    using mat = std::vector<vec>;
    size_t r, c;
    mat M;

    Matrix(size_t r, size_t c) : r(r), c(c), M(r, vec(c, K())) {}
    Matrix(mat A) : M(A), r(A.size()), c(A[0].size()) {}

    vec &operator[](size_t k) { return M[k]; }
    const vec &operator[](size_t k) const { return M[k]; }

    Matrix &operator+=(const Matrix &A) {
        assert(r == A.r && c == A.c);
        REP_(i, r) REP_(j, c) M[i][j] += A[i][j];
        return *this;
    }
    Matrix &operator-=(const Matrix &A) {
        assert(r == A.r && c == A.c);
        REP_(i, r) REP_(j, c) M[i][j] -= A[i][j];
        return *this;
    }
    Matrix operator+(const Matrix &A) { return Matrix(M) += A; }
    Matrix operator-(const Matrix &A) { return Matrix(M) -= A; }

    friend Matrix operator*(const Matrix &A, const Matrix &B) {
        assert(A.c == B.r);
        Matrix res(A.r, B.c);
        REP_(i, A.r) REP_(k, A.c) REP_(j, B.c) res[i][j] += A[i][k] * B[k][j];
        return res;
    }
    Matrix &operator*=(const Matrix &A) {
        M = ((*this) * A).M;
        return *this;
    }

    bool operator==(const Matrix &A) {
        if (r != A.r || c != A.c)
            return false;
        REP_(i, r) REP_(j, c) if (M[i][j] != A[i][j]) return false;
        return true;
    }
    bool operator!=(const Matrix &A) { return !((*this) == A); }

    static Matrix I(size_t n) {
        Matrix res(n, n);
        REP_(i, n) res[i][i] = K(1);
        return res;
    }

    Matrix pow(long long n) const {
        assert(n >= 0 && r == c);
        Matrix A(M), res = I(r);
        while (n) {
            if (n & 1)
                res *= A;
            A *= A;
            n >>= 1;
        }
        return res;
    }

    std::pair<int, int> GaussJordan() {
        int rnk = 0, cnt = 0;
        REP_(k, c) {
            if (M[rnk][k] == 0)
                REP2_(i, rnk + 1, r)
            if (M[i][k] != 0) {
                std::swap(M[i], M[rnk]);
                cnt ^= 1;
                break;
            }
            if (M[rnk][k] == 0)
                continue;
            REP_(i, r) if (i != rnk) {
                K x = M[i][k] / M[rnk][k];
                REP_(j, c) M[i][j] -= M[rnk][j] * x;
            }
            if (++rnk == r)
                break;
        }
        return {rnk, cnt};
    }

    K det() const {
        assert(r == c);
        Matrix A(M);
        const auto &[rnk, cnt] = A.GaussJordan();
        if (rnk != r)
            return 0;
        K res = 1;
        REP_(i, r) res *= A[i][i];
        return (cnt ? -res : res);
    }

    std::optional<Matrix> inv() const {
        assert(r == c);
        Matrix A(r, c + c);
        REP_(i, r) REP_(j, c) A[i][j] = M[i][j];
        REP_(i, r) REP_(j, c) A[i][c + j] = K(i == j);
        A.GaussJordan();
        REP_(i, r) if (A[i][i] == 0) return std::nullopt;
        Matrix res(r, c);
        REP_(i, r) REP_(j, c) res[i][j] = A[i][c + j] / A[i][i];
        return res;
    }

    friend std::ostream &operator<<(std::ostream &os, const Matrix &M) {
        os << M.M;
        return os;
    }
    friend std::istream &operator>>(std::istream &is, Matrix &M) {
        REP_(i, M.r) REP_(j, M.c) is >> M.M[i][j];
        return is;
    }
};
#undef REP_
#undef REP2_
#line 2 "library/math/ExtraGCD.hpp"
using ll = long long;
std::pair<ll, ll> ext_gcd(ll a, ll b) {
    if (b == 0)
        return {1, 0};
    auto [X, Y] = ext_gcd(b, a % b);
    // bX + (a%b)Y = gcd(a,b)
    // a%b = a - b(a/b)
    // ∴ aY + b(X-(a/b)Y) = gcd(a,b)
    ll x = Y, y = X - (a / b) * Y;
    return {x, y};
}
#line 3 "library/mod/Modint.hpp"
template <typename T, T MOD = 998244353> struct Mint {
    inline static constexpr T mod = MOD;
    T v;
    Mint() : v(0) {}
    Mint(signed v) : v(v) {}
    Mint(long long t) {
        v = t % MOD;
        if (v < 0)
            v += MOD;
    }

    static Mint raw(int v) {
        Mint x;
        x.v = v;
        return x;
    }

    Mint pow(long long k) const {
        Mint res(1), tmp(v);
        while (k) {
            if (k & 1)
                res *= tmp;
            tmp *= tmp;
            k >>= 1;
        }
        return res;
    }

    static Mint add_identity() { return Mint(0); }
    static Mint mul_identity() { return Mint(1); }

    // Mint inv()const{return pow(MOD-2);}
    Mint inv() const { return Mint(ext_gcd(v, mod).first); }

    Mint &operator+=(Mint a) {
        v += a.v;
        if (v >= MOD)
            v -= MOD;
        return *this;
    }
    Mint &operator-=(Mint a) {
        v += MOD - a.v;
        if (v >= MOD)
            v -= MOD;
        return *this;
    }
    Mint &operator*=(Mint a) {
        v = 1LL * v * a.v % MOD;
        return *this;
    }
    Mint &operator/=(Mint a) { return (*this) *= a.inv(); }

    Mint operator+(Mint a) const { return Mint(v) += a; }
    Mint operator-(Mint a) const { return Mint(v) -= a; }
    Mint operator*(Mint a) const { return Mint(v) *= a; }
    Mint operator/(Mint a) const { return Mint(v) /= a; }
#define FRIEND(op)                                                             \
    friend Mint operator op(int a, Mint b) { return Mint(a) op b; }
    FRIEND(+);
    FRIEND(-);
    FRIEND(*);
    FRIEND(/);
#undef FRIEND
    Mint operator+() const { return *this; }
    Mint operator-() const { return v ? Mint(MOD - v) : Mint(v); }

    bool operator==(const Mint a) const { return v == a.v; }
    bool operator!=(const Mint a) const { return v != a.v; }

    static Mint comb(long long n, int k) {
        Mint num(1), dom(1);
        for (int i = 0; i < k; i++) {
            num *= Mint(n - i);
            dom *= Mint(i + 1);
        }
        return num / dom;
    }

    friend std::ostream &operator<<(std::ostream &os, const Mint &m) {
        os << m.v;
        return os;
    }
    friend std::istream &operator>>(std::istream &is, Mint &m) {
        is >> m.v;
        m.v %= MOD;
        if (m.v < 0)
            m.v += MOD;
        return is;
    }
};
#line 6 "test/library-checker/Matrix/Inverse.test.cpp"

using mint = Mint<long long, 998244353>;
using M = Matrix<mint>;

#define REP(i, n) for (int i = 0; i < (n); i++)

int main() {
    int n;
    std::cin >> n;
    M A(n, n);
    std::cin >> A;
    auto B = A.inv();
    if (B.has_value())
        REP (i, n)
            REP (j, n)
                std::cout << B.value()[i][j] << "\n "[j + 1 < n];
    else
        std::cout << -1 << "\n";
}