test/Matrix/inverse.test.cpp

• View this file on GitHub
• Last update: 2025-07-04 15:57:54+09:00

Depends on

• Matrix.cpp
• modint.cpp
• template.cpp

Code

#include "../../Matrix.cpp"
#include "../../modint.cpp"
#include "../../template.cpp"
int main() {
    int n;
    cin >> n;
    Matrix<modint<998244353>> a(n);
    cin >> a;
    if (a.determinant() == 0)
        cout << -1 << endl;
    else
        cout << a.inverse() << endl;
}

#line 1 "Matrix.cpp"
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
template <class T>
struct Matrix {
    vector<vector<T>> m;
    Matrix() = default;
    Matrix(int x) : m(vector(x, vector<T>(x, 0))) {}
    Matrix(const vector<vector<T>> &a) : m(a) {}
    vector<T> &operator[](int i) { return m[i]; }
    const vector<T> &operator[](int i) const { return m[i]; }
    static Matrix identity(int x) {
        Matrix res(x);
        for (int i = 0; i < x; i++) res[i][i] = 1;
        return res;
    }
    static Matrix zero(int x) { return Matrix(x); }

    Matrix operator+() const { return (*this); }
    Matrix operator-() const { return Matrix(this->m.size()) - (*this); }
    Matrix operator+(const Matrix &a) const {
        Matrix x = (*this);
        return x += a;
    }
    Matrix operator-(const Matrix &a) const {
        Matrix x = (*this);
        return x -= a;
    }
    Matrix operator*(const Matrix &a) const {
        Matrix x = (*this);
        return x *= a;
    }
    Matrix operator*(const T &a) const {
        Matrix x = (*this);
        return x *= a;
    }
    Matrix &operator+=(const Matrix &a) {
        int n = m.size();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                m[i][j] += a[i][j];
            }
        }
        return *this;
    }
    Matrix &operator-=(const Matrix &a) {
        int n = m.size();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                m[i][j] -= a[i][j];
            }
        }
        return *this;
    }
    Matrix &operator*=(const Matrix &a) {
        int n = m.size();
        Matrix res(n);
        for (int i = 0; i < n; i++) {
            for (int k = 0; k < n; k++) {
                for (int j = 0; j < n; j++) {
                    res[i][j] += m[i][k] * a[k][j];
                }
            }
        }
        m = res.m;
        return *this;
    }
    Matrix &operator*=(const T &a) {
        int n = m.size();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                m[i][j] *= a;
            }
        }
        return *this;
    }
    Matrix pow(ll b) const {
        Matrix x = *this, res = identity(x.m.size());
        while (b) {
            if (b & 1) {
                res *= x;
            }
            x *= x;
            b >>= 1;
        }
        return res;
    }
    T determinant() {
        int n = m.size();
        Matrix A(*this);
        T ret = 1;
        for (int i = 0; i < n; i++) {
            int pivot = -1;
            for (int j = i; j < n; j++) {
                if (A[j][i] != 0) pivot = j;
            }
            if (pivot == -1) return T(0);
            if (i != pivot) {
                ret *= -1;
                swap(A[i], A[pivot]);
            }
            ret *= A[i][i];
            T tmp = T(1) / A[i][i];
            for (int j = 0; j < n; j++) {
                A[i][j] *= tmp;
            }
            for (int j = i + 1; j < n; j++) {
                T a = A[j][i];
                for (int k = 0; k < n; k++) {
                    A[j][k] -= A[i][k] * a;
                }
            }
        }
        return ret;
    }
    Matrix inverse() {
        assert(determinant() != 0);
        int n = m.size();
        Matrix ret = identity(n);
        Matrix A(*this);
        for (int i = 0; i < n; i++) {
            int pivot = -1;
            for (int j = i; j < n; j++) {
                if (A[j][i] != 0) pivot = j;
            }
            if (i != pivot) {
                swap(ret[i], ret[pivot]);
                swap(A[i], A[pivot]);
            }
            T tmp = T(1) / A[i][i];
            for (int j = 0; j < n; j++) {
                A[i][j] *= tmp;
                ret[i][j] *= tmp;
            }
            for (int j = 0; j < n; j++) {
                if (j == i) continue;
                T a = A[j][i];
                for (int k = 0; k < n; k++) {
                    A[j][k] -= A[i][k] * a;
                    ret[j][k] -= ret[i][k] * a;
                }
            }
        }
        return ret;
    }
    vector<T> characteristic_polynomial() {
        int n = m.size();
        if (n == 0) return {1};
        Matrix A(*this);
        for (int i = 1; i < n; i++) {
            int pivot = -1;
            for (int j = i; j < n; j++) {
                if (A[j][i - 1] != 0) pivot = j;
            }
            if (pivot == -1) continue;
            if (i != pivot) {
                swap(A[i], A[pivot]);
                for (int j = 0; j < n; j++) swap(A[j][i], A[j][pivot]);
            }
            T tmp = T(1) / A[i][i - 1];
            vector<T> c(n);
            for (int j = i + 1; j < n; j++) c[j] = tmp * A[j][i - 1];
            for (int j = i + 1; j < n; j++) {
                for (int k = i - 1; k < n; k++) {
                    A[j][k] -= A[i][k] * c[j];
                }
            }
            for (int j = 0; j < n; j++) {
                for (int k = i + 1; k < n; k++) {
                    A[j][i] += A[j][k] * c[k];
                }
            }
        }
        vector dp(n + 1, vector<T>(n + 1));
        dp[0][0] = 1;
        for (int i = 0; i < n; i++) {
            for (int j = 0; j <= i; j++) {
                dp[i + 1][j] -= dp[i][j] * A[i][i];
                dp[i + 1][j + 1] += dp[i][j];
            }
            T p = 1;
            for (int k = i + 1; k < n; k++) {
                p *= A[k][k - 1];
                for (int j = 0; j <= i; j++) dp[k + 1][j] -= dp[i][j] * p * A[i][k];
            }
        }
        return dp[n];
    }
    friend ostream &operator<<(ostream &os, const Matrix &a) {
        int n = a.m.size();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                os << a[i][j];
                if (j != n - 1) os << " ";
            }
            os << "\n";
        }
        return os;
    }
    friend istream &operator>>(istream &is, Matrix &a) {
        int n = a.m.size();
        for (int i = 0; i < n; i++) {
            for (int j = 0; j < n; j++) {
                is >> a[i][j];
            }
        }
        return is;
    }
};
#line 1 "modint.cpp"

#line 3 "modint.cpp"
using namespace std;
template <int modulo>
struct modint {
    int x;
    modint() : x(0) {}
    modint(int64_t y) : x(y >= 0 ? y % modulo : (modulo - (-y) % modulo) % modulo) {}
    modint &operator+=(const modint &p) {
        if ((x += p.x) >= modulo) x -= modulo;
        return *this;
    }
    modint &operator-=(const modint &p) {
        if ((x += modulo - p.x) >= modulo) x -= modulo;
        return *this;
    }
    modint &operator*=(const modint &p) {
        x = (int)(1LL * x * p.x % modulo);
        return *this;
    }
    modint &operator/=(const modint &p) {
        *this *= p.inv();
        return *this;
    }
    modint operator-() const { return modint(-x); }
    modint operator+(const modint &p) const { return modint(*this) += p; }
    modint operator-(const modint &p) const { return modint(*this) -= p; }
    modint operator*(const modint &p) const { return modint(*this) *= p; }
    modint operator/(const modint &p) const { return modint(*this) /= p; }
    bool operator==(const modint &p) const { return x == p.x; }
    bool operator!=(const modint &p) const { return x != p.x; }
    modint inv() const {
        int a = x, b = modulo, u = 1, v = 0, t;
        while (b > 0) {
            t = a / b;
            swap(a -= t * b, b);
            swap(u -= t * v, v);
        }
        return modint(u);
    }
    modint pow(int64_t n) const {
        modint ret(1), mul(x);
        while (n > 0) {
            if (n & 1) ret *= mul;
            mul *= mul;
            n >>= 1;
        }
        return ret;
    }
    friend ostream &operator<<(ostream &os, const modint &p) { return os << p.x; }
    friend istream &operator>>(istream &is, modint &a) {
        int64_t t;
        is >> t;
        a = modint<modulo>(t);
        return (is);
    }
    int val() const { return x; }
    static constexpr int mod() { return modulo; }
    static constexpr int half() { return (modulo + 1) >> 1; }
};
#line 2 "template.cpp"
#include <ext/pb_ds/assoc_container.hpp>
#include <ext/pb_ds/tree_policy.hpp>
using namespace std;
using namespace __gnu_pbds;
using ll = long long;
template <class T>
using pbds_set = tree<T, null_type, less<T>, rb_tree_tag, tree_order_statistics_node_update>;
template <class T>
using pbds_mset = tree<T, null_type, less_equal<T>, rb_tree_tag, tree_order_statistics_node_update>;
using pbds_trie = trie<string, null_type, trie_string_access_traits<>, pat_trie_tag, trie_prefix_search_node_update>;
#define rep(i, n) for (int i = 0; i < n; i++)
#define all(v) v.begin(), v.end()
template <class T, class U>
inline bool chmax(T &a, U b) {
    if (a < b) {
        a = b;
        return true;
    }
    return false;
}
template <class T, class U>
inline bool chmin(T &a, U b) {
    if (a > b) {
        a = b;
        return true;
    }
    return false;
}
constexpr int INF = 1000000000;
constexpr int64_t llINF = 3000000000000000000;
constexpr double eps = 1e-10;
const double pi = acos(-1);
template <class T>
inline void compress(vector<T> &a) {
    sort(a.begin(), a.end());
    a.erase(unique(a.begin(), a.end()), a.end());
}
struct linear_sieve {
    vector<int> least_factor, prime_list;
    linear_sieve(int n) : least_factor(n + 1, 0) {
        for (int i = 2; i <= n; i++) {
            if (least_factor[i] == 0) {
                least_factor[i] = i;
                prime_list.push_back(i);
            }
            for (int p : prime_list) {
                if (ll(i) * p > n || p > least_factor[i]) break;
                least_factor[i * p] = p;
            }
        }
    }
};
ll extgcd(ll a, ll b, ll &x, ll &y) {
    // ax+by=gcd(|a|,|b|)
    if (a < 0 || b < 0) {
        ll d = extgcd(abs(a), abs(b), x, y);
        if (a < 0) x = -x;
        if (b < 0) y = -y;
        return d;
    }
    if (b == 0) {
        x = 1;
        y = 0;
        return a;
    }
    ll d = extgcd(b, a % b, y, x);
    y -= a / b * x;
    return d;
}
ll modpow(ll a, ll b, ll m) {
    ll res = 1;
    while (b) {
        if (b & 1) {
            res *= a;
            res %= m;
        }
        a *= a;
        a %= m;
        b >>= 1;
    }
    return res;
}
template <typename T, typename U>
inline istream &operator>>(istream &is, pair<T, U> &rhs) {
    return is >> rhs.first >> rhs.second;
}
template <typename T>
inline istream &operator>>(istream &is, vector<T> &v) {
    for (auto &e : v) is >> e;
    return is;
}
template <typename T, typename U>
inline ostream &operator<<(ostream &os, const pair<T, U> &rhs) {
    return os << rhs.first << " " << rhs.second;
}
template <typename T>
inline ostream &operator<<(ostream &os, const vector<T> &v) {
    for (auto itr = v.begin(), end_itr = v.end(); itr != end_itr;) {
        os << *itr;
        if (++itr != end_itr) os << " ";
    }
    return os;
}
#line 4 "test/Matrix/inverse.test.cpp"
int main() {
    int n;
    cin >> n;
    Matrix<modint<998244353>> a(n);
    cin >> a;
    if (a.determinant() == 0)
        cout << -1 << endl;
    else
        cout << a.inverse() << endl;
}

Back to top page