cp-library
This documentation is automatically generated by competitive-verifier/competitive-verifier
matrix.hpp
- View this file on GitHub
- Last update: 2024-09-22 02:06:53+09:00
- Include:
#include "matrix.hpp"
Code
template <typename T>
struct Matrix : public vector<vector<T>> {
using Base = vector<vector<T>>;
using Base::front;
using Base::size;
using Base::vector;
Matrix(int h, int w) : Matrix(h, vector<T>(w)) {}
int sizeH() const { return size(); }
int sizeW() const { return front().size(); }
Matrix operator+() const { return *this; }
Matrix operator-() const {
Matrix res(*this);
for (auto &e : res) e = -e;
return res;
}
Matrix operator*(T k) const { return Mat(*this) *= k; }
Matrix operator/(T k) const { return Mat(*this) /= k; }
Matrix operator+(const Matrix &mat) const { return Mat(*this) += mat; }
Matrix operator-(const Matrix &mat) const { return Mat(*this) -= mat; }
Matrix operator*(const Matrix &mat) const {
int n = sizeH(), m = sizeW(), l = mat.sizeW();
Matrix res(n, l);
for (int i = 0; i < n; ++i) {
for (int k = 0; k < m; ++k) {
for (int j = 0; j < l; ++j) {
res[i][j] += (*this)[i][k] * mat[k][j];
}
}
}
return res;
}
Matrix operator*=(T k) {
for (auto &e : *this) e *= k;
return *this;
}
Matrix operator/=(T k) {
assert(k != 0);
*this *= 1 / k;
return *this;
}
Matrix &operator+=(const Matrix &mat) {
int h = sizeH(), w = sizeW();
for (int i = 0; i < h; ++i) {
for (int j = 0; j < w; ++j) {
(*this)[i][j] += mat[i][j];
}
}
return *this;
}
Matrix &operator-=(const Matrix &mat) { return *this += -mat; }
Matrix &operator*=(const Matrix &mat) { return *this = (*this) * mat; }
static Matrix identity(int n) {
Matrix res(n, n);
for (int i = 0; i < n; ++i) {
res[i][i] = 1;
}
return res;
}
Matrix pow(long long k) const {
Matrix res = identity(sizeH()), tmp = *this;
while (k > 0) {
if (k & 1) {
res *= tmp;
}
tmp *= tmp;
k >>= 1;
}
return res;
}
T det() const {
Matrix mat(*this);
int n = mat.sizeH();
T prod = 1;
for (int j0 = 0; j0 < n; ++j0) {
bool found = false;
for (int i0 = j0; i0 < n; ++i0) {
if (mat[i0][j0] != 0) {
if (i0 != j0) {
swap(mat[i0], mat[j0]);
prod *= -1;
}
prod *= mat[j0][j0];
T inv = 1 / mat[j0][j0];
for (auto &e : mat[j0]) e *= inv;
for (int i1 = 0; i1 < n; ++i1) {
if (j0 != i1) {
T mul = mat[i1][j0];
for (int j1 = j0; j1 < n; ++j1) {
mat[i1][j1] -= mul * mat[j0][j1];
}
}
}
found = true;
break;
}
}
if (!found) {
return 0;
}
}
return prod;
}
optional<Matrix> inv() const {
int n = (*this).sizeH();
Matrix mat(n, 2 * n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
mat[i][j] = (*this)[i][j];
}
mat[i][i + n] = 1;
}
for (int j0 = 0; j0 < n; ++j0) {
bool found = false;
for (int i0 = j0; i0 < n; ++i0) {
if (mat[i0][j0] != 0) {
if (i0 != j0) {
swap(mat[i0], mat[j0]);
}
T inv = 1 / mat[j0][j0];
for (auto &e : mat[j0]) e *= inv;
for (int i1 = 0; i1 < n; ++i1) {
if (j0 != i1) {
T mul = mat[i1][j0];
for (int j1 = j0; j1 < 2 * n; ++j1) {
mat[i1][j1] -= mul * mat[j0][j1];
}
}
}
found = true;
break;
}
}
if (!found) {
return nullopt;
}
}
Matrix res(n, n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
res[i][j] = mat[i][j + n];
}
}
return res;
}
friend ostream &operator<<(ostream &os, const Matrix &mat) {
for (auto itr1 = mat.begin(), end_itr1 = mat.end(); itr1 != end_itr1;) {
const auto &v = *itr1;
for (auto itr2 = v.begin(), end_itr2 = v.end(); itr2 != end_itr2;) {
os << *itr2;
if (++itr2 != end_itr2) os << " ";
}
if (++itr1 != end_itr1) os << "\n";
}
return os;
}
};#line 1 "matrix.hpp"
template <typename T>
struct Matrix : public vector<vector<T>> {
using Base = vector<vector<T>>;
using Base::front;
using Base::size;
using Base::vector;
Matrix(int h, int w) : Matrix(h, vector<T>(w)) {}
int sizeH() const { return size(); }
int sizeW() const { return front().size(); }
Matrix operator+() const { return *this; }
Matrix operator-() const {
Matrix res(*this);
for (auto &e : res) e = -e;
return res;
}
Matrix operator*(T k) const { return Mat(*this) *= k; }
Matrix operator/(T k) const { return Mat(*this) /= k; }
Matrix operator+(const Matrix &mat) const { return Mat(*this) += mat; }
Matrix operator-(const Matrix &mat) const { return Mat(*this) -= mat; }
Matrix operator*(const Matrix &mat) const {
int n = sizeH(), m = sizeW(), l = mat.sizeW();
Matrix res(n, l);
for (int i = 0; i < n; ++i) {
for (int k = 0; k < m; ++k) {
for (int j = 0; j < l; ++j) {
res[i][j] += (*this)[i][k] * mat[k][j];
}
}
}
return res;
}
Matrix operator*=(T k) {
for (auto &e : *this) e *= k;
return *this;
}
Matrix operator/=(T k) {
assert(k != 0);
*this *= 1 / k;
return *this;
}
Matrix &operator+=(const Matrix &mat) {
int h = sizeH(), w = sizeW();
for (int i = 0; i < h; ++i) {
for (int j = 0; j < w; ++j) {
(*this)[i][j] += mat[i][j];
}
}
return *this;
}
Matrix &operator-=(const Matrix &mat) { return *this += -mat; }
Matrix &operator*=(const Matrix &mat) { return *this = (*this) * mat; }
static Matrix identity(int n) {
Matrix res(n, n);
for (int i = 0; i < n; ++i) {
res[i][i] = 1;
}
return res;
}
Matrix pow(long long k) const {
Matrix res = identity(sizeH()), tmp = *this;
while (k > 0) {
if (k & 1) {
res *= tmp;
}
tmp *= tmp;
k >>= 1;
}
return res;
}
T det() const {
Matrix mat(*this);
int n = mat.sizeH();
T prod = 1;
for (int j0 = 0; j0 < n; ++j0) {
bool found = false;
for (int i0 = j0; i0 < n; ++i0) {
if (mat[i0][j0] != 0) {
if (i0 != j0) {
swap(mat[i0], mat[j0]);
prod *= -1;
}
prod *= mat[j0][j0];
T inv = 1 / mat[j0][j0];
for (auto &e : mat[j0]) e *= inv;
for (int i1 = 0; i1 < n; ++i1) {
if (j0 != i1) {
T mul = mat[i1][j0];
for (int j1 = j0; j1 < n; ++j1) {
mat[i1][j1] -= mul * mat[j0][j1];
}
}
}
found = true;
break;
}
}
if (!found) {
return 0;
}
}
return prod;
}
optional<Matrix> inv() const {
int n = (*this).sizeH();
Matrix mat(n, 2 * n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
mat[i][j] = (*this)[i][j];
}
mat[i][i + n] = 1;
}
for (int j0 = 0; j0 < n; ++j0) {
bool found = false;
for (int i0 = j0; i0 < n; ++i0) {
if (mat[i0][j0] != 0) {
if (i0 != j0) {
swap(mat[i0], mat[j0]);
}
T inv = 1 / mat[j0][j0];
for (auto &e : mat[j0]) e *= inv;
for (int i1 = 0; i1 < n; ++i1) {
if (j0 != i1) {
T mul = mat[i1][j0];
for (int j1 = j0; j1 < 2 * n; ++j1) {
mat[i1][j1] -= mul * mat[j0][j1];
}
}
}
found = true;
break;
}
}
if (!found) {
return nullopt;
}
}
Matrix res(n, n);
for (int i = 0; i < n; ++i) {
for (int j = 0; j < n; ++j) {
res[i][j] = mat[i][j + n];
}
}
return res;
}
friend ostream &operator<<(ostream &os, const Matrix &mat) {
for (auto itr1 = mat.begin(), end_itr1 = mat.end(); itr1 != end_itr1;) {
const auto &v = *itr1;
for (auto itr2 = v.begin(), end_itr2 = v.end(); itr2 != end_itr2;) {
os << *itr2;
if (++itr2 != end_itr2) os << " ";
}
if (++itr1 != end_itr1) os << "\n";
}
return os;
}
};