bzoj4162 shlw loves matrix II

传送门

对于一个$k\ast k$的矩阵$A$，它的特征多项式是一个$k$次多项式。
下面我们把多项式的自变量$x$看做一个矩阵。
首先我们假设已经求出了特征多项式$f(x)$
特征多项式满足：$f(A)=0$。(Cayley-Hamilton定理)
我们令$G(x)=x^n$，那么$A^n$即为$G(A)$，也就是要求的东西。
考虑$G(x)$和$f(x)$取个模，得到：
$$G(x)=Q(x)\ast f(x)+R(x)$$
注意到$f(A)=0$，所以当我们把$A$代入之后可以得到：
$$G(A)=Q(A)\ast f(A)+R(A)=R(A)$$
于是我们只要求出$R$后，把矩阵$A$代入$R$，即可计算出结果。
$R$相当于是一个一次多项式$x$在模$f$意义下的$n$次幂（$x^n$）。那么在快速幂的时候不断取模就能算出$R$了。
由于模数并不友好，又$k⩽50$，所以可以暴力大除法$O(k^2)$取模，暴力$O(k^2)$多项式乘法。或者如果有时间写个任意模数NTT？
所以知道$f$后，求出$R$的复杂度就是$O(k^{2}logn)$的。

那么如何求$f$呢？
由于$f$是一个$k$次多项式，所以找到$k+1$个点然后拉格朗日插值即可。
$f(\lambda )=|\lambda I-A|$
那么取$\lambda =0,1...k-1,k$，然后代进去之后可以得到$k+1$个矩阵。
用高斯消元求出它们的行列式，然后你就得到了$k+1$个点。
高斯消元复杂度$O(k^3)$，枚举$k$个点$O(k)$，拉格朗日插值复杂度$O(k^2)$，所以这里复杂度$O(k^4+k^2)$
于是总复杂度$O(k^4+k^2+k^{2}logn)$

#include<bits/stdc++.h>
#define ll long long
#define re register
#define cs const
typedef std::vector<int> poly;
cs int N=1e5+10,mod=1e9+7,K=51;
namespace IO{
	cs int Rlen=1<<22|1;
	char buf[Rlen],*p1,*p2;
	inline char gc(){return (p1==p2)&&(p2=(p1=buf)+fread(buf,1,Rlen,stdin),p1==p2)?EOF:*p1++;}
	template <typename T>
	inline T get(){
		char ch=gc();T x=0,f=1;
		while(!isdigit(ch)){if(ch=='-')f=0;ch=gc();}
		while(isdigit(ch)) x=(x+(x<<2)<<1)+(ch^48),ch=gc();
		return f?x:-x;
	}
	inline int get(int *s){
		memset(s,0,sizeof s);
		char ch=gc();int len=0;
		while(!isdigit(ch)) ch=gc();
		while(isdigit(ch)) s[len++]=ch-'0',ch=gc();
		return len;
	}
	inline int gi(){return get<int>();}
	inline ll gl(){return get<ll>();}
}
using IO::gi;
using IO::gl;
using IO::get;
int x[N],y[N],S[N],k,Len;poly all(1,1),P[N];
inline int max(int x,int y){return x>y?x:y;}
inline int min(int x,int y){return x<y?x:y;}
inline void Max(int &x,int y){if(x<y)x=y;}
inline void Min(int &x,int y){if(x>y)x=y;}
inline int add(int x,int y){return x+y>=mod?x+y-mod:x+y;}
inline int dec(int x,int y){return x-y<0?x-y+mod:x-y;}
inline int mul(int x,int y){return 1ll*x*y%mod;}

inline int Add(int &x,int y){x=(x+y>=mod)?x+y-mod:x+y;}
inline int Dec(int &x,int y){x=(x-y<0)?x-y+mod:x-y;}
inline int Mul(int &x,int y){x=1ll*x*y%mod;}

inline int quickpow(int a,int b,int ret=1){for(;b;b>>=1,Mul(a,a))if(b&1)Mul(ret,a);return ret;}
inline int inv(int x){return quickpow(x,mod-2);}

inline poly operator*(poly A,poly B){
	int len=A.size()+B.size()-2;poly C(len+1,0);
	for(int re i=0;i<A.size();++i)
		for(int re j=0;j<B.size();++j)
			Add(C[i+j],mul(A[i],B[j]));
	return C;
}

inline poly operator/(poly A,poly B){
	int lenA=A.size()-1,lenB=B.size()-1,Inv=inv(B[lenB]);
	if(lenA<lenB) return poly(1,0);poly C(lenA-lenB+1,0);
	for(int re i=lenA;i>=lenB;--i) if(A[i]){
		C[i-lenB]=mul(A[i],Inv);
		for(int re j=i;j>=i-lenB;--j)
			Dec(A[j],mul(C[i-lenB],B[lenB-i+j]));
	}return C;
}

inline poly operator%(poly A,poly B){
	int lenA=A.size()-1,lenB=B.size()-1,Inv=inv(B[lenB]);
	if(lenA<lenB) return A;
	for(int re i=lenA;i>=lenB;--i) if(A[i]){
		int t=mul(A[i],Inv);
		for(int re j=i;j>=i-lenB;--j)
			Dec(A[j],mul(t,B[lenB-i+j]));
	}return A.resize(lenB+1),A;
}

inline poly operator+(poly A,poly B){
	int lenA=A.size()-1,lenB=B.size()-1,lenC=max(lenA,lenB);
	A.resize(lenC+1),B.resize(lenC+1);poly C(lenC+1,0);
	for(int re i=0;i<=lenC;++i) C[i]=add(A[i],B[i]);
	return C;
}

inline poly operator*(poly A,int v){
	for(int re i=0;i<A.size();++i) Mul(A[i],v);
	return A;
}

inline poly quickpow(poly A,int *b,cs poly &M){
	poly ret(1,1);
	for(int re i=Len;i;--i,A=(A*A)%M) if(b[i]) ret=(ret*A)%M;
	return ret;
}

inline poly get_poly(int pos,int *x,int *y,int num,int MUL=1){
	for(int re i=1;i<=num;++i) if(i!=pos)
		Mul(MUL,dec(x[pos],x[i]));
	return ((all/P[pos])*inv(MUL))*y[pos];
}

inline poly poly_fast_inter(int *x,int *y,int num){
	poly ans;
	for(int re i=1;i<=num;++i) P[i]={dec(0,x[i]),1},all=all*P[i];
	for(int re i=1;i<=num;++i) ans=ans+get_poly(i,x,y,num);
	return ans;
}

struct matrix{
	int a[K][K];
	matrix(int t=0){
		memset(a,0,sizeof a);
		for(int re i=1;i<=k;++i) a[i][i]=t;
	}
	inline int det(int ret=1,int flag=1){
		for(int re i=1;i<=k;++i){
			int h=i;
			for(int re j=i+1;j<=k;++j)
				if(a[h][i]<a[j][i]){h=j;break;}
			if(h!=i) std::swap(a[h],a[i]),flag^=1;
			if(!a[i][i]) return 0;
			int Inv=inv(a[i][i]);
			for(int re j=i+1;j<=k;++j){
				for(int re l=i+1;l<=k;++l)
					Dec(a[j][l],mul(mul(a[i][l],a[j][i]),Inv));
			}
		}for(int re i=1;i<=k;++i) Mul(ret,a[i][i]);
		return flag?ret:dec(0,ret);
	}
	friend inline matrix operator+(cs matrix &a,cs matrix &b){
		matrix ret;
		for(int re i=1;i<=k;++i)
			for(int re j=1;j<=k;++j)
				ret.a[i][j]=add(a.a[i][j],b.a[i][j]);
		return ret;
	}
	friend inline matrix operator-(cs matrix &a,cs matrix &b){
		matrix ret;
		for(int re i=1;i<=k;++i)
			for(int re j=1;j<=k;++j)
				ret.a[i][j]=dec(a.a[i][j],b.a[i][j]);
		return ret;
	}
	friend inline matrix operator*(cs matrix &a,int b){
		matrix ret;
		for(int re i=1;i<=k;++i)
			for(int re j=1;j<=k;++j)
				ret.a[i][j]=mul(a.a[i][j],b);
		return ret;
	}
	friend inline matrix operator*(cs matrix &a,cs matrix &b){
		matrix ret;
		for(int re i=1;i<=k;++i)
			for(int re l=1;l<=k;++l)
				for(int re j=1;j<=k;++j)
					Add(ret.a[i][j],mul(a.a[i][l],b.a[l][j]));
		return ret;
	}
	friend inline matrix operator^(matrix a,int b){
		matrix ret(1);
		for(;b;b>>=1,a=a*a) if(b&1) ret=ret*a;
		return ret;
	}
	inline void print(){
		for(int re i=1;i<=k;++i){
			for(int re j=1;j<=k;++j)
				printf("%d ",a[i][j]);
			puts("");
		}puts("");
	}
}A,E;

inline matrix evaluation(cs poly &P,cs matrix &M){
	matrix ret;
	for(int re i=0;i<P.size();++i)
		ret=ret+(M^i)*P[i];
	return ret;
}
int main(){
//	freopen("4162.in","r",stdin);
//	freopen("4162.out","w",stdout);
	Len=get(S+1),k=gi(),E=matrix(1);
	for(int re i=1;i<=k;++i)
		for(int re j=1;j<=k;++j)
			A.a[i][j]=gi();
	for(int re lmb=0;lmb<=k;++lmb)
		x[lmb+1]=lmb,y[lmb+1]=(E*lmb-A).det();
	poly charac_poly=poly_fast_inter(x,y,k+1),ans={0,1};
	ans=quickpow(ans,S,charac_poly),evaluation(ans,A).print();
}