poj 3233 Matrix Power Series 二分求和式

##代码

#include<cstdio>
#include<cstdlib>
#include<cstring>
struct Mat{
	int ma[35][35];
	Mat unit(int n){
		Mat res;
		memset(res.ma,0,sizeof res.ma);
		for(int i=0;i<n;++i)
			res.ma[i][i]=1;	
		return res;
	}		       
};
int n,m,k;
Mat add(Mat a,Mat b)  
{  
    Mat res;  
    for(int i=0; i<n; i++)    
        for(int j=0; j<n; j++)    
            res.ma[i][j] =(a.ma[i][j] + b.ma[i][j])%m;       
    return res;  
}  
Mat multi(Mat a,Mat b){
	Mat res;
	for(int i=0; i<n; i++)  
    {  
        for(int j=0; j<n; j++)  
        {  
            res.ma[i][j] = 0;  
            for(int k=0; k<n; k++)  
                res.ma[i][j] =(res.ma[i][j]+a.ma[i][k] * b.ma[k][j]%m)%m;  
        }  
    }  
	return res;	
}
Mat power(Mat a,int b){
	Mat ans;
	ans=ans.unit(n);  
    while(b)  
    {  
        if(b & 1)    
            ans = multi(ans,a);
		a = multi(a,a);	           
        b >>= 1;  
    }  
    return ans;
}
Mat solve(Mat a,int b){
	if(b==1) return a;
	Mat t=solve(a,b/2);
	if(b&1){
		Mat cur=power(a,b/2+1);
		t=add(t,multi(t,cur));
		t=add(t,cur);
	}
	else{
		Mat cur=power(a,b/2);
		t=add(t,multi(t,cur));	
	}	
	return t;
} 
int main(){    
	Mat A,ans;;
    scanf("%d%d%d",&n,&k,&m);
    for(int i=0;i<n;++i){
		for(int j=0;j<n;++j){
			scanf("%d",&A.ma[i][j]);	
		}	
	}
	ans=solve(A,k);
	for(int i=0;i<n;++i){
		for(int j=0;j<n;++j){
			printf("%d%c",ans.ma[i][j],j==n-1?'\n':' ');	
		}	
	}
    return 0;
}