bzoj 4162 shlw loves matrix II

Posted on | Views:

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

$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)$ .

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
#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;
}