C++の独習がてら行列計算ライブラリを自作してみました。
なんでこんなことしたかと言うと、普段競プロしてたら他人の作ってくれた便利なmintライブラリを使ったりするのですが、演算子オーバーロードの文法とかをよくわかってないままコピペして使ってたので、理解するがてら自分でも何か自作しようと思ったからです。
実務で使うならコードをヘッダファイルに分けるべきなんですが、競プロで使えたらなーという期待で.cppファイルにベタ書きしてます。
調べたらC++には行列計算ライブラリにEigenというのがあるので、実際に実用するならそっちを使ったほうがいいでしょう。あと後述してますが行列式の計算は重いです。改善の余地ありですね。
コード
#include <bits/stdc++.h>
using namespace std;
namespace MatrixLib {
class print_format_error: public std::exception {
virtual const char* what() const noexcept { return "Error: print() method cannot print these values because of unsupported!"; }
};
class matrix_product_error: public std::exception {
virtual const char* what() const noexcept { return "Error: Impossible to calculate matrix product!"; }
};
class determinant_error: public std::exception {
virtual const char* what() const noexcept { return "Error: Impossible to calculate determinant!"; }
};
class inversion_error: public std::exception {
virtual const char* what() const noexcept { return "Error: Impossible to calculate inversion!"; }
};
template <typename T>
class Matrix {
private:
vector<vector<T>> A;
T __cofactor(const vector<vector<T>>& coA) const {
if (coA.size() == 1) return coA[0][0];
if (coA.size() == 2) {
return coA[0][0]*coA[1][1]-coA[0][1]*coA[1][0];
}
T res = 0;
for (int col2=0; col2<coA.size(); col2++) {
vector<vector<T>> cocoA(coA.size()-1);
for (int row=1; row<coA.size(); row++) {
for (int col=0; col<coA.size(); col++) {
if (col2==col) continue;
cocoA[row-1].emplace_back(coA[row][col]);
}
}
if (!(col2&1)) res += coA[0][col2] * __cofactor(cocoA);
else res -= coA[0][col2] * __cofactor(cocoA);
}
return res;
}
public:
size_t row;
size_t col;
Matrix(size_t row=1, size_t col=1) {
this->row = row;
this->col = col;
this->A.resize(row);
for(size_t i=0; i<row; i++) this->A[i].resize(col);
}
vector<T>& operator[](const size_t i) {
return this->A[i];
}
Matrix& operator+=(Matrix other) {
for(size_t i=0; i<this->row; i++) {
for(size_t j=0; j<this->col; j++) {
this->A[i][j] += other[i][j];
}
}
return *this;
}
Matrix operator+(const Matrix other) const {
return Matrix(*this) += other;
}
Matrix& operator+=(const double a) {
for (size_t i=0; i<this->row; i++) {
for (size_t j=0; j<this->col; j++) {
this->A[i][j] += a;
}
}
return *this;
}
Matrix operator+(const double a) const {
return Matrix(*this) += a;
}
Matrix& operator-=(Matrix other) {
for(size_t i=0; i<this->row; i++) {
for(size_t j=0; j<this->col; j++) {
this->A[i][j] -= other[i][j];
}
}
return *this;
}
Matrix operator-(const Matrix other) const {
return Matrix(*this) -= other;
}
Matrix operator-() const {
Matrix res(this->row, this->col);
for (size_t i=0; i<this->row; i++) {
for (size_t j=0; j<this->col; j++) {
res[i][j] = -this->A[i][j];
}
}
return res;
}
Matrix& operator-=(const double a) {
for (size_t i=0; i<this->row; i++) {
for (size_t j=0; j<this->col; j++) {
this->A[i][j] -= a;
}
}
return *this;
}
Matrix operator-(const double a) const {
return Matrix(*this) -= a;
}
Matrix& operator*=(Matrix other) {
if (!(this->row==this->col && other.row==other.col && this->row==other.row)) throw matrix_product_error();
vector<vector<T>> res(this->row, vector<T>(this->col, 0));
for(size_t i=0; i<other.row; i++) {
for (size_t j=0; j<this->col; j++) {
for (size_t k=0; k<this->col; k++) {
res[i][j] += this->A[i][k]*other[k][j];
}
}
}
this->A = res;
return *this;
}
Matrix operator*(Matrix other) const {
if (this->col!=other.row) throw matrix_product_error();
Matrix<T> res(this->row, other.col);
for (size_t i=0; i<this->row; i++) {
for (size_t j=0; j<other.col; j++) {
for (size_t k=0; k<this->col; k++) {
res[i][j] += this->A[i][k]*other[k][j];
}
}
}
return res;
}
Matrix& operator*=(const double a) {
for (size_t i=0; i<this->row; i++) {
for (size_t j=0; j<this->col; j++) {
this->A[i][j] *= a;
}
}
return *this;
}
Matrix operator*(const double a) const {
return Matrix(*this) *= a;
}
void print() {
string format;
if (typeid(T)==typeid(int)) format = "%*d"s;
else if (typeid(T)==typeid(unsigned int) || typeid(T)==typeid(unsigned short)) format = "%*u"s;
else if (typeid(T)==typeid(long)) format = "%*ld"s;
else if (typeid(T)==typeid(unsigned long)) format = "%*lu"s;
else if (typeid(T)==typeid(long long)) format = "%*lld"s;
else if (typeid(T)==typeid(unsigned long long)) format = "%*llu"s;
else if (typeid(T)==typeid(short)) format = "%*f"s;
else if (typeid(T)==typeid(double)) format = "%*lf"s;
else if (typeid(T)==typeid(long double)) format = "%*LF"s;
else throw print_format_error();
int len = 0;
for (size_t i=0; i<this->row; i++) {
for (size_t j=0; j<this->col; j++) {
string str = to_string(this->A[i][j]);
len = max(len, (int)str.size());
}
}
printf("[[");
for (size_t i=0; i<this->row; i++) {
for (size_t j=0; j<this->col; j++) {
printf(format.c_str(), len, this->A[i][j]);
if (j!=this->col-1) printf(", ");
else printf("]");
}
if (i!=this->row-1) printf(",\n ");
else printf("]\n");
}
}
T det() const {
if (this->row!=this->col) throw determinant_error();
return this->__cofactor(this->A);
}
Matrix dot(const Matrix B) const {
return Matrix(*this) * B;
}
Matrix inv() const {
if (!(this->row==this->col)) throw inversion_error();
size_t N = this->row;
Matrix<T> A = *this;
Matrix<T> invA(N, N);
for (size_t i=0; i<N; i++) {
for (size_t j=0; j<N; j++) {
invA[i][j] = (i==j) ? 1 : 0;
}
}
for (size_t i=0; i<N; i++) {
T buf = 1/A[i][i];
for (size_t j=0; j<N; j++) {
A[i][j] *= buf;
invA[i][j] *= buf;
}
for (size_t j=0; j<N; j++) {
if (i!=j) {
buf = A[j][i];
for (size_t k=0; k<N; k++) {
A[j][k] -= A[i][k]*buf;
invA[j][k] -= invA[i][k]*buf;
}
}
}
}
return invA;
}
};
}
GitHubにあげたコードは以下。
github.com
使い方
要素がdouble型の行列
void test(){
cout << "===test===" << endl;
MatrixLib::Matrix<double> mat1(3,3);
MatrixLib::Matrix<double> mat2(3,3);
MatrixLib::Matrix<double> res;
for (int i=0; i<3; i++) {
for (int j=0; j<3; j++) {
mat1[i][j] = j+i*3;
mat2[i][j] = 1;
}
}
mat1.print();
mat2.print();
cout << "mat1.row: " << mat1.row << endl;
cout << "mat1.col: " << mat1.col << endl;
mat2[0][0] = 20;
mat2.print();
mat1 = mat1 + mat2;
mat1.print();
mat1 += mat2;
mat1.print();
mat1 = mat1 - mat2;
mat1.print();
mat1 -= mat2;
mat1.print();
res = -mat1;
res.print();
res = mat1 * mat2;
res.print();
mat1 *= mat2;
mat1.print();
res = mat1 + 1;
res.print();
mat1 += 1;
mat1.print();
res = mat1 - 1;
res.print();
mat1 -= 1;
mat1.print();
res = mat1 * 2;
res.print();
mat1 *= 2;
mat1.print();
mat1 *= 2.1;
mat1.print();
}
int main(int argc, char const *argv[]){
test();
return 0;
}
要素がint型の行列
void test2(){
cout << "===test2===" << endl;
MatrixLib::Matrix<int> mat1(2,3);
MatrixLib::Matrix<int> mat2(3,2);
mat1[0][0]=2; mat1[0][1]=1; mat1[0][2]=3;
mat1[1][0]=5; mat1[1][1]=6; mat1[1][2]=4;
mat2[0][0]=1; mat2[0][1]=2;
mat2[1][0]=3; mat2[1][1]=2;
mat2[2][0]=4; mat2[2][1]=3;
mat1.print();
mat2.print();
MatrixLib::Matrix<int> mat = mat1 * mat2;
mat.print();
mat *= 2;
mat.print();
mat *= 2.1;
mat.print();
}
行列式の計算は余因子の計算を再帰的にしているのですが、計算量的に重いです。
調べたらもっと効率的で高速なやり方が以下のサイトに書いてありました。
Tech Tips: 行列式の計算
今回は余因子の再帰的な計算のままで放置してます。
行列式の計算確認は以下のサイトが便利でした。
行列式（n次元） - 高精度計算サイト
2x2の正方行列の行列式の計算。
void test_det2() {
cout << "===test_det2===" << endl;
MatrixLib::Matrix<int> mat(2,2);
mat[0][0]=3; mat[0][1]=1;
mat[1][0]=5; mat[1][1]=1;
mat.print();
int d = mat.det();
cout << "Det: " << d << endl;
}
3x3の正方行列の行列式の計算。
void test_det3() {
cout << "===test_det3===" << endl;
MatrixLib::Matrix<int> mat(3,3);
mat[0][0]=3; mat[0][1]=1; mat[0][2]=1;
mat[1][0]=5; mat[1][1]=1; mat[1][2]=3;
mat[2][0]=2; mat[2][1]=0; mat[2][2]=1;
mat.print();
int d = mat.det();
cout << "Det: " << d << endl;
}
4x4の正方行列の行列式の計算。
void test_det4() {
cout << "===test_det4===" << endl;
MatrixLib::Matrix<int> mat(4,4);
mat[0][0]=3; mat[0][1]=1; mat[0][2]=1; mat[0][3]=2;
mat[1][0]=5; mat[1][1]=1; mat[1][2]=3; mat[1][3]=4;
mat[2][0]=2; mat[2][1]=0; mat[2][2]=1; mat[2][3]=0;
mat[3][0]=1; mat[3][1]=3; mat[3][2]=2; mat[3][3]=1;
mat.print();
int d = mat.det();
cout << "Det: " << d << endl;
}
5x5の正方行列の行列式の計算。
void test_det5() {
cout << "===test_det5===" << endl;
MatrixLib::Matrix<int> mat(5,5);
mat[0][0]=3; mat[0][1]=1; mat[0][2]=1; mat[0][3]=2; mat[0][4]=1;
mat[1][0]=5; mat[1][1]=1; mat[1][2]=3; mat[1][3]=4; mat[1][4]=1;
mat[2][0]=2; mat[2][1]=0; mat[2][2]=1; mat[2][3]=0; mat[2][4]=1;
mat[3][0]=1; mat[3][1]=3; mat[3][2]=2; mat[3][3]=1; mat[3][4]=1;
mat[4][0]=1; mat[4][1]=1; mat[4][2]=1; mat[4][3]=1; mat[4][4]=1;
mat.print();
int d = mat.det();
cout << "Det: " << d << endl;
}
ドット積の計算
行列のドット積の計算。やってることは*演算子と同じ。Pythonのnumpyライクに使えたら便利かなと思って作成。
void test_dot() {
cout << "===test_dot===" << endl;
MatrixLib::Matrix<int> mat1(3,3);
MatrixLib::Matrix<int> mat2(3,3);
MatrixLib::Matrix<int> res;
mat1[0][0]=2; mat1[0][1]=1; mat1[0][2]=3;
mat1[1][0]=5; mat1[1][1]=6; mat1[1][2]=4;
mat1[2][0]=2; mat1[2][1]=1; mat1[2][2]=3;
mat2[0][0]=1; mat2[0][1]=2; mat2[0][2]=3;
mat2[1][0]=3; mat2[1][1]=2; mat2[1][2]=2;
mat2[2][0]=4; mat2[2][1]=3; mat2[2][2]=1;
mat1.print();
mat2.print();
res = mat1.dot(mat2);
res.print();
mat1.print();
mat2.print();
}
逆行列の求め方は以下のサイトを参考にしました（掃き出し法）。
行列式の値の求め方、逆行列の作り方の C 言語プログラム | Blue Sky Lab
void test_inv() {
cout << "===test_inv===" << endl;
MatrixLib::Matrix<double> mat(4,4);
mat[0][0]=3; mat[0][1]=1; mat[0][2]=1; mat[0][3]=2;
mat[1][0]=5; mat[1][1]=1; mat[1][2]=3; mat[1][3]=4;
mat[2][0]=2; mat[2][1]=0; mat[2][2]=1; mat[2][3]=0;
mat[3][0]=1; mat[3][1]=3; mat[3][2]=2; mat[3][3]=1;
mat.print();
MatrixLib::Matrix<double> inv_mat = mat.inv();
inv_mat.print();
}
