verify/yosupo/yosupo-inverse_matrix.test.cpp

• View this file on GitHub
• Last update: 2022-07-21 16:12:33-04:00
• Problem: https://judge.yosupo.jp/problem/inverse_matrix

Depends on

• library/contest/template-minimal.hpp
• library/modular-arithmetic/mod-int2.hpp
• library/numerical/matrix2.hpp

Code

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

#include "../../library/contest/template-minimal.hpp"
#include "../../library/numerical/matrix2.hpp"
#include "../../library/modular-arithmetic/mod-int2.hpp"

using mi = Mint<998244353, 5>;

int main() {
	using namespace MatrixOperations;
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n;
	cin >> n;
	auto a = make_matrix<mi>(n, n);
	for (int i = 0; i < n; ++i) {
		for (int j = 0; j < n; ++j) {
			cin >> a[i][j];
		}
	}
	auto tmp = a;
	mi det = gauss(tmp).first;
	if (det == 0) {
		cout << -1 << '\n';
	} else {
		auto ans = inv(a);
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < n; ++j) {
				cout << ans[i][j] << ' ';
			}
			cout << '\n';
		}
	}
	return 0;
}

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


#include <algorithm>
#include <array>
#include <bitset>
#include <cassert>
#include <chrono>
#include <cmath>
#include <complex>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <ctime>
#include <deque>
#include <iostream>
#include <iomanip>
#include <list>
#include <map>
#include <numeric>
#include <queue>
#include <random>
#include <set>
#include <stack>
#include <string>
#include <unordered_map>
#include <vector>

using namespace std;

namespace MatrixOperations {

template <class T> using Matrix = std::vector<std::vector<T>>;

template <class T> Matrix<T> make_matrix(int r, int c) { return Matrix<T>(r, std::vector<T>(c)); }

template <class T> Matrix<T>& operator+=(Matrix<T>& a, const Matrix<T>& b) {
	for (int i = 0; i < (int)a.size(); i++) {
		for (int j = 0; j < (int)a[0].size(); j++) {
			a[i][j] += b[i][j];
		}
	}
	return a;
}

template <class T> Matrix<T>& operator-=(Matrix<T>& a, const Matrix<T>& b) {
	for (int i = 0; i < (int)a.size(); i++) {
		for (int j = 0; j < (int)a[0].size(); j++) {
			a[i][j] -= b[i][j];
		}
	}
	return a;
}

template <class T> Matrix<T> operator*(const Matrix<T>& a, const Matrix<T>& b) {
	assert(a[0].size() == b.size());
	int x = (int)a.size();
	int y = (int)a[0].size();
	int z = (int)b[0].size();
	Matrix<T> c = make_matrix<T>(x, z);
	for (int i = 0; i < x; i++) {
		for (int j = 0; j < y; j++) {
			for (int k = 0; k < z; k++) {
				c[i][k] += a[i][j] * b[j][k];
			}
		}
	}
	return c;
}

template <class T> Matrix<T> operator+(Matrix<T> a, const Matrix<T>& b) { return a += b; }
template <class T> Matrix<T> operator-(Matrix<T> a, const Matrix<T>& b) { return a -= b; }
template <class T> Matrix<T>& operator*=(Matrix<T>& a, const Matrix<T>& b) { return a = a * b; }

template <class T> Matrix<T> pow(Matrix<T> m, long long p) {
	int n = (int)m.size();
	assert(n == (int)m[0].size() && p >= 0);
	Matrix<T> res = make_matrix<T>(n, n);
	for (int i = 0; i < n; i++) 
		res[i][i] = 1;
	for (; p; p >>= 1, m *= m) {
		if (p & 1) {
			res *= m;
		}
	}
	return res;
}

template <class T> int get_row(Matrix<T>& m, int r, int i, int nxt) {
	for (int j = nxt; j < r; j++) {
		if (m[j][i] != 0) {
			return j;
		}
	}
	return -1;
}

const long double EPS = 1e-12;

template <> int get_row<long double>(Matrix<long double>& m, int r, int i, int nxt) {
	std::pair<long double, int> best = {0, -1};
	for (int j = nxt; j < r; j++) {
		best = std::max(best, {abs(m[j][i]), j});
	}
	return best.first < EPS ? -1 : best.second;
}

// returns a pair of determinant, rank, while doing Gaussian elimination to m

template <class T> std::pair<T, int> gauss(Matrix<T>& m) {
	int r = (int)m.size();
	int c = (int)m[0].size();
	int rank = 0, nxt = 0;
	T prod = 1;
	for (int i = 0; i < r; i++) {
		int row = get_row(m, r, i, nxt);
		if (row == -1) {
			prod = 0;
			continue;
		}
		if (row != nxt) {
			prod *= -1;
			m[row].swap(m[nxt]);
		}
		prod *= m[nxt][i];
		rank++;
		T x = 1 / m[nxt][i];
		for (int k = i; k < c; k++) 
			m[nxt][k] *= x;
		for (int j = 0; j < r; j++) {
			if (j != nxt) {
				T v = m[j][i];
				if (v == 0) continue;
				for (int k = i; k < c; k++) {
					m[j][k] -= v * m[nxt][k];
				}
			}
		}
		nxt++;
	}
	return {prod, rank};
}

template <class T> Matrix<T> inv(Matrix<T> m) {
	int r = (int)m.size();
	assert(r == (int)m[0].size());
	Matrix<T> x = make_matrix<T>(r, 2 * r);
	for (int i = 0; i < r; i++) {
		x[i][i + r] = 1;
		for (int j = 0; j < r; j++) {
			x[i][j] = m[i][j];
		}
	}
	if (gauss(x).second != r) return Matrix<T>();
	Matrix<T> res = make_matrix<T>(r, r);
	for (int i = 0; i < r; i++) {
		for (int j = 0; j < r; j++) {
			res[i][j] = x[i][j + r];
		}
	}
	return res;
}

} // namespace MatrixOperations

// 5 is a root of both mods

template <int MOD, int RT> struct Mint {
	static const int mod = MOD;
	static constexpr Mint rt() { return RT; } // primitive root for FFT

	static constexpr int md() { return MOD; } // primitive root for FFT

	int v; 
	explicit operator int() const { return v; } // explicit -> don't silently convert to int

	explicit operator bool() const { return v != 0; }
	Mint() { v = 0; }
	Mint(long long _v) { v = int((-MOD <= _v && _v < MOD) ? _v : _v % MOD); if (v < 0) v += MOD; }
	friend bool operator==(const Mint& a, const Mint& b) { return a.v == b.v; }
	friend bool operator!=(const Mint& a, const Mint& b) { return !(a == b); }
	friend bool operator<(const Mint& a, const Mint& b) { return a.v < b.v; }
	friend bool operator>(const Mint& a, const Mint& b) { return a.v > b.v; }
	friend bool operator<=(const Mint& a, const Mint& b) { return a.v <= b.v; }
	friend bool operator>=(const Mint& a, const Mint& b) { return a.v >= b.v; }
	friend std::istream& operator >> (std::istream& in, Mint& a) { 
		long long x; std::cin >> x; a = Mint(x); return in; }
	friend std::ostream& operator << (std::ostream& os, const Mint& a) { return os << a.v; }
	Mint& operator+=(const Mint& m) { 
		if ((v += m.v) >= MOD) v -= MOD; 
		return *this; }
	Mint& operator-=(const Mint& m) { 
		if ((v -= m.v) < 0) v += MOD; 
		return *this; }
	Mint& operator*=(const Mint& m) { 
		v = (long long)v * m.v % MOD; return *this; }
	Mint& operator/=(const Mint& m) { return (*this) *= inv(m); }
	friend Mint pow(Mint a, long long p) {
		Mint ans = 1; assert(p >= 0);
		for (; p; p /= 2, a *= a) if (p & 1) ans *= a;
		return ans; 
	}
	friend Mint inv(const Mint& a) { assert(a.v != 0); return pow(a, MOD - 2); }
	Mint operator-() const { return Mint(-v); }
	Mint& operator++() { return *this += 1; }
	Mint& operator--() { return *this -= 1; }
	friend Mint operator+(Mint a, const Mint& b) { return a += b; }
	friend Mint operator-(Mint a, const Mint& b) { return a -= b; }
	friend Mint operator*(Mint a, const Mint& b) { return a *= b; }
	friend Mint operator/(Mint a, const Mint& b) { return a /= b; }
};

using mi = Mint<998244353, 5>;

int main() {
	using namespace MatrixOperations;
	ios::sync_with_stdio(false);
	cin.tie(nullptr);
	int n;
	cin >> n;
	auto a = make_matrix<mi>(n, n);
	for (int i = 0; i < n; ++i) {
		for (int j = 0; j < n; ++j) {
			cin >> a[i][j];
		}
	}
	auto tmp = a;
	mi det = gauss(tmp).first;
	if (det == 0) {
		cout << -1 << '\n';
	} else {
		auto ans = inv(a);
		for (int i = 0; i < n; ++i) {
			for (int j = 0; j < n; ++j) {
				cout << ans[i][j] << ' ';
			}
			cout << '\n';
		}
	}
	return 0;
}

Back to top page