BZOJ 4162 Shlw loves matrix II

最新推荐文章于 2021-07-19 16:35:27 发布

Defener

最新推荐文章于 2021-07-19 16:35:27 发布

阅读量192

点赞数

分类专栏：数学文章标签：矩阵

本文链接：https://blog.csdn.net/Defener/article/details/111464828

版权

数学专栏收录该内容

15 篇文章 1 订阅

订阅专栏

题面
此处为计算矩阵幂，直接快速幂不够快，考虑用数学方法优化。
对于一个 $\times k$ 的方阵 $M$ ，都有一个特征多项式 $f (x)$ 满足 $f (M) = 0$ 。 $f (x)$ 的构造方法如下：
$f (x) = d e t (M - I x)$
其中， $I$ 为单位矩阵，显然 $f (M) = 0$ 成立。
对于这个多项式，直接求行列式不好求，但可以知道的是 $f (x)$ 为 $k$ 次的，因此利用拉格朗日插值，求得 $f(x)=\sum_{i=0}^k a_i M^i$ ，其中 $a_i$ 为已知常数。因此，可以得到：
$M^{k+c}=\sum_{i=0}^{k-1}-\frac{a_i}{a_k}M^{i+c}$
利用这个方程降次。对于 $M$ 的任意次幂可以用一个 $k - 1$ 次多项式表达出来。这时，模拟之前的矩阵快速幂，改为多项式快速幂。设 $M^a=g_1(M),M^b=g_2(M)$ ，则 $M^{a+b}=g_1(M)g_2(M)$ 。这时 $g_1(x)g_2(x)$ 结果为一个 $2 k - 2$ 次多项式，将次幂从高到低依次降次，替换为更低次幂的多项式。设最后将答案表达为 $h (x)$ ，则最后再还原，求出 $h (M)$ 即可。
时间复杂度 $O(k^4+log_2n k^2)$ ，空间复杂度 $O(k^2)$

#include<stdio.h>
#define R register int
#define L long long
#define I inline
#define P 1000000007
char s[10001];
I void Swap(int&x,int&y){
	int tem=x;
	x=y;
	y=tem;
}
I int Read(){
	int len=-1;
	do{
		len++;
		s[len]=getchar();
	}while(s[len]!=' ');
	return len;
}
I int PowMod(int x,int y){
	int res=1;
	while(y!=0){
		if((y&1)==1){
			res=(L)res*x%P;
		}
		x=(L)x*x%P;
		y>>=1;
	}
	return res;
}
typedef int Matrix[50][50],Poly[50];
I int Calc(Matrix B,int&n,const int x){
	Matrix A;
	for(R i=0;i!=n;i++){
		for(R j=0;j!=n;j++){
			A[i][j]=B[i][j];
		}
		A[i][i]+=P-x;
		if(A[i][i]>=P){
			A[i][i]-=P;
		}
	}
	int res=1;
	for(R i=0;i!=n;i++){
		if(A[i][i]==0){
			int p=-1;
			for(R j=i+1;j!=n;j++){
				if(A[j][i]!=0){
					p=j;
					break;
				}
			}
			if(p==-1){
				return 0;
			}
			res=P-res;
			for(R j=i;j!=n;j++){
				Swap(A[i][j],A[p][j]);
			}
		}
		res=(L)res*A[i][i]%P;
		int tem=PowMod(A[i][i],P-2);
		for(R j=i+1;j!=n;j++){
			int tem2=(L)tem*A[j][i]%P;
			for(R k=i;k!=n;k++){
				A[j][k]-=(L)tem2*A[i][k]%P;
				if(A[j][k]<0){
					A[j][k]+=P;
				}
			}
		}
	}
	return res;
}
int a[51],f[51],g[51];
I void Init(int&n){
	g[0]=1;
	f[0]=a[0]; 
	for(R i=1;i<=n;i++){
		int p=1,c=0;
		for(R j=i-1;j!=-1;j--){
			c=((L)c*i+f[j])%P;
			p=(L)p*(i-j)%P;
		}
		p=(L)(a[i]-c)*PowMod(p,P-2)%P;
		if(p<0){
			p+=P;
		}
		for(R j=1;j<=i;j++){
			f[j]=((L)p*g[j-1]+f[j])%P;
		}
		for(R j=i;j!=0;j--){
			g[j]=((L)g[j]*(P-i)+g[j-1])%P;
		}
		g[0]=(L)g[0]*(P-i)%P;
	}
}
Poly T;
I void PolyMul(Poly A,Poly B,int&n){
	static int t[100];
	for(R i=0;i!=n;i++){
		for(R j=0;j!=n;j++){
			t[i+j]=((L)A[i]*B[j]+t[i+j])%P;
		}
	}
	for(R i=n-1<<1;i>=n;i--){
		for(R j=0;j!=n;j++){
			t[i-n+j]=((L)t[i]*T[j]+t[i-n+j])%P;
		}
		t[i]=0;
	}
	for(R i=0;i!=n;i++){
		A[i]=t[i];
		t[i]=0;
	}
}
I void MatrixAdd(Matrix A,Matrix B,int x,int&n){
	for(R i=0;i!=n;i++){
		for(R j=0;j!=n;j++){
			A[i][j]=((L)x*B[i][j]+A[i][j])%P;
		}
	}
}
I void MatrixMul(Matrix A,Matrix B,int&n){
	Matrix C;
	for(R i=0;i!=n;i++){
		for(R j=0;j!=n;j++){
			C[i][j]=0;
		}
	}
	for(R i=0;i!=n;i++){
		for(R j=0;j!=n;j++){
			register unsigned L t=A[i][j];
			for(R k=0;k!=n;k++){
				C[i][k]=(t*B[j][k]+C[i][k])%P;
			}
		}
	}
	for(R i=0;i!=n;i++){
		for(R j=0;j!=n;j++){
			A[i][j]=C[i][j];
		}
	}
}
int main(){
	int l=Read(),n;
	scanf("%d",&n);
	Matrix A,S,B;
	for(R i=0;i!=n;i++){
		for(R j=0;j!=n;j++){
			scanf("%d",A[i]+j);
			S[i][j]=B[i][j]=0;
		}
		B[i][i]=1;
	}
	for(R i=0;i<=n;i++){
		a[i]=Calc(A,n,i);
	}
	Init(n);
	Poly D,F;
	int tem=PowMod(P-f[n],P-2);
	for(R i=0;i!=n;i++){
		D[i]=F[i]=0;
		T[i]=(L)f[i]*tem%P;
	}
	D[1]=F[0]=1;
	for(R i=l-1;i!=-1;i--){
		if(s[i]=='1'){
			PolyMul(F,D,n);
		}
		PolyMul(D,D,n);
	}
	for(R i=0;i!=n;i++){
		MatrixAdd(S,B,F[i],n);
		MatrixMul(B,A,n);
	}
	for(R i=0;i!=n;i++){
		for(R j=0;j!=n;j++){
			printf("%d ",S[i][j]);
		}
		puts("");
	}
	return 0;
}