POJ 3233 Matrix Power Series(矩阵快速幂+二分求和)

#include <cstdio>
using namespace std;
const int maxn=100+5;
typedef struct
{
int m[maxn][maxn];
}Matrix;
Matrix a,per;
int n,M;

void init()
{
for(int i=0;i<n;i++)
for(int j=0;j<n;j++){
scanf("%d",&a.m[i][j]);
a.m[i][j]%=M;
per.m[i][j]=(i==j);
}
}
Matrix add(Matrix a,Matrix b)
{
Matrix z;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
z.m[i][j]=(a.m[i][j]+b.m[i][j])%M;
return z;
}
Matrix Multi(Matrix a,Matrix b)
{
Matrix z;
for(int i=0;i<n;i++)
for(int j=0;j<n;j++){
z.m[i][j]=0;
for(int k=0;k<n;k++)
z.m[i][j]+=a.m[i][k]*b.m[k][j];
z.m[i][j]%=M;
}
return z;
}
Matrix power(int k)
{
Matrix r=per,p=a;
while(k)
{
if(k&1) r=Multi(r,p);
p=Multi(p,p);
k>>=1;
}
return r;
}
Matrix MatrixSum(int k)///二分求和(等比数列的形式)
{
if(k==1) return a;
Matrix tmp,b;
tmp=MatrixSum(k/2);
if(k&1)
{
b=power(k/2+1);
tmp=add(tmp,Multi(tmp,b));
tmp=add(tmp,b);
}
else
{
b=power(k/2);
tmp=add(tmp,Multi(tmp,b));
}
return tmp;
}
void output(Matrix a)
{
for(int i=0;i<n;i++)
for(int j=0;j<n;j++)
printf("%d%c",a.m[i][j],j==n-1?'\n':' ');
}
int main()
{
int k;
while(~scanf("%d%d%d",&n,&k,&M))
{
init();
Matrix ans=MatrixSum(k);
output(ans);
}
return 0;
}