再加上快速冪演算法和就好了

#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;
}