POJ3233 矩阵快速幂

这道题考查的点是如何把这个公式给化简了，或者就要超时：

有矩阵A，则the sum S = A + A2 + A3 + … + Ak.求矩阵A的幂
 公式可以化简成：若k为奇数=(k&1!=0)  sum(k)=(1+A^(k/2+1))*sum(k/2)+A^(k/2+1);  
若k为偶数： sum(k)=(1+A^(k/2))*sum(k/2)  

代码如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#define ma 10000
using namespace std;
struct mat{
	int x[35][35];
}a,b,c,s;
int num,k1;
mat multiply(mat a,mat b,int n) {
	mat temp;
	int i,j,k;
	memset(temp.x,0,sizeof(temp.x));
	for(i=0; i<n; i++)
		for(j=0; j<n; j++)
			for(k=0; k<n; k++) {
				temp.x[i][j]+=a.x[i][k]*b.x[k][j];
				temp.x[i][j]%=k1;
			}
	return temp;			
}
mat cal(int n) {
	mat c1=c,b1=b;
	while(n) {
		if(n&1) b1=multiply(b1,c1,num);
		n>>=1;
		c1=multiply(c1,c1,num);
	}
	return b1;
}
mat add(mat a,mat b) {
	mat t;
	memset(t.x,0,sizeof(t.x));
	for(int i=0; i<num; i++) {
		for(int j=0; j<num; j++) {
			t.x[i][j]=(a.x[i][j]+b.x[i][j])%k1;
		}
	}
	return t;
}
mat sum(int k) {
	if(k==1) return a;
	mat mul,t;
	mul=sum(k/2);
	//cout<<k<<endl;
	if(k%2) {
		t=cal(k/2+1);
		mat p = add(mul,multiply(mul,t,num));
		t = add(p,t);
	}
	else {
		//cout<<k<<endl;
		t=cal(k/2);
		//cout<<k<<endl;
		t=add(mul,multiply(mul,t,num));
	}
	return t;
}
int main() {
	int m,i,j;
	scanf("%d%d%d",&num,&m,&k1);
	for(i=0; i<num; i++)
		for(j=0; j<num; j++) {
			scanf("%d",&a.x[i][j]);
			a.x[i][j]=a.x[i][j]%k1;
			c.x[i][j]=a.x[i][j];
			if(i==j)
			  b.x[i][j]=1;
			  else b.x[i][j] = 0;
		}
		//for(j=0; j<num; j++) b.x[j][j] = 1;
	s=sum(m);
	for(i=0; i<num; i++)	 {
		printf("%d",s.x[i][0]);
		for(j=1; j<num; j++) 
			printf(" %d",s.x[i][j]);
			printf("\n");
	}	
		return 0;
}