(Relax 矩阵快速幂 1.2)POJ 3233 Matrix Power Series(用矩阵加法+矩阵快速幂来求sum= A + A2 + A3 + … + Ak)

#include <iostream>
#include <cstring>
#include <cstdio>

using namespace std;

#define maxn 35
int n, k, m;//A是n*n阶的矩阵,k是要求的sum = A^1 + A^2 +.....A^k; m 是被模的元素
struct Mat {
	int val[maxn][maxn];
	void unit() { //单位矩阵
		for (int i = 0; i < maxn; i++)
			val[i][i] = 1;
	}
	void zero() { //零矩阵
		memset(val, 0, sizeof(val));
	}
} x;

Mat operator *(const Mat &a, const Mat &b) { //矩阵乘法
	Mat tmp;
	tmp.zero();
	for (int k = 1; k <= n; k++) {
		for (int i = 1; i <= n; i++)
			if (a.val[i][k])
				for (int j = 1; j <= n; j++) {
					tmp.val[i][j] += a.val[i][k] * b.val[k][j];
					if (tmp.val[i][j] >= m)
						tmp.val[i][j] %= m;
				}
	}
	return tmp;
}

Mat operator ^(Mat x, int n) { //矩阵快速幂
	Mat tmp;
	tmp.zero();
	tmp.unit();
	while (n) {
		if (n & 1)
			tmp = tmp * x;
		x = x * x;
		n >>= 1;
	}
	return tmp;
}

Mat operator +(const Mat &a, const Mat &b) { //矩阵加法
	Mat tmp;
	for (int i = 1; i <= n; i++)
		for (int j = 1; j <= n; j++)
			tmp.val[i][j] = (a.val[i][j] + b.val[i][j]) % m;
	return tmp;
}

Mat sum(int k) {//求sum S = A + A2 + A3 + … + Ak.
	if (k == 1)
		return x;
	else {
		Mat tmp = sum(k >> 1);
		if (k & 1) {
			Mat tmp2 = x ^ ((k >> 1) + 1);
			return tmp + tmp2 + tmp * tmp2;
		} else {
			Mat tmp2 = x ^ (k >> 1);
			return tmp + tmp * tmp2;
		}
	}
}

int main(){
	while(scanf("%d%d%d",&n,&k,&m)!=EOF){
		int i,j;
		for(i = 1 ; i <= n ; ++i){
			for(j = 1 ; j <= n ; ++j){
				scanf("%d",&x.val[i][j]);

				x.val[i][j] %= m;
			}
		}

		Mat ans = sum(k);
		for(i = 1 ; i <= n ; ++i){
			for(j = 1 ; j < n ; ++j){
				printf("%d ",ans.val[i][j]);
			}
			printf("%d\n",ans.val[i][n]);
		}
	}

	return 0;
}