【BZOJ4162】shlw loves matrix II（特征多项式）

一般看到这种求某个矩阵的多项式的题就有可能是利用其特征多项式。

给定矩阵 \(M\)，求 \(M^n\)。

求出 \(M\) 的特征多项式 \(f(x)\)，那么 \(f(M)=0\)。

所以我们可以让 \(M^n\) 一直减 \(f(M)\) 直到次数低于 \(f(M)\) 为止。

意思就是我们先求出 \(g(x)=x^n\bmod f(x)\)，然后再代入 \(x=M\)。

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

#include<bits/stdc++.h>

#define N 55

using namespace std;

namespace modular
{
	const int mod=1000000007;
	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;}
}using namespace modular;

inline int poww(int a,int b)
{
	int ans=1;
	while(b)
	{
		if(b&1) ans=mul(ans,a);
		a=mul(a,a);
		b>>=1;
	}
	return ans;
}

inline int read()
{
	int x=0,f=1;
	char ch=getchar();
	while(ch<'0'||ch>'9')
	{
		if(ch=='-') f=-1;
		ch=getchar();
	}
	while(ch>='0'&&ch<='9')
	{
		x=(x<<1)+(x<<3)+(ch^'0');
		ch=getchar();
	}
	return x*f;
}

typedef vector<int> poly;

poly add(poly &a,poly &b,int d)
{
	poly c;
	int sa=a.size(),sb=b.size();
	for(int i=0;i<min(sa,sb);i++)
		c.push_back(add(a[i],mul(d,b[i])));
	if(sa>sb)
	{
		for(int i=sb;i<sa;i++)
			c.push_back(a[i]);
	}
	else
	{
		for(int i=sa;i<sb;i++)
			c.push_back(mul(d,b[i]));
	}
	return c;
}

poly mul(poly &a,poly &b)
{
	int sa=a.size(),sb=b.size();
	poly c(sa+sb-1);
	for(int i=0;i<sa;i++)
		for(int j=0;j<sb;j++)
			c[i+j]=add(c[i+j],mul(a[i],b[j]));
	return c;
}

int n,a[N][N];
int invFn;
char s[10010];
poly f[N],F,G;

struct Matrix
{
	int a[N][N];
	Matrix(){memset(a,0,sizeof(a));}
	void unit(){for(int i=1;i<=n;i++) a[i][i]=1;}
}M;

Matrix mul(Matrix &a,Matrix &b)
{
	Matrix c;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			for(int k=1;k<=n;k++)
				c.a[i][j]=add(c.a[i][j],mul(a.a[i][k],b.a[k][j]));
	return c;
}

Matrix mul(Matrix &a,int b)
{
	Matrix c;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			c.a[i][j]=mul(a.a[i][j],b);
	return c;
}

Matrix add(Matrix &a,Matrix b)
{
	Matrix c;
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			c.a[i][j]=add(a.a[i][j],b.a[i][j]);
	return c;
}

void Gauss()
{
	for(int i=1;i<n;i++)
	{
		int p=i+1;
		for(int j=i+1;j<=n;j++)
			if(a[j][i]>a[p][i]) p=j;
		if(p!=i+1)
		{
			for(int j=1;j<=n;j++)
				swap(a[i+1][j],a[p][j]);
			for(int j=1;j<=n;j++)
				swap(a[j][i+1],a[j][p]);
		}
		for(int j=i+2;j<=n;j++)
		{
			int div=mul(a[j][i],poww(a[i+1][i],mod-2));
			for(int k=1;k<=n;k++)
				a[j][k]=dec(a[j][k],mul(a[i+1][k],div));
			for(int k=1;k<=n;k++)
				a[k][i+1]=add(a[k][i+1],mul(a[k][j],div));
		}
	}
}

void getpoly()
{
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			a[i][j]=dec(0,a[i][j]);
	f[n+1].push_back(1);
	for(int i=n;i>=1;i--)
	{
		poly tmp;
		tmp.push_back(a[i][i]);
		tmp.push_back(1);
		f[i]=mul(tmp,f[i+1]);
		int now=a[i+1][i],c=dec(0,1);
		for(int j=i+1;j<=n;j++)
		{
			f[i]=add(f[i],f[j+1],mul(mul(c,now),a[i][j]));
			now=mul(now,a[j+1][j]);
			c=dec(0,c);
		}
	}
}

poly polymul(poly a,poly b)
{
	int sa=a.size(),sb=b.size();
	poly c(sa+sb-1);
	for(int i=0;i<sa;i++)
		for(int j=0;j<sb;j++)
			c[i+j]=add(c[i+j],mul(a[i],b[j]));
	for(int i=c.size()-1;i>=n;i--)
	{
		int div=mul(c[i],invFn);
		for(int j=i,l=n;l>=0;j--,l--)
			c[j]=dec(c[j],mul(F[l],div));
	}
	c.resize(n);
	return c;
}

void Mod()
{
	poly a;
	a.push_back(0),a.push_back(1);
	G.push_back(1);
	for(int len=strlen(s+1),i=len;i>=1;i--)
	{
		if(s[i]=='1') G=polymul(G,a);
		a=polymul(a,a);
	}
}

int main()
{
	scanf("%s",s+1);
	n=read();
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			a[i][j]=read();
	for(int i=1;i<=n;i++)
		for(int j=1;j<=n;j++)
			M.a[i][j]=a[i][j];
	Gauss();
	getpoly();
	F=f[1];
	invFn=poww(F[n],mod-2);
	Mod();
	Matrix now,ans;
	now.unit();
	for(int i=0,t=min(n,(int)G.size());i<t;i++)
	{
		ans=add(ans,mul(now,G[i]));
		now=mul(now,M);
	}
	for(int i=1;i<=n;i++)
	{
		for(int j=1;j<=n;j++)
			printf("%d ",ans.a[i][j]);
		puts("");
	}
	return 0;
}
/*
000000 2
1 2
3 4
*/