POJ 3233 Matrix Power Series (矩阵乘法+快速幂+等比二分求和) -

#include<iostream>
#include<cstring>
#include<string>
#include<cstdio>
#include<algorithm>
using namespace std;
struct Matrix{
int a[32][32];
int r; // r*r
Matrix(int r):r(r){memset(a,0,sizeof(a));} //初始化长
void MakeI(int x){             //变为x*E的单位矩阵
memset(a,0,sizeof(a));
for(int i=0;i<r;i++) a[i][i]=x;
}
Matrix operator * (const Matrix& m) {  //矩阵乘法
Matrix result(r);
for(int i=0;i<r;i++)
for(int j=0;j<r;j++)
{
for(int k=0;k<r;k++)
result.a[i][j]+=a[i][k]*m.a[k][j];
}
return result;
}
Matrix operator % (const int p){ //矩阵对每个元素取余
Matrix result(*this);
for(int i=0;i<r;i++)
for(int j=0;j<r;j++)
result.a[i][j]%=p;
return result;
}
Matrix operator + (const Matrix& m){  //矩阵互相相加
Matrix result(r);
for(int i=0;i<r;i++)
for(int j=0;j<r;j++)
result.a[i][j]=a[i][j]+m.a[i][j];
return result;
}
Matrix operator + (int n){     //加上个常数
Matrix result(r);
result.MakeI(n);
return *this+result;
}
};
Matrix PowMod(const Matrix& m,int k,int p)  //矩阵快速幂
{   //m^k mod p
int r=m.r;
Matrix result(r);
result.MakeI(1); //变成单位矩阵
Matrix base=m;
while(k){
if(k&1) result=result*base%p;
base=base*base%p;
k>>=1;
}
return result;
}
Matrix PowSumMod(const Matrix& m,int k,int p)   //等比二分求和
{
if(k==1) return Matrix(m)%p;
else if(k%2==0) return (PowMod(m,k/2,p)+1)*PowSumMod(m,k/2,p)%p;
else return ((PowMod(m,(k-1)/2,p)+1)*PowSumMod(m,(k-1)/2,p)%p + PowMod(m,k,p))%p;
}
int main()
{
int r,k,p;
cin>>r>>k>>p;
Matrix m(r);
for(int i=0;i<r;i++)
for(int j=0;j<r;j++)
cin>>m.a[i][j],m.a[i][j]%=p;  //一定这里就要取模，不然会WR!!
Matrix result=PowSumMod(m,k,p);
for(int i=0;i<r;i++)
{
for(int j=0;j<r;j++)
{
if(j) cout<<' ';
cout<<result.a[i][j];
}
cout<<endl;
}
return 0;
}