行列
(Mathematics/Matrix.hpp)

View this file on GitHub
Last update: 2021-05-02 14:09:01+09:00
Include: #include "Mathematics/Matrix.hpp"
Link: https://tkmst201.github.io/Library/Mathematics/Matrix.hpp

概要

行列を管理する構造体です。
基本的な四則演算やスカラー倍、std::cout に対応しています。

Matrix()
- $\Theta(1)$ 大きさ $(0, 0)$ の空の行列で初期化
Matrix(size_t h, size_t w, const T & x = 0)
- $\Theta(hw)$ 大きさ $(h, w)$ の行列で初期化
Matrix(const vector<vector<T>> & val)
- val の大きさを $(h, w)$ として $\Theta(hw)$ $2$ 次元 vector val で初期化
Matrix(initializer_list<vector<T>> init)
- init の大きさを $(h, w)$ として $\Theta(hw)$ initializer-list init で初期化
static Matrix identity(size_t n)
- $\Theta(n^2)$ 大きさ $(n, n)$ の単位行列を返す
(const) T & at(size_t i, size_t j)
- $\Theta(1)$ $(i, j)$ 成分の参照を返す
bool empty()
- $\Theta(1)$ 空行列であるか判定
pair<size_t, size_t> type()
- $\Theta(1)$ 行列の大きさ $(h, w)$ を返す
bool match_type(const Matrix & A)
- $\Theta(1)$ 行列 A と大きさが等しいか判定
bool is_square()
- $\Theta(1)$ 行列が正方行列であるか判定
const vector<vector<T>> & get()
- $\Theta(1)$ 内部で扱っている $2$ 次元 vector 配列を返す
Matrix pow(long long n)
- $\Theta(h^3 \log{n})$ 現在の行列を $n$ 乗した行列を返す
Matrix trans()
- $\Theta(hw)$ 現在の行列を転置した行列を返す
Matrix vstack(const Matrix & A)
- $\Theta((h + A_h) w)$ 現在の行列の下に行列 A (大きさ $(A_h, A_w)$ ) をつなげた行列を返す
Matrix hstack(const Matrix & A)
- $\Theta(h (w + A_w))$ 現在の行列の右に行列 A (大きさ $(A_h, A_w)$ ) をつなげた行列を返す
Matrix submat(size_t i1, size_t i2, size_t j1, size_t j2)
- $\Theta((i_2 - i_1)(j_2 - j_1))$ 0-indexed で現在の行列の $i_1$ 行目から $i_2-1$ 行目まで、$j_1$ 列目から $j_2-1$ 列目までの部分行列を返す
Matrix inv()
- $\mathcal{O}(h^3)$ 現在の行列の逆行列を返す
Matrix det()
- $\mathcal{O}(h^3)$ 現在の行列の行列式を返す
pair<uint32_t, T> gauss_jordan(uint32_t col)
- $\mathcal{O}(hw$ col$)$ 破壊的な操作で $0$ 列目から $col - 1$ 列目まで行基本変形でガウスの消去法を行い、行列の rank と行列式の組を返す

コンストラクタ

要素の型を T に指定します。

制約

T は double, ModInt

Matrix()

大きさ $(0, 0)$ の空の行列で初期化します。

計算量

$\Theta(1)$

Matrix(size_t h, size_t w, const T & x = 0)

大きさ $(h, w)$ の行列で初期化します。はじめ、すべての要素は $x$ です。

計算量

$\Theta(hw)$

Matrix(const vector<vector<T>> & val)

$2$ 次元 vector val で初期化します。

計算量

val の大きさを $(h, w)$ として $\Theta(hw)$

Matrix(initializer_list<vector<T>> init)

initializer-list init で初期化します。

計算量

init の大きさを $(h, w)$ として $\Theta(hw)$

演算

以下の演算が使用可能です。

+Matrix
-Matrix

Matrix += Matrix
Matrix -= Matrix
Matrix *= Matrix
Matrix /= Matrix

Matrix + Matrix
Matrix - Matrix
Matrix * Matrix
Matrix / Matrix

Matrix * T
T * Matrix
Matrix / T

Matrix == Matrix
Matrix != Matrix

Matrix[i][j] // (i, j) 成分
Matrix(i, j) // (i, j) 成分

std::cout << Matrix

各演算は、その演算が適用できるような行列の大きさの関係式を満たす必要があります。例えば Matrix(lhs) += Matrix(rhs) の場合、lhs 、rhs の大きさは一致している必要があります。

計算量

左の項、右の項の行列の大きさをそれぞれ $(h_1, w_1), (h_2, w_2)$ とします。
また、要素 T の演算に $f$ の計算量がかかる場合は f を掛け合わせてください。

$\Theta(h_1 w_1 w_2)$ ( Matrix *= Matrix または Matrix * Matrix )
$\mathcal{O}((h_1 + w_1) w_1^2)$ ( Matrix /= Matrix または Matrix / Matrix )
$\Theta(1)$ ( Matrix[i][j] または Matrix(i, j) )
$\Theta(h_1 w_1)$ (それ以外)

メンバ関数

以下、行列の大きさを $(h, w)$ とします。

static Matrix identity(size_t n)

大きさ $(n, n)$ の単位行列を返します。

制約

$0 \leq i < n$

計算量

$\Theta(n^2)$

(const) T & at(size_t i, size_t j)

$(i, j)$ 成分の参照を返します。

制約

$0 \leq i < h$
$0 \leq j < w$

計算量

$\Theta(1)$

bool empty()

空行列であるか判定します。空なら $true$ を、空でないなら $false$ を返します。

計算量

$\Theta(1)$

pair<size_t, size_t> type()

行列の大きさ $(h, w)$ を返します。

計算量

$\Theta(1)$

bool match_type(const Matrix & A)

行列 A と大きさが等しいか判定します。等しいなら $true$ を、等しくないなら $false$ を返します。

計算量

$\Theta(1)$

bool is_square()

行列が正方行列であるか判定します。正方行列なら $true$ を、そうでないなら $false$ を返します。

計算量

$\Theta(1)$

const vector<vector<T>> & get()

内部で扱っている $2$ 次元 vector 配列を返します。

計算量

$\Theta(1)$

Matrix pow(long long n)

現在の行列を $n$ 乗した行列を返します。

制約

行列は正方行列 ($h = w$)
$0 \leq n$

計算量

$\Theta(h^3 \log{n})$

Matrix trans()

現在の行列を転置した行列を返します。

計算量

$\Theta(hw)$

Matrix vstack(const Matrix & A)

現在の行列の下に行列 A (大きさ $(A_h, A_w)$ ) をつなげた行列を返します。

制約

$w = A_w$

計算量

$\Theta((h + A_h) w)$

Matrix hstack(const Matrix & A)

現在の行列の右に行列 A (大きさ $(A_h, A_w)$ ) をつなげた行列を返します。

制約

$h = A_h$

計算量

$\Theta(h (w + A_w))$

Matrix submat(size_t i1, size_t i2, size_t j1, size_t j2)

0-indexed で現在の行列の $i_1$ 行目から $i_2-1$ 行目まで、$j_1$ 列目から $j_2-1$ 列目まで部分行列を返します。

制約

$0 \leq i_1 < i_2 \leq h$
$0 \leq j_1 < j_2 \leq w$

計算量

$\Theta((i_2 - i_1)(j_2 - j_1))$

Matrix inv()

現在の行列の逆行列を返します。正則でない (逆行列が存在しない) 場合は空の行列を返します。

制約

行列は正方行列 ($h = w$)

計算量

$\mathcal{O}(h^3)$

Matrix det()

現在の行列の行列式を返します。

制約

行列は正方行列 ($h = w$)

計算量

$\mathcal{O}(h^3)$

pair gauss_jordan(uint32_t col)

$0$ 列目から $col - 1$ 列目まで行基本変形でガウスの消去法を行い、行列の rank と行列式の組を返します。行列が正方行列でないときの行列式の値は $0$ です。
行列の値を変更する破壊的な操作であることに注意してください。

制約

$0 \leq$ col $ \leq w$

計算量

$\mathcal{O}(hw$ col$)$

使用例

#include <bits/stdc++.h>
#include "Mathematics/Matrix.hpp"
#include "Mathematics/ModInt.hpp"
using namespace std;

using mint = ModInt<103>;

int main() {
	Matrix<mint> mat1(2, 2);
	mat1[0][0] = 1;
	mat1.at(0, 1) = 2;
	mat1(1, 0) = 3;
	mat1[1][1] = 4;
	// (2, 2) [[1, 2], [3, 4]]
	cout << mat1 << endl;
	Matrix<mint> mat2({ {1, 2}, {3, 4} });
	// (2, 2) [[1, 2], [3, 4]]
	cout << mat2 << endl;
	cout << "mat1 == mat2 : " << boolalpha << (mat1 == mat2) << endl; // true
	
	cout << mat1 + mat2 << endl; // (2, 2) [[2, 4], [6, 8]]
	cout << mat1 - mat2 << endl; // (2, 2) [[0, 0], [0, 0]]
	cout << mat1 * mat2 << endl; // (2, 2) [[7, 10], [15, 22]]
	cout << mat1 / mat2 << endl; // (2, 2) [[1, 0], [0, 1]]
	cout << mat1 * 2 << endl; // (2, 2) [[2, 4], [6, 8]]
	cout << 3 * mat1 << endl; // (2, 2) [[3, 6], [9, 12]]
	cout << mat1 / 2 << endl; // (2, 2) [[52, 1], [53, 2]]
	cout << mat1[1][1] << endl; // 4
	
	cout << mat1.pow(3) << endl; // (2, 2) [[37, 54], [81, 15]]
	cout << mat1.det() << endl; // 101
	cout << mat1.trans() << endl; // (2, 2) [[1, 3], [2, 4]]
	
	cout << Matrix<mint>::identity(3) << endl; // (3, 3) [[1, 0, 0], [0, 1, 0], [0, 0, 1]]
	cout << "mat1 empty() = " << boolalpha << mat1.empty() << endl; // false
	cout << "empty() = " << boolalpha << Matrix<mint>().empty() << endl; // true
	
	vector v(2, vector<mint>(3, 2));
	Matrix mat3(v);
	cout << mat3 << endl; // (2, 3) [[2, 2, 2], [2, 2, 2]]
	cout << "mat3 type() = (" << mat3.type().first << ", " << mat3.type().second << ")" << endl; // (2, 3)
	cout << "mat1 match_type(mat3) = " << boolalpha << mat1.match_type(mat3) << endl; // false
	cout << "mat1 match_type(mat2) = " << boolalpha << mat1.match_type(mat2) << endl; // true
	cout << "mat3 is_square() = " << boolalpha << mat3.is_square() << endl; // false
	
	mat3[1][2] = 0;
	cout << mat3.trans() << endl; // (3, 2) [[2, 2], [2, 2], [2, 0]]
	
	auto v2 = mat3.get();
	cout << "v == v2 : " << boolalpha << (v == v2) << endl; // false
	
	// (5, 2) [[1, 2], [3, 4], [2, 2], [2, 2], [2, 0]]
	cout << "mat1 vstack(mat3.trans())" << mat1.vstack(mat3.trans()) << endl;
	// (2, 5) [[1, 2, 2, 2, 2], [3, 4, 2, 2, 0]]
	cout << "mat1 hstack(mat3)" << mat1.hstack(mat3) << endl;
	// (2, 2) [[2, 2], [2, 0]]
	cout << "mat3 submat(0, 2, 1, 3) = " << mat3.submat(0, 2, 1, 3) << endl;
	cout << "mat1 inv() = " << mat1.inv() << endl; // (2, 2) [[101, 1], [53, 51]]
	cout << mat1.inv() * mat1 << endl; // (2, 2) [[1, 0], [0, 1]]
	cout << Matrix<mint>(2, 2, 2).inv() << endl; // (0, 0) []
	
	cout << "mat1 gauss_jordan(1)" << endl;
	mat1.gauss_jordan(1); // (2, 2) [[1, 2], [0, 101]]
	cout << mat1 << endl;
	mat1.gauss_jordan(2);
	cout << mat1 << endl; // (2, 2) [[1, 0], [0, 1]]
	
	Matrix<double> matd({ {1e-15, 0}, {0, 0} });
	cout << (Matrix<double>(2, 2, 0) == matd) << endl; // true (誤差を許容する)
}

参考

2021/03/20: https://drken1215.hatenablog.com/entry/2019/03/20/202800

Verified with

Code

#ifndef INCLUDE_GUARD_MATRIX_HPP
#define INCLUDE_GUARD_MATRIX_HPP

#include <iostream>
#include <cassert>
#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 << "]";
	}
};

#endif // INCLUDE_GUARD_MATRIX_HPP

#line 1 "Mathematics/Matrix.hpp"



#include <iostream>
#include <cassert>
#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 << "]";
	}
};

行列 (Mathematics/Matrix.hpp)

概要

コンストラクタ

Matrix()

Matrix(size_t h, size_t w, const T & x = 0)

Matrix(const vector<vector<T>> & val)

Matrix(initializer_list<vector<T>> init)

演算

メンバ関数

static Matrix identity(size_t n)

(const) T & at(size_t i, size_t j)

bool empty()

pair<size_t, size_t> type()

bool match_type(const Matrix & A)

bool is_square()

const vector<vector<T>> & get()

Matrix pow(long long n)

Matrix trans()

Matrix vstack(const Matrix & A)

Matrix hstack(const Matrix & A)

Matrix submat(size_t i1, size_t i2, size_t j1, size_t j2)

Matrix inv()

Matrix det()

pair gauss_jordan(uint32_t col)

使用例

参考

Verified with

Code

行列
(Mathematics/Matrix.hpp)