数论 快速矩阵幂 POJ 3233 Matrix Power Series 二分和

题意：给出矩阵A，求S = A + A^2 + A^3 + … + A^k ，n为矩阵

题解：

矩阵快速幂  写出矩阵 * 与+的功能

二分求和

#include <cstdio>
#include <cstring>
#include <iostream>
using namespace std;
const int N=33;
int n,mod;
struct mat{
	//矩阵 
	int n,at[N][N];
	mat(int n=0){
		this->n=n;
		memset(at,0,sizeof(at));
	}
	mat operator *(mat b){
		mat tmp=mat(n);
		for(int i=0;i<n;i++){
			for(int k=0;k<n;k++){
				if(at[i][k]){
					for(int j=0;j<n;j++){
						tmp.at[i][j]=(tmp.at[i][j]+(at[i][k] * b.at[k][j])%mod)%mod;
					}
				}
			}
		}
		return tmp;
	}
	mat operator +(mat b){
		mat tmp=mat(n);
		for(int i=0;i<n;i++){
			for(int j=0;j<n;j++){
				tmp.at[i][j]=(at[i][j]+b.at[i][j])%mod;
			}
		}
		return tmp;
	} 
};
mat expo(mat d,int k){
	if(k==1) return d;
	mat e(n);
	memset(e.at,0,sizeof(e.at));
	for(int i=0;i<n;i++) 
		e.at[i][i]=1;
	if(k==0) return e;
	while(k)
	{
		if(k&1) e=d*e;
		d=d*d;
		k>>=1;
	}
	return e;
}
mat sum(mat d,int k){
	if(k==1) return d;
	if(k&1){//奇数  
//	printf("=%d\n",k);	
		// d^k+上一步累计下来的和 
		return expo(d,k)+sum(d,k-1); 
	}
	else {  //偶数 
//	printf("-%d\n",k);
		mat s=sum(d,k>>1);
		// d^(k/2)*（上一步累计的和）+上一步累计的和
		// 举例： 
		// d^2 * (d^1 + d^2) + (d^1 + d^2)
		//        d^3 + d^4  +  d^1 + d^2
		return expo(d,k>>1)*s+s;
	}
}
int main()
{
	int k;
	while(~scanf("%d%d%d",&n,&k,&mod)){
		mat a(n);
		for(int i=0;i<n;i++)
			for(int j=0;j<n;j++){
				scanf("%d",&a.at[i][j]);
				a.at[i][j]=a.at[i][j]%mod;
			}
		mat ans=sum(a,k);
		for(int i=0;i<n;i++){
			for(int j=0;j<n-1;j++){
				printf("%d ",ans.at[i][j]);
			}
			printf("%d\n",ans.at[i][n-1]);
		}
			
	}
	return 0;
}