转载自大神脑钊文:http://blog.csdn.net/wzw1376124061/article/details/52492665
#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
int n,k;
int mod; //n*n阶的矩阵;k是sum=A^1+A^2+.....+A^k;m是被模的元素
struct Matrix{
int matrix[50][50];
Matrix(int a=0){
memset(matrix,0,sizeof(matrix)); //矩阵清零
if(a==1) //单位矩阵
for(int i=0;i<50;i++)
matrix[i][i]=1;
}
}m;
Matrix operator * (Matrix a,Matrix b){ //矩阵乘法
Matrix res;
memset(res.matrix,0,sizeof(res.matrix));
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++)
res.matrix[i][j]=(res.matrix[i][j]+(a.matrix[i][k]*b.matrix[k][j])%mod)%mod;
return res;
}
Matrix operator + (Matrix a,Matrix b){ //矩阵加法
Matrix res;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
res.matrix[i][j]=(a.matrix[i][j]+b.matrix[i][j])%mod;
return res;
}
Matrix operator ^ (Matrix a,int k){ //矩阵快速幂
bool flag=false;
Matrix ans=1;
while(k){
if(k&1){
if(flag)ans=ans*a;
else ans=a;
flag=true;
}
a=a*a;
k>>=1;
}
return ans;
}
Matrix sum(int k){ //求sum( S(k) )= A + A2 + A3 + …+ Ak
if(k==1)return m; //递归回去
else{
Matrix tmp=sum(k>>1); //sum( S(k/2) )
if(k&1){ //k为奇数
Matrix tmp2=m^((k>>1)+1);
return tmp+tmp2+tmp*tmp2;
}
else{
Matrix tmp2=m^(k>>1);
return tmp+tmp*tmp2; //sum( S(k/2) )*(tmp2+1)
}
}
}
int main(){
int i,j;
Matrix ans;
scanf("%d%d%d",&n,&k,&mod);
for(i=1;i<=n;i++)
for(j=1;j<=n;j++){
scanf("%d",&m.matrix[i][j]);
m.matrix[i][j]%=mod;
}
ans=sum(k);
for(i=1;i<=n;i++,puts(""))
for(j=1;j<=n;j++)
printf("%d ",ans.matrix[i][j]);
return 0;
}