题意

给你一个 $n\times n$的矩阵 $A$，求 $A+A^2+...+A^k$

运算结果对 $m$取模

$n\le30,m\lt10^4,k\le10^9$

题解

设 $E$为单位矩阵

方法一 分治

令 $f(k)=A+A^2+...+A^k$

那么

$f(k)= \begin{cases} f(\frac k2)(E+A^{\frac k2})& \text{$2\mid k$} \\ f(\frac{k-1}2)(E+A^{\frac{k-1}2})+A^k& \text{$2\nmid k $} \end{cases}$

分治计算 $f(k)$即可

时间复杂度 $O(n^3\log^2k)$

#include<cstdio>
#include<cstring>
#define fp(i,a,b) for(register int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(register int i=a,I=b-1;i>I;--i)
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
template<class T>inline bool cmax(T&a,const T&b){return a<b?a=b,1:0;}
template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;}
using namespace std;
const int N=30+5,inf=~0u>>1;
typedef int arr[N];
int n,Pow,Mod;
inline int plus(int a,int b){return (a+=b)>=Mod?a-Mod:a;}
struct Matrix{
	int a[N][N];
	Matrix(int d=0){memset(a,0,sizeof a);if(d)fp(i,1,N-1)a[i][i]=1;}
	inline int*operator[](int x){return a[x];}
	Matrix operator+(Matrix b){
		Matrix c;
		fp(i,1,n)fp(j,1,n)
			c[i][j]=plus(a[i][j],b[i][j]);
		return c;
	}
	inline Matrix operator*(Matrix b)const{
		Matrix c;
		fp(k,1,n)fp(i,1,n)if(a[i][k])fp(j,1,n)
			c[i][j]=plus(c[i][j],a[i][k]*b[k][j]%Mod);
		return c;
	}
	Matrix operator^(int b){
		Matrix x(1),A=*this;
		for(;b;b>>=1,A=A*A)if(b&1)x=x*A;
		return x;
	}
}A,S,E(1),Z;
Matrix f(int k){
	if(k==1)return A;
	return f(k>>1)*(E+(A^(k>>1)))+(k&1?A^k:Z);
}
int main(){
	#ifndef ONLINE_JUDGE
		file("s");
	#endif
	scanf("%d%d%d",&n,&Pow,&Mod);
	fp(i,1,n)fp(j,1,n)
		scanf("%d",&A[i][j]);
	S=f(Pow);
	fp(i,1,n){
		fp(j,1,n)printf("%d ",S[i][j]);
		puts("");
	}
return 0;
}

方法二 构造

令
$S^2= \begin{bmatrix} A&E\\ 0&E \end{bmatrix}$
可以发现
$S^2= \begin{bmatrix} A^2&A+E\\ 0&E \end{bmatrix} \ \ \ S^3= \begin{bmatrix} A^3&A^2+A+E\\ 0&E \end{bmatrix}\ ...\ \ S^{k+1}= \begin{bmatrix} A^{k+1}&A^k+A^{k-1}+...+A+E\\ 0&E \end{bmatrix}$
然后根据分块矩阵理论可以把 $S$扩充成一个 $2n\times 2n$的新矩阵做矩阵快速幂

记得最后减去单位矩阵

时间复杂度 $O(n^3\log k)$

#include<cstdio>
#include<cstring>
#define fp(i,a,b) for(register int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(register int i=a,I=b-1;i>I;--i)
#define file(s) freopen(s".in","r",stdin),freopen(s".out","w",stdout)
template<class T>inline bool cmax(T&a,const T&b){return a<b?a=b,1:0;}
template<class T>inline bool cmin(T&a,const T&b){return a>b?a=b,1:0;}
using namespace std;
const int N=30*2+5,inf=~0u>>1;
typedef int arr[N];
int n,Pow,Mod;
struct Matrix{
	int a[N][N];
	Matrix(int d=0){memset(a,0,sizeof a);if(d)fp(i,1,n)a[i][i]=1;}
	inline int*operator[](int x){return a[x];}
	inline Matrix operator*(Matrix b)const{
		Matrix c;
		fp(k,1,n)fp(i,1,n)if(a[i][k])fp(j,1,n)
			c[i][j]=(c[i][j]+a[i][k]*b[k][j])%Mod;
		return c;
	}
	Matrix operator^(int b){
		Matrix x(1),A=*this;
		for(;b;b>>=1,A=A*A)if(b&1)x=x*A;
		return x;
	}
}S;
int main(){
	#ifndef ONLINE_JUDGE
		file("s");
	#endif
	scanf("%d%d%d",&n,&Pow,&Mod);
	fp(i,1,n)fp(j,1,n){
		scanf("%d",&S[i][j]);
		if(i==j)S[i+n][j+n]=S[i][j+n]=1;
	}
	n<<=1;
	S=S^(Pow+1);
	n>>=1;
	fp(i,1,n){
		fp(j,1,n)printf("%d ",(S[i][j+n]-(i==j)+Mod)%Mod);
		puts("");
	}
return 0;
}