[Poj3233]Matrix Power Series

题意

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

运算结果对$m$取模

$n\le 30,m<10^4,k\le 10^9$

---

题解

设$E$为单位矩阵

方法一 分治

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

那么

$$f(k)=\{\begin{array}{ll}f(\frac{k}{2})(E+A^{\frac{k}{2}})& \text{2∣k}\\ f(\frac{k-1}{2})(E+A^{\frac{k-1}{2}})+A^k& \text{2∤k}\end{array}$$

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

时间复杂度$O(n^3{\log }^{2}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+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=\left[\begin{array}{cc}A& E\\ 0& E\end{array}\right]$$
可以发现
$$S^2=\left[\begin{array}{cc}{A}^2& A+E\\ 0& E\end{array}\right]S^3=\left[\begin{array}{cc}{A}^3& A^2+A+E\\ 0& E\end{array}\right]...S^{k+1}=\left[\begin{array}{cc}{A}^{k+1}& A^k+A^{k-1}+...+A+E\\ 0& E\end{array}\right]$$
然后根据分块矩阵理论可以把$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;
}