bzoj 4162 shlw loves matrix II

拉格朗日插值求特征多项式 + 多项式取模.

$k\times k$ 的方阵 $M$ 的特征多项式 $f(x)=|xI-M|$ ,它是个 $k$ 次多项式,首项系数为 $1$ .

可以随便选 $k+1$ 个 $x$ ,代到多项式里,用高斯消元算出对应的 $f(x)​$ .

然后利用拉格朗日插值得到特征多项式 $f(x)$ 的系数.

根据 $Cayley-Hamilton$ 定理, $f(M)=0$ ,

那么我们给 $M^n$ 加上任意多个 $f(M)$ 都不会影响答案,即 $M^n=M^n\bmod f(M)$ .

只需计算 $M^n\bmod f(M)$ ,两者都是关于 $M$ 的多项式,模数不是 $NTT$ 模数,多项式快速幂 + 暴力取模.

最终得到的 $M^n\bmod f(M)$ 次数显然小于 $k$ ,将 $M$ 代进去暴力算即可.

时间复杂度 $O(k^4+k^2\log n)$ .

Copy

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
inline int read()
{
	int out=0,fh=1;
	char jp=getchar();
	while ((jp>'9'||jp<'0')&&jp!='-')
		jp=getchar();
	if (jp=='-')
		fh=-1,jp=getchar();
	while (jp>='0'&&jp<='9')
		out=out*10+jp-'0',jp=getchar();
	return out*fh;
}
const int P=1e9+7;
int add(int a,int b)
{
	return (a+b>=P)?(a+b-P):(a+b);
}
int mul(int a,int b)
{
	return 1LL * a * b % P;
}
int fpow(int a,int b)
{
	int res=1;
	while(b)
	{
		if(b&1)
			res=mul(res,a);
		a=mul(a,a);
		b>>=1;
	}
	return res;
}
const int MAXN=51,L=1e4+10;
char buf[L];
int n,m,len,tmp[MAXN<<1],p[MAXN];
void Mul(int *a,int *b,int *f)
{
	memset(tmp,0,sizeof tmp);
	for(int i=0;i<=n;++i)
		for(int j=0;j<=n;++j)
			tmp[i+j]=add(tmp[i+j],mul(a[i],b[j]));
	for(int i=m;i>=n;--i)
	{
		int inv=mul(tmp[i],fpow(p[n],P-2));
		for(int j=0;j<=n;++j)
			tmp[i-j]=add(tmp[i-j],P-mul(p[n-j],inv));
	}		
	for(int i=0;i<=n;++i)
		f[i]=tmp[i];
}
int a[MAXN][MAXN],b[MAXN][MAXN];
int Det()
{
	int res=1,sgn=1;
	for(int i=1;i<=n;++i)
	{
		for(int j=i;j<=n;++j)
			if(b[j][i])
			{
				if(i!=j)
					swap(b[i],b[j]),sgn*=-1;
				break;
			}
		int inv=fpow(b[i][i],P-2);
		for(int j=i+1;j<=n;++j)
		{
			if(!b[j][i])
				continue;
			int t=mul(b[j][i],inv);
			for(int k=i;k<=n;++k)
				b[j][k]=add(b[j][k],P-mul(b[i][k],t));
		}
		res=mul(res,b[i][i]);
	}
	if(sgn==1)
		return res;
	else
		return add(P-res,0);
}
int y[MAXN],c[MAXN],d[MAXN],mat[MAXN][MAXN],rmat[MAXN][MAXN];
int main()
{
	scanf("%s%d",buf+1,&n);
	len=strlen(buf+1);
	m=n<<1;
	for(int i=1;i<=n;++i)
		for(int j=1;j<=n;++j)
			a[i][j]=read();
	for(int i=0;i<=n;++i)
	{
		memset(b,0,sizeof b);
		for(int j=1;j<=n;++j)
			for(int k=1;k<=n;++k)
				b[j][k]=(j==k)?(add(i,P-a[j][k])):(add(0,P-a[j][k]));
		y[i]=Det();
	}
	for(int i=0;i<=n;++i)
	{
		memset(tmp,0,sizeof tmp);
		tmp[0]=y[i];
		for(int j=0;j<=n;++j)
			if(i!=j)
			{
				for(int k=n;k;--k)
					tmp[k]=add(tmp[k-1],P-mul(tmp[k],j));
				tmp[0]=add(0,P-mul(tmp[0],j));
				int inv=fpow(add(i,P-j),P-2);
				for(int k=0;k<=n;++k)
					tmp[k]=mul(tmp[k],inv);
			}
		for(int j=0;j<=n;++j)
			p[j]=add(p[j],tmp[j]);
	}
	c[0]=d[1]=1;
	for(int i=len;i;--i)
	{
		if(buf[i]=='1')
			Mul(c,d,c);
		Mul(d,d,d);
	}
	memset(b,0,sizeof b);
	for(int i=1;i<=n;++i)
		mat[i][i]=1;
	for(int l=0;l<=n;++l)
	{
		for(int i=1;i<=n;++i)
			for(int j=1;j<=n;++j)
				b[i][j]=add(b[i][j],mul(c[l],mat[i][j]));
		memset(rmat,0,sizeof rmat);
		for(int i=1;i<=n;++i)
			for(int k=1;k<=n;++k)
				for(int j=1;j<=n;++j)	
					rmat[i][j]=add(rmat[i][j],mul(mat[i][k],a[k][j]));
		memcpy(mat,rmat,sizeof mat);
	}
	for(int i=1;i<=n;++i)
	{
		for(int j=1;j<=n;++j)
			printf("%d ",b[i][j]);
		puts("");
	}
	return 0;
}