ACM矩阵相关模板；

矩阵快速幂： 因为结构体已经重载了 * 所以只需在原来普通的快速幂模板的基础上稍微修改一下。

其他操作主要是模拟。

#define  _rep(i, a, b) for(int i = a; i <= b; i++)
#define  _for(i, a, b) for(int i = a; i < b; i++)
#define  M_sel *this
using namespace std;

typedef vector<int> vec_int;
typedef vector<vec_int> mat_int;
typedef long long ll;
const int Maxn = 1010;
const int Maxm = 1010;
const int Mod = 10000;

struct Matrix {
	int n, m;
	mat_int a;
	Matrix(int n = 0, int  m= 0) {
		this->n = n; this->m = m;
		vec_int tmp;
		tmp.assign(m, 0);
		_for(i, 0, n)
			a.push_back(tmp);
	}
	Matrix operator + (const Matrix &b) const {
		Matrix tmp(n, m);
		_for(i, 0, n) 
			_for(j, 0, m)
				tmp.a[i][j] = a[i][j] + b.a[i][j];
		return tmp;
	}
	Matrix operator - (const Matrix &b) const {
		Matrix tmp(n, m);
		_for(i, 0, n)
			_for(j, 0, m)
				tmp.a[i][j] = a[i][j] - b.a[i][j];
		return tmp;
	}
	Matrix operator * (const Matrix &b) const {
		Matrix tmp(n, b.m);
		_for(i, 0, n) 
			_for(j, 0, b.m) 
				_for(k, 0, m) 
					tmp.a[i][j] = (tmp.a[i][j] + a[i][k] * b.a[k][j] % Mod) % Mod;//Mod在外部设置，主要是为了方便求快速幂;
		return tmp;
	}
	Matrix pow(ll b) {
		Matrix tmp(n, m);
		_for(i, 0, n)
			tmp.a[i][i] = 1;
		while (b) {
			//if (b & 1) tmp = (*this) * tmp;
			//(*this) = (*this) * (*this);
			if (b & 1) tmp = M_sel * tmp; //指针和乘号放在一起很难受， 在上面开了个宏定义
			M_sel = M_sel * M_sel;
			b >>= 1;
		}
		return tmp;
	}
};