行列構造体を使ってみる

競技プログラミングにおいて、行列演算が活躍する場面はそれなりにあります。
もちろんその都度配列を用意しても良いのですが特に行列の積の計算などは面倒で、やっぱり構造体で定義したくなりました。
二次元の行列で簡単な演算ができるように、構造体を書いてみたので以下に貼ります。
私は普段C++を使っていますがそこまで詳しくはないので、何かおかしな点に気付いた場合は教えてくださると幸いです。
(こちらの記事も参考になると思います)
競技プログラミング用にC++で行列（2次元平面）ライブラリを作成する - Qiita
作ったばかりなのでもしかしたらバグが潜んでいるかもしれないです。もし使用する場合は挙動の確認をお願いします。

ソースコード

template<typename T>
struct matrix{
	int row, column;
	vector<vector<T>> mat;

	matrix(int r = 0, int c = 0, T val = 0) : row(r), column(c) {
		mat = vector<vector<T>>(r, vector<T>(c, val));
	}
	vector<T>& operator [] (const int num) {
		return mat[num];
	}
	matrix& operator += (const matrix A) {
		assert(row == A.row && column == A.column);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < column; j++) {
				mat[i][j] += A.mat[i][j];
			}
		}
		return *this;
	}
	matrix& operator -= (const matrix A) {
		assert(row == A.row && column == A.column);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < column; j++) {
				mat[i][j] -= A.mat[i][j];
			}
		}
		return *this;
	}
	matrix& operator *= (const matrix A) {
		assert(column == A.row);
		matrix B(row, A.column);
		for(int i = 0; i < row; i++) {
			for(int k = 0; k < column; k++) {
				for(int j = 0; j < A.column; j++) {
					B.mat[i][j] += mat[i][k] * A.mat[k][j];
				}
			}
		}
		return *this = B;
	}
	matrix& operator *= (long long k) {
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < column; j++) {
				mat[i][j] *= k;
			}
		}
		return *this;
	}
	matrix operator + (const matrix A) const {
		assert(row == A.row && column == A.column);
		matrix B = *this;
		return B += A;
	}
	matrix operator - (const matrix A) const {
		assert(row == A.row && column == A.column);
		matrix B = *this;
		return B -= A;
	}
	matrix operator * (const matrix A) const {
		assert(column == A.row);
		matrix B = *this;
		return B *= A;
	}
	matrix operator * (long long k) const {
		matrix B = *this;
		return B *= k;
	}
	vector<T> operator * (const vector<int> A) const {
		int vector_size = A.size();
		assert(column == vector_size);
		vector<T> B(row);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < column; j++) {
				B[i] += mat[i][j] * A[j];
			}
		}
		return B;
	}
	matrix pow(long long n) const {
		assert(row == column);
		if(n == 0){
			matrix E(row, column);
			for(int i = 0; i < row; i++) {
				E.mat[i][i] = 1;
			}
			return E;
		}
		matrix A = pow(n >> 1);
		A *= A;
		if(n & 1) A *= *this;
		return A;
	}
};

簡単に説明すると、行列同士の和、差、積の計算、行列の定数倍やベクトルとの積、累乗計算などを用意しました。
行列同士の演算では、行数や列数の関係で定義されない演算をしようとすると実行時エラーが発生するようになっているはずです。
使用例として、この構造体を使ってABC189E - Rotate and Flipを解いてみます。
解説にもあるように各変換はアフィン変換ですので、3 x 3 行列で表せます。こういう時に構造体があると楽そうですね。

ACコード

#include <bits/stdc++.h>
#define rep(i,n) for(int i = 0; i < (int)(n); i++)

using namespace std;
using LL = long long;

template<typename T>
struct matrix{
	int row, column;
	vector<vector<T>> mat;

	matrix(int r = 0, int c = 0, T val = 0) : row(r), column(c) {
		mat = vector<vector<T>>(r, vector<T>(c, val));
	}
	vector<T>& operator [] (const int num) {
		return mat[num];
	}
	matrix& operator += (const matrix A) {
		assert(row == A.row && column == A.column);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < column; j++) {
				mat[i][j] += A.mat[i][j];
			}
		}
		return *this;
	}
	matrix& operator -= (const matrix A) {
		assert(row == A.row && column == A.column);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < column; j++) {
				mat[i][j] -= A.mat[i][j];
			}
		}
		return *this;
	}
	matrix& operator *= (const matrix A) {
		assert(column == A.row);
		matrix B(row, A.column);
		for(int i = 0; i < row; i++) {
			for(int k = 0; k < column; k++) {
				for(int j = 0; j < A.column; j++) {
					B.mat[i][j] += mat[i][k] * A.mat[k][j];
				}
			}
		}
		return *this = B;
	}
	matrix& operator *= (long long k) {
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < column; j++) {
				mat[i][j] *= k;
			}
		}
		return *this;
	}
	matrix operator + (const matrix A) const {
		assert(row == A.row && column == A.column);
		matrix B = *this;
		return B += A;
	}
	matrix operator - (const matrix A) const {
		assert(row == A.row && column == A.column);
		matrix B = *this;
		return B -= A;
	}
	matrix operator * (const matrix A) const {
		assert(column == A.row);
		matrix B = *this;
		return B *= A;
	}
	matrix operator * (long long k) const {
		matrix B = *this;
		return B *= k;
	}
	vector<T> operator * (const vector<T> A) const {
		int vector_size = A.size();
		assert(column == vector_size);
		vector<T> B(row);
		for(int i = 0; i < row; i++) {
			for(int j = 0; j < column; j++) {
				B[i] += mat[i][j] * A[j];
			}
		}
		return B;
	}
	matrix pow(long long n) const {
		assert(row == column);
		if(n == 0){
			matrix E(row, column);
			for(int i = 0; i < row; i++) {
				E.mat[i][i] = 1;
			}
			return E;
		}
		matrix A = pow(n >> 1);
		A *= A;
		if(n & 1) A *= *this;
		return A;
	}
};

int main(){
	int N;
	cin >> N;
	vector<LL> X(N), Y(N);
	rep(i,N) cin >> X[i] >> Y[i];
	int M;
	cin >> M;
	vector<matrix<LL>> cum(M + 1);
	rep(i,M+1) cum[i] = matrix<LL>(3, 3);
	rep(i,3) cum[0][i][i] = 1;
	rep(i,M){
		matrix<LL> A(3, 3);
		int t;
		cin >> t;
		if(t == 1){
			A[0][1] = 1, A[1][0] = -1, A[2][2] = 1;
		}
		else if(t == 2){
			A[0][1] = -1, A[1][0] = 1, A[2][2] = 1;
		}
		else if(t == 3){
			int p;
			cin >> p;
			A[0][0] = -1, A[0][2] = 2 * p, A[1][1] = 1, A[2][2] = 1;
		}
		else{
			int p;
			cin >> p;
			A[0][0] = 1, A[1][1] = -1, A[1][2] = 2 * p, A[2][2] = 1;			
		}
		cum[i+1] = A * cum[i];
	}
	int Q;
	cin >> Q;
	rep(i,Q){
		int a, b;
		cin >> a >> b;
		b--;
		vector<LL> p(3);
		p[0] = X[b], p[1] = Y[b], p[2] = 1;
		vector<LL> ans = cum[a] * p;
		cout << ans[0] << " " << ans[1] << endl;
	}

	return 0;
}