A为一个n*n的矩阵，求A+A^2+A^3+...+A^n
S = A + A + A + … + A
=(1+A)*(A + A + A +
… + A)+{A}
=(1+A)*(S)+{A}//
k为偶数时无 {A}
A
可用二分迭代求出
因此,只要求出 上面的三部分就可以求出 S
设f(n)=A+A^2+A^3+...+A^n
n%2==1时，f(n)=f(n-1)+A^n
n%2==0时，f(n)=f(n/2)+f(n/2)*A^(n/2)
由于矩阵乘法满足结合律，计算A^n时，也可以二分
#include<iostream>
#include<cstdio>
using namespace std;
int n,m;
struct Mat
{
int mat[31][31];
Mat operator*(const Mat &x)
{
Mat tmp;
int i,j,k;
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
{
tmp.mat[i][j]=0;
for(k=0;k<n;k++)
{
tmp.mat[i][j]+=(mat[i][k]*x.mat[k][j])%m;
tmp.mat[i][j]%=m;
} }
}
return tmp;
}
Mat operator+(const Mat &x)
{
Mat tmp;
int i,j;
for(i=0;i<n;i++)
for(j=0;j<n;j++)
{
tmp.mat[i][j]=(mat[i][j]+x.mat[i][j])%m;
}
return tmp;
}
}p;
Mat _pow(int k)
{
if(k==1)
return p;
if(k&1)
return _pow(k-1)*p;
else
{
Mat tmp=_pow(k/2); return tmp*tmp;
}
}
Mat cal(int k)
{
if(k==1)
return p;
else
{
if(k&1)
return cal(k-1)+_pow(k);
else
{
Mat tmp=cal(k/2);
return tmp+tmp*_pow(k/2);
}
}
}
int main()
{
int i,j,k;
scanf("%d%d%d",&n,&k,&m);
for(i=0;i<n;i++)
{
for(j=0;j<n;j++)
scanf("%d",&p.mat[i][j]);
}
Mat res=cal(k);
for(i=0;i<n;i++)
{
for(j=0;j<n-1;j++)
{
printf("%d ",res.mat[i][j]);
}
printf("%d\n",res.mat[i][j]);
}
return 0;
}