tkmst201's Library

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

View the Project on GitHub tkmst201/Library

:heavy_check_mark: Test/matrix.det.test.cpp

Depends on

Code

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

#include "Mathematics/ModInt.hpp"
#include "Mathematics/Matrix.hpp"

#include <iostream>

int main() {
	std::cin.tie(0);
	std::ios::sync_with_stdio(false);
	
	int N;
	std::cin >> N;
	Matrix<ModInt<998244353>> mat(N, N);
	for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) std::cin >> mat[i][j];
	std::cout << mat.det() << std::endl;
}
#line 1 "Test/matrix.det.test.cpp"
#define PROBLEM "https://judge.yosupo.jp/problem/matrix_det"

#line 1 "Mathematics/ModInt.hpp"



#include <cassert>
#include <iostream>
#include <cstdint>

/**
 * @brief https://tkmst201.github.io/Library/Mathematics/ModInt.hpp
 */
template<int M>
struct ModInt {
	static_assert(M > 0);
	
	using value_type = int;
	using calc_type = std::int_fast64_t;
	
private:
	value_type val_;
	
public:
	constexpr ModInt(calc_type val = 0) : val_(val % M + (val >= 0 ? 0 : M)) {}
	constexpr value_type val() const noexcept { return val_; }
	constexpr static decltype(M) mod() noexcept { return M; }
	
	explicit constexpr operator bool() const noexcept { return val_; }
	constexpr bool operator !() const noexcept { return !static_cast<bool>(*this); }
	constexpr ModInt operator +() const noexcept { return *this; }
	constexpr ModInt operator -() const noexcept { return ModInt(val_ == 0 ? 0 : M - val_); }
	constexpr ModInt operator ++(int) noexcept { ModInt res = *this; ++*this; return res; }
	constexpr ModInt operator --(int) noexcept { ModInt res = *this; --*this; return res; }
	constexpr ModInt & operator ++() noexcept { val_ = val_ + 1 == M ? 0 : val_ + 1; return *this; }
	constexpr ModInt & operator --() noexcept { val_ = val_ == 0 ? M - 1 : val_ - 1; return *this; }
	constexpr ModInt & operator +=(const ModInt & rhs) noexcept { val_ += val_ < M - rhs.val_ ? rhs.val_ : rhs.val_ - M; return *this; }
	constexpr ModInt & operator -=(const ModInt & rhs) noexcept { val_ += val_ >= rhs.val_ ? -rhs.val_ : M - rhs.val_; return *this; }
	constexpr ModInt & operator *=(const ModInt & rhs) noexcept { val_ = static_cast<calc_type>(val_) * rhs.val_ % M; return *this; }
	constexpr ModInt & operator /=(const ModInt & rhs) noexcept { return *this *= rhs.inv(); }
	friend constexpr ModInt operator +(const ModInt & lhs, const ModInt & rhs) noexcept { return ModInt(lhs) += rhs; }
	friend constexpr ModInt operator -(const ModInt & lhs, const ModInt & rhs) noexcept { return ModInt(lhs) -= rhs; }
	friend constexpr ModInt operator *(const ModInt & lhs, const ModInt & rhs) noexcept { return ModInt(lhs) *= rhs; }
	friend constexpr ModInt operator /(const ModInt & lhs, const ModInt & rhs) noexcept { return ModInt(lhs) /= rhs; }
	friend constexpr bool operator ==(const ModInt & lhs, const ModInt & rhs) noexcept { return lhs.val_ == rhs.val_; }
	friend constexpr bool operator !=(const ModInt & lhs, const ModInt & rhs) noexcept { return !(lhs == rhs); }
	friend std::ostream & operator <<(std::ostream & os, const ModInt & rhs) { return os << rhs.val_; }
	friend std::istream & operator >>(std::istream & is, ModInt & rhs) { calc_type x; is >> x; rhs = ModInt(x); return is; }
	
	constexpr ModInt pow(calc_type n) const noexcept {
		ModInt res = 1, x = val_;
		if (n < 0) { x = x.inv(); n = -n; }
		while (n) { if (n & 1) res *= x; x *= x; n >>= 1; }
		return res;
	}
	
	constexpr ModInt inv() const noexcept {
		value_type a = val_, a1 = 1, b = M, b1 = 0;
		while (b > 0) {
			const value_type q = a / b;
			value_type tmp = a - q * b; a = b; b = tmp;
			tmp = a1 - q * b1; a1 = b1; b1 = tmp;
		}
		assert(a == 1);
		if (a1 < 0) a1 += M;
		return a1;
	}
};


#line 1 "Mathematics/Matrix.hpp"



#line 6 "Mathematics/Matrix.hpp"
#include <vector>
#include <utility>
#include <type_traits>
#include <initializer_list>
#include <algorithm>
#include <cmath>

/**
 * @brief https://tkmst201.github.io/Library/Mathematics/Matrix.hpp
 */
template<typename T>
struct Matrix {
	using value_type = T;
	using size_type = std::size_t;
	using uint32 = std::uint32_t;
	
private:
	size_type h = 0, w = 0;
	std::vector<std::vector<value_type>> val;
	constexpr static value_type eps = std::is_floating_point<value_type>::value ? 1e-8 : 0;
	
public:
	Matrix() = default;
	
	Matrix(size_type h, size_type w, const value_type & x = 0)
		: h(h), w(w), val(h, std::vector<value_type>(w, x)) {}
	
	explicit Matrix(const std::vector<std::vector<value_type>> & val)
		: h(val.size()), w(val.empty() ? 0 : val[0].size()), val(val) {
		for (size_type i = 1; i < h; ++i) assert(val[i].size() == w);
	}
	
	explicit Matrix(std::initializer_list<std::vector<value_type>> init) : val(init.begin(), init.end()) {
		h = val.size();
		w = val.empty() ? 0 : val[0].size();
		for (size_type i = 1; i < h; ++i) assert(val[i].size() == w);
	}
	
	static Matrix identity(size_type n) {
		Matrix res(n, n);
		for (size_type i = 0; i < n; ++i) res(i, i) = 1;
		return res;
	}
	
	Matrix operator +() const noexcept {
		return *this;
	}
	
	Matrix operator -() const {
		return Matrix(h, w, 0) -= *this;
	}
	
	Matrix & operator +=(const Matrix & rhs) noexcept {
		assert(match_type(rhs));
		for (size_type i = 0; i < h; ++i) for (size_type j = 0; j < w; ++j) val[i][j] += rhs.val[i][j];
		return *this;
	}
	
	Matrix & operator -=(const Matrix & rhs) noexcept {
		assert(match_type(rhs));
		for (size_type i = 0; i < h; ++i) for(size_type j = 0; j < w; ++j) val[i][j] -= rhs.val[i][j];
		return *this;
	}
	
	Matrix & operator *=(const Matrix & rhs) {
		assert(w == rhs.h);
		Matrix mat(h, rhs.w);
		for (size_type i = 0; i < h; ++i) for (size_type j = 0; j < rhs.w; ++j) for (size_type k = 0; k < w; ++k)
			mat.val[i][j] += val[i][k] * rhs.val[k][j];
		h = mat.h;
		w = mat.w;
		val = std::move(mat.val);
		return *this;
	}
	
	Matrix & operator /=(const Matrix & rhs) {
		assert(rhs.is_square());
		assert(!rhs.empty());
		assert(w == rhs.h);
		const Matrix mat = rhs.inv();
		assert(!mat.empty());
		*this *= mat;
		return *this;
	}
	
	friend Matrix operator +(const Matrix & lhs, const Matrix & rhs) {
		return Matrix(lhs) += rhs;
	}
	
	friend Matrix operator -(const Matrix & lhs, const Matrix & rhs) {
		return Matrix(lhs) -= rhs;
	}
	
	friend Matrix operator *(const Matrix & lhs, const Matrix & rhs) {
		return Matrix(lhs) *= rhs;
	}
	
	friend Matrix operator /(const Matrix & lhs, const Matrix & rhs) {
		return Matrix(lhs) /= rhs;
	}
	
	friend Matrix operator *(const value_type & lhs, const Matrix & rhs) {
		Matrix res(rhs.val);
		for (size_type i = 0; i < res.h; ++i) for (size_type j = 0; j < res.w; ++j) res.val[i][j] = lhs * res.val[i][j];
		return res;
	}
	
	friend Matrix operator *(const Matrix & lhs, const value_type & rhs) {
		Matrix res(lhs.val);
		for (size_type i = 0; i < res.h; ++i) for (size_type j = 0; j < res.w; ++j) res.val[i][j] *= rhs;
		return res;
	}
	
	friend Matrix operator /(const Matrix & lhs, const value_type & rhs) {
		if constexpr (std::is_floating_point<value_type>::value) assert(std::abs(rhs) >= eps);
		else assert(rhs != 0);
		Matrix res(lhs.val);
		for (size_type i = 0; i < res.h; ++i) for (size_type j = 0; j < res.w; ++j) res.val[i][j] /= rhs;
		return res;
	}
	
	friend bool operator ==(const Matrix & lhs, const Matrix & rhs) noexcept {
		if (!lhs.match_type(rhs)) return false;
		if constexpr (!std::is_floating_point<value_type>::value) return lhs.val == rhs.val;
		else {
			for (size_type i = 0; i < lhs.h; ++i) for (size_type j = 0; j < lhs.w; ++j) {
				if (std::abs(lhs.val[i][j] - rhs.val[i][j]) >= eps) return false;
			}
			return true;
		}
	}
	
	friend bool operator !=(const Matrix & lhs, const Matrix & rhs) noexcept {
		return !(lhs == rhs);
	}
	
	std::vector<value_type> & operator [](size_type i) noexcept {
		return val[i];
	}
	
	const std::vector<value_type> & operator [](size_type i) const noexcept {
		return val[i];
	};
	
	value_type & at(size_type i, size_type j) noexcept {
		assert(i < h);
		assert(j < w);
		return val[i][j];
	}
	
	const value_type & at(size_type i, size_type j) const noexcept {
		assert(i < h);
		assert(j < w);
		return val[i][j];
	}
	
	value_type & operator ()(size_type i, size_type j) noexcept {
		assert(i < h);
		assert(j < w);
		return val[i][j];
	};
	
	const value_type & operator ()(size_type i, size_type j) const noexcept {
		assert(i < h);
		assert(j < w);
		return val[i][j];
	}
	
	bool empty() const noexcept {
		return h == 0 && w == 0;
	}
	
	std::pair<size_type, size_type> type() const noexcept {
		return {h, w};
	}
	
	bool match_type(const Matrix & A) const noexcept {
		return h == A.h && w == A.w;
	}
	
	bool is_square() const noexcept {
		return h == w;
	}
	
	const std::vector<std::vector<value_type>> & get() const noexcept {
		return val;
	}
	
	Matrix pow(long long n) const {
		assert(h == w);
		assert(n >= 0);
		Matrix res = identity(h), x(*this);
		while (n > 0) { if (n & 1) res *= x; x *= x; n >>= 1; }
		return res;
	}
	
	Matrix trans() const {
		Matrix res(w, h);
		for (size_type i = 0; i < h; ++i) for (size_type j = 0; j < w; ++j) res.val[j][i] = val[i][j];
		return res;
	}
	
	Matrix vstack(const Matrix & A) const {
		assert(w == A.w);
		Matrix res(h + A.h, w);
		std::copy(val.begin(), val.end(), res.val.begin());
		std::copy(A.val.begin(), A.val.end(), res.val.begin() + h);
		return res;
	}
	
	Matrix hstack(const Matrix & A) const {
		assert(h == A.h);
		Matrix res(h, w + A.w);
		for (int i = 0; i < h; ++i) {
			std::copy(val[i].begin(), val[i].end(), res.val[i].begin());
			std::copy(A.val[i].begin(), A.val[i].end(), res.val[i].begin() + w);
		}
		return res;
	}
	
	Matrix submat(size_type i1, size_type i2, size_type j1, size_type j2) const {
		assert(i1 < i2 && i2 <= h);
		assert(j1 < j2 && j2 <= w);
		Matrix res(i2 - i1, j2 - j1);
		for (size_type i = 0; i < i2 - i1; ++i) std::copy(val[i + i1].begin() + j1, val[i + i1].begin() + j2, res.val[i].begin());
		return res;
	}
	
	Matrix inv() const {
		assert(is_square());
		Matrix aug_mat = this->hstack(identity(h));
		if (aug_mat.gauss_jordan(h).first != h) return Matrix();
		return aug_mat.submat(0, h, h, 2 * h);
	}
	
	value_type det() const {
		assert(is_square());
		return Matrix(*this).gauss_jordan(w).second;
	}
	
	std::pair<uint32, value_type> gauss_jordan(uint32 col) noexcept {
		assert(col <= w);
		uint32 rank = 0;
		value_type det = empty() || !is_square() ? 0 : 1;
		bool rflag = false;
		for (uint32 k = 0; k < col; ++k) {
			int pivot = -1;
			if constexpr (std::is_floating_point<value_type>::value) {
				value_type mx = eps;
				for (uint32 i = rank; i < h; ++i) {
					const value_type cur = std::abs(val[i][k]);
					if (mx < cur) mx = cur, pivot = i;
				}
			}
			else {
				for (uint32 i = rank; i < h; ++i) if (val[i][k] != 0) { pivot = i; break; }
			}
			if (pivot == -1) continue;
			if (static_cast<uint32>(pivot) != rank) {
				rflag ^= true;
				std::swap(val[pivot], val[rank]);
			}
			det *= val[rank][k];
			const value_type div = static_cast<value_type>(1) / val[rank][k];
			val[rank][k] = 1;
			for (uint32 j = k + 1; j < w; ++j) val[rank][j] *= div;
			for (uint32 i = 0; i < rank; ++i) {
				for (uint32 j = k + 1; j < w; ++j) val[i][j] -= val[rank][j] * val[i][k];
				val[i][k] = 0;
			}
			for (uint32 i = std::max<uint32>(rank + 1, pivot); i < h; ++i) {
				for (uint32 j = k + 1; j < w; ++j) val[i][j] -= val[rank][j] * val[i][k];
				val[i][k] = 0;
			}
			++rank;
		}
		if (rank != h) det = 0;
		if (rflag) det = -det;
		return {rank, det};
	}
	
	friend std::ostream & operator <<(std::ostream & os, const Matrix & rhs) {
		os << "type = (" << rhs.h << "," << rhs.w << ") [\n";
		for (size_type i = 0; i < rhs.h; ++i) {
			os << ' ';
			for (size_type j = 0; j < rhs.w; ++j) os << rhs.val[i][j] << " \n"[j + 1 == rhs.w];
		}
		return os << "]";
	}
};


#line 5 "Test/matrix.det.test.cpp"

#line 7 "Test/matrix.det.test.cpp"

int main() {
	std::cin.tie(0);
	std::ios::sync_with_stdio(false);
	
	int N;
	std::cin >> N;
	Matrix<ModInt<998244353>> mat(N, N);
	for (int i = 0; i < N; ++i) for (int j = 0; j < N; ++j) std::cin >> mat[i][j];
	std::cout << mat.det() << std::endl;
}
Back to top page