「Luogu P5293」[HNOI2019]白兔之舞

发表于 |
字数统计: 2.3k | 阅读时长 ≈ 12

题面有点长，而且还不好解释，所以不放简要题面了。

链接

Luogu P5293

前言

2019 年 4 月 7 日，HNOI2019 Day 2 赛场上，菜鸡 Karry5307 打了这题的 50 分部分分做法，可是下午出成绩的时候发现本题 boom 0。

2019 年 7 月，菜鸡 Karry5307 看到了 memset0 的题解，但是由于太菜，并没有看懂。

2019 年 10 月，菜鸡 Karry5307 学习了单位根反演，可是并没有学出什么所以然。

2019 年 11 月，菜鸡 Karry5307 学习了 MTT，第一次失败了，但是最终还是成功的通过了 Luogu 的模板题。

2020 年 1 月 21 日，菜鸡 Karry5307 成功 AC CF1286F，多项式 16 题只剩下 CF 901E 了，但是 CF901E 由于需要用到 Bluestein’s Algorithm，所以他放弃了。

2020 年 2 月 27 日，菜鸡 Karry5307 看到神仙 x义x 讨论性能优化那题，发现 x义x 会 Bluestein’s Algorithm 而自己不会。

2020 年 3 月 2 日，菜鸡 Karry5307 最终还是学习了 Bluestein’s Algorithm，并且在 30min 切掉此题。

一年前我在 HNOI2019 的赛场上，折戟沉沙。一年后，我从倒下的地方爬起。

我成功了。我不再是以前的那个我了。

题解

首先来看一个很傻逼的 $\texttt{dp}$。

考虑 $f_{i,j}$ 表示走 $i$ 步到第 $j$ 行的任意合法位置的方案数，$g_{i,j}$ 表示走 $i$ 步到第 $j$ 行的方案数（注意这里并没有钦定在哪个位置上）。

转移比较显然，枚举一下上一次在哪一行，有

然后注意到行数比较小，所以可以使用矩阵优化。

一开始初始矩阵是 $G_0=\begin{pmatrix}1&0&0\end{pmatrix}$（相当于白兔在第一行），转移矩阵是

所以 $G_i=\begin{pmatrix}g_{i,1}&\cdots&g_{i,n}\end{pmatrix}=G_0S^i$。

那么 $f_{i,j}$ 等价于白兔在 $1$ 到 $L$ 中选择 $i$ 个经过点进行跳舞，所以有

所以我们要求的是

发现这个不太好搞定，所以移项。

然后就是套路单位根反演。

然后把 $k$ 提出来，单位根的括号拆掉，有

交换一下求和符号就可以把单位根提出来了，顺便把之前的 $\texttt{dp}$ 式子代进来所以

也就是说，我们要求

的第 $y$ 项，也即

然后发现是个二项式定理。更加形式化的来讲，是这样的

所以

后面这一坨可以矩阵快速幂。

设 $h_{j}=G_0(\omega_{k}^{j}S+I)^L$ 的第 $y$ 项，那么

如果熟悉任意长度循环卷积的套路的话可以考虑 $\texttt{Bluestein’s Algorithm}$。

但是如果把 $ij$ 拆成 $\frac{(i+j)^2}{2}-\frac{j^2}{2}-\frac{j^2}{2}$ 的话会出现一个问题，在模意义下没有二次剩余。

所以考虑把 $ij$ 拆成 $\binom{i+j}{2}-\binom{i}{2}-\binom{j}{2}$。

比较显然，拆开就可以了，所以有

把一些奇奇怪怪的东西提到外面去

然后设 $f_j=h_j\omega_{k}^{\binom{j}{2}},g_j=\omega_k^{-\binom{t+j}{2}}$，有

这个随便搞搞就可以了，可以倍长也可以先 reverse 一下卷完 reverse 回来。

注意到，在模意义下的单位根是原根，所以要先求出原根，最后 MTT 即可。

代码

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
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
#include<bits/stdc++.h>
using namespace std;
typedef long long int ll;
typedef long double db;
const ll MAXN=524291;
const db PI=acos(-1);
struct Complex{
	db re,im;
	Complex(db x=0,db y=0)
	{
		this->re=x,this->im=y;
	}
	inline Complex conj()
	{
		return Complex(re,-im);
	}
};
struct Matrix{
	ll num[5][5];
	Matrix()
	{
		memset(num,0,sizeof(num));
	}
	inline ll* operator [](const ll &x)
	{
		return num[x];
	}
	inline const ll* operator [](const ll &x)const
	{
		return num[x];
	}
};
Matrix p,m;
ll n,kk,l,x,y,md,g,cnt,limit,zkdxl;
ll w[5][5];
ll rev[MAXN],pw[MAXN],f[MAXN],gg[MAXN],h[MAXN];
inline ll read()
{
    register ll num=0,neg=1;
    register char ch=getchar();
    while(!isdigit(ch)&&ch!='-')
    {
        ch=getchar();
    }
    if(ch=='-')
    {
        neg=-1;
        ch=getchar();
    }
    while(isdigit(ch))
    {
        num=(num<<3)+(num<<1)+(ch-'0');
        ch=getchar();
    }
    return num*neg;
} 
inline ll qpow(ll base,ll exponent,ll md)
{
	ll res=1;
	while(exponent)
	{
		if(exponent&1)
		{
			res=res*base%md;	
		}	
		base=base*base%md,exponent>>=1;
	}	
	return res;
} 
inline Complex operator +(Complex x,Complex y)
{
	return Complex(x.re+y.re,x.im+y.im);
}
inline Complex operator -(Complex x,Complex y)
{
	return Complex(x.re-y.re,x.im-y.im);
}
inline Complex operator *(Complex x,Complex y)
{
	return Complex(x.re*y.re-x.im*y.im,x.re*y.im+x.im*y.re);
}
inline Complex operator /(Complex x,db y)
{
	return Complex(x.re/y,x.im/y);
}
inline Matrix id()
{
	Matrix res;
	for(register int i=1;i<=n;i++)
	{
		res[i][i]=1;
	}
	return res;
}
inline Matrix operator *(Matrix x,Matrix y)
{
	Matrix res;
	for(register int i=1;i<=n;i++)
	{
		for(register int j=1;j<=n;j++)
		{
			for(register int k=1;k<=n;k++)
			{
				res[i][j]=(res[i][j]+x[i][k]*y[k][j]%md)%md;
			}
		}
	}
	return res;
}
inline Matrix qpow(Matrix base,ll exponent)
{
	Matrix res=id();
	while(exponent)
	{
		if(exponent&1)
		{
			res=res*base;
		}
		base=base*base,exponent>>=1;
	}
	return res;
}
inline ll getG(ll x)
{
	static ll q[10051];
	ll top=0,flg=1;
	for(register int i=2;i<=sqrt(x-1);i++)
	{
		if((x-1)%i==0)
		{
			q[++top]=i;
			if(i*i!=x-1)
			{
				q[++top]=(x-1)/i;
			}
		}
	}	
	for(register int i=2;;i++)
	{
		flg=1;
		for(register int j=1;j<=top;j++)
		{
			if(qpow(i,q[j],x)==1)
			{
				flg=0;
				break;
			}
		}
		if(flg)
		{
			return i;
		}
	}
	return -1;
}
inline void calcP(ll x)
{
	for(register int i=1;i<=n;i++)
	{
		for(register int j=1;j<=n;j++)
		{
			p[i][j]=w[i][j]*x%md;
		}
		p[i][i]++;
	}
}
inline ll comb(ll x)
{
	return x<=1?0:(x-1)*x/2;
}
inline void FFT(Complex *cp,ll cnt,ll inv)
{
	ll cur=0;
	Complex res,omg;
	for(register int i=0;i<cnt;i++)
	{
		if(i<rev[i])
		{
			swap(cp[i],cp[rev[i]]);
		}
	}
	for(register int i=2;i<=cnt;i<<=1)
	{
		cur=i>>1,res=Complex(cos(2*PI/i),inv*sin(2*PI/i));
		for(register Complex *p=cp;p!=cp+cnt;p+=i)
		{
			omg=Complex(1,0);
			for(register int j=0;j<cur;j++)
			{
				Complex t=omg*p[j+cur],t2=p[j];
				p[j+cur]=t2-t,p[j]=t+t2;
				omg=omg*res;
			}
		}
	}
	if(inv==-1)
	{
		for(register int i=0;i<cnt;i++)
		{
			cp[i]=cp[i]/cnt;
		}
	}
}
inline void conv(ll cnt,ll *f,ll *g,ll *res)
{
	Complex p,q,r,s;
	ll t;
	ll t0,t1,t2,t3;
	static Complex tmpf[MAXN],tmpg[MAXN],p0[MAXN],p1[MAXN],p2[MAXN],p3[MAXN];
	for(register int i=0;i<cnt;i++)
	{
		tmpf[i]=Complex(f[i]&32767,f[i]>>15);
		tmpg[i]=Complex(g[i]&32767,g[i]>>15);
	}
	FFT(tmpf,cnt,1),FFT(tmpg,cnt,1);
	for(register int i=0;i<cnt;i++)
	{
		t=(cnt-i)&(cnt-1);
		p=(tmpf[i]+tmpf[t].conj())*Complex(0.5,0);
		q=(tmpf[i]-tmpf[t].conj())*Complex(0,-0.5);
		r=(tmpg[i]+tmpg[t].conj())*Complex(0.5,0);
		s=(tmpg[i]-tmpg[t].conj())*Complex(0,-0.5);
		p0[i]=p*r,p1[i]=p*s,p2[i]=q*r,p3[i]=q*s;
	}
	for(register int i=0;i<cnt;i++)
	{
		tmpf[i]=p0[i]+p1[i]*Complex(0,1);
		tmpg[i]=p2[i]+p3[i]*Complex(0,1);
	}
	FFT(tmpf,cnt,-1),FFT(tmpg,cnt,-1);
	for(register int i=0;i<cnt;i++)
	{
		t0=(ll)(tmpf[i].re+0.5)%md,t1=(ll)(tmpf[i].im+0.5)%md;
		t2=(ll)(tmpg[i].re+0.5)%md,t3=(ll)(tmpg[i].im+0.5)%md;
		res[i]=(t0+(t1+t2)*32768%md+t3*1073741824ll%md)%md;
	}
}
int main()
{
	n=read(),kk=read(),l=read(),x=read(),y=read(),md=read();
	for(register int i=1;i<=n;i++)
	{
		for(register int j=1;j<=n;j++)
		{
			w[i][j]=read();
		}
	}
	g=qpow(getG(md),(md-1)/kk,md),pw[0]=1;
	for(register int i=1;i<kk;i++)
	{
		pw[i]=pw[i-1]*g%md;
	}
	for(register int i=0,omg=1;i<kk;i++,omg=omg*g%md)
	{
		calcP(omg),p=qpow(p,l),gg[i]=p[x][y]*pw[comb(i)%kk]%md;
	}
	for(register int i=0;i<=(kk<<1);i++)
	{
		f[i]=pw[(kk-comb(i)%kk)%kk];
	}
	reverse(f,f+((kk<<1)|1)),cnt=1,limit=-1;
	while(cnt<=((kk<<1)+kk+1))
	{
		cnt<<=1,limit++;
	}
	for(register int i=0;i<cnt;i++)
	{
		rev[i]=(rev[i>>1]>>1)|((i&1)<<limit);
	}
	conv(cnt,f,gg,h),zkdxl=qpow(kk,md-2,md);
	for(register int i=0;i<kk;i++)
	{
		printf("%lld\n",h[(kk<<1)-i]*zkdxl%md*pw[comb(i)%kk]%md);
	}
}