「CodeForces 1054H」Epic Convolution

给两个下标从$0$开始且长度分别为$n,m$的序列$A,B$和一个常数$C$,求:

答案对$490019$取模。

$\texttt{Data Range:}n,m\leq 10^5,A_i,B_i\leq 10^3,c\leq 490019$

链接

Luogu

vjudge

题解

从昨天晚上写到今天晚上,现在终于 AC 了,写篇题解纪念一下。

MNz29I.md.png

不得不说这是我做过的最有意思并且最难的多项式题,实在配得上出题人给的 Epic 的称号和 CF 给的 $3400$ 的评分。

考虑到$\texttt{\color{black}y\color{red}wy}$的题解写的有一点点乱,我来写一篇比较好懂的(当然还是要感谢他的啦qwq)

这个题目与古代猪文那个题有点像,因为$490019$是质数,由欧拉定理,我们只需要求

发现这个指数只有$490018$中可能,于是我们可以有一个想法

看这个乘法挺讨厌的,显然可以取个离散对数化乘为加,然后发现不晓得以什么为底,考虑换个思路。

诶,我们发现这个$490018=2\times 491\times 499$是个$\texttt{square-free number}$,好啊。

很明显我们可以$\texttt{CRT}$合并,因为模数互质,我们就可以不用上拓中了。

于是我们可以自然地设一个这样的东西

然后我们发现一个有意思的事情:$491$有$g=2$,$499$有$g=7$,所以我们可以以原根为底啊,这就非常好了。

但是这个$2$没有原根啊,这怎么办?

考虑到这一维取值只有$2$种,可以暴力讨论,剩下两个东西可以离散对数搞搞。

于是我们可以设两个东西:(以下用$\log_{x,y}p$表示$p$在模$x$意义下以$y$为底的离散对数)

然后我们就会发现我们的要求的东西就可以写成这个东西:

二元多项式卷积了解一下

我们定义两个矩阵$A,B$的卷积$C$为$C_{i,j}=\sum\limits_{i=0}^{n}\sum\limits_{j=0}^{m}A_{i,j}B_{n-i,m-j}$。

考虑到第一维模数为$2$的特殊性质,我们可以大力讨论简单容斥一下,有

接着就是二维$\texttt{FFT}$的表演啦qwq。

你可以去网上百度这个东西然后发现所谓二维$\texttt{FFT}$其实就是把每一行做一遍$\texttt{FFT}$再每一列也做一遍$\texttt{FFT}$,也可以反过来。

证明由于篇幅原因所限其实是我不会略去。

然后发现某一些数可能是$491$或者$499$的倍数导致没有离散对数,这些数可以暴力。

于是我们就做完了qwq。

然后还有一个附加环节,就是卡常。

能少模的就少模,开个$\texttt{long long}$在一次模掉!

一定记得吸氧!不要用万能头!

然后你就能从 TLE 13 到 AC 了qwq

代码

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
#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<vector>
#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize("Ofast")
using namespace std;
typedef int ll;
typedef long long int li;
typedef long double db;
const ll MAXN=1e5+51,G499=7,G491=2,MOD=490019;
const db PI=acos(-1.0);
struct Complex{
db re,im;
Complex(db x=0,db y=0)
{
this->re=x,this->im=y;
}
};
Complex f0,f1,g0,g1;
Complex omgs[1025],invo[1025],tmp[1024];
Complex f[2][1024][1024],g[2][1024][1024];
vector<ll>xx,yy;
ll n,m,c,kk=250,kk2=399,cur,p,q,res,h,xp;
li kd,qq;
ll rev[1024],x[MAXN],y[MAXN];
ll lg499[500],lg491[500],pw499[500],pw491[500];
li fx[2][1024][1024],fy[2][1024][1024];
li pw[490021];
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 CRT(ll x,ll y,ll z)
{
ll r=kk*(y-x);
r=((r%499)+499)%499;
ll p=2*r+x,r2=kk2*(z-p);
r2=((r2%491)+491)%491;
return r2*998+p;
}
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 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=inv==1?omgs[i]:invo[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 FFT2D(ll n,ll m,Complex cp[][1024],ll inv)
{
for(register int i=0;i<m;i++)
{
for(register int j=0;j<n;j++)
{
tmp[j]=cp[j][i];
}
FFT(tmp,n,inv);
for(register int j=0;j<n;j++)
{
cp[j][i]=tmp[j];
}
}
for(register int i=0;i<n;i++)
{
FFT(cp[i],m,inv);
}
}
int main()
{
n=read(),m=read(),c=read(),cur=1;
for(register int i=0;i<490;i++)
{
lg491[cur]=i,pw491[i]=cur,cur=(cur*G491)%491;
}
cur=1;
for(register int i=0;i<498;i++)
{
lg499[cur]=i,pw499[i]=cur,cur=(cur*G499)%499;
}
for(register int i=0;i<1024;i++)
{
rev[i]=(rev[i>>1]>>1)|((i&1)<<9);
}
for(register int i=0;i<n;i++)
{
x[i]=read();
if(i%499&&i%491)
{
fx[i&1][lg499[(li)i*i%499]][lg491[(li)i*i%491]]+=x[i];
continue;
}
xx.push_back(i);
}
for(register int i=0;i<m;i++)
{
y[i]=read();
if(i%499&&i%491)
{
fy[i&1][lg499[(li)i*i*i%499]][lg491[(li)i*i*i%491]]+=y[i];
continue;
}
yy.push_back(i);
}
for(register int i=0;i<499;i++)
{
for(register int j=0;j<491;j++)
{
f[0][i][j].re=fx[0][i][j]%MOD,f[1][i][j].re=fx[1][i][j]%MOD;
g[0][i][j].re=fy[0][i][j]%MOD,g[1][i][j].re=fy[1][i][j]%MOD;
}
}
for(register int i=1;i<=1024;i<<=1)
{
omgs[i]=Complex(cos(2*PI/i),sin(2*PI/i));
invo[i]=Complex(cos(2*PI/i),-sin(2*PI/i));
}
FFT2D(1024,1024,f[0],1),FFT2D(1024,1024,f[1],1);
FFT2D(1024,1024,g[0],1),FFT2D(1024,1024,g[1],1);
for(register int i=0;i<1024;i++)
{
for(register int j=0;j<1024;j++)
{
f0=f[0][i][j],f1=f[1][i][j],g0=g[0][i][j],g1=g[1][i][j];
f[0][i][j]=(f0+f1)*(g0+g1)-f1*g1,f[1][i][j]=f1*g1;
}
}
FFT2D(1024,1024,f[0],-1),FFT2D(1024,1024,f[1],-1),pw[0]=1;
for(register int i=1;i<MOD;i++)
{
pw[i]=(li)pw[i-1]*c%MOD;
}
for(register int i=0;i<1024;i++)
{
for(register int j=0;j<1024;j++)
{
for(register int k=0;k<2;k++)
{
h=((li)(f[k][i][j].re+0.5))%MOD;
xp=CRT(k,pw499[i%498],pw491[j%490]);
res=(res+h*pw[xp])%MOD;
}
}
}
for(register int i=0;i<xx.size();i++)
{
for(register int j=0;j<m;j++)
{
kd=(li)xx[i]*xx[i]*j%(MOD-1)*j*j%(MOD-1);
qq=x[xx[i]]*y[j]*pw[kd]%MOD;
res=(li)(res+qq)%MOD;
}
}
for(register int i=0;i<yy.size();i++)
{
for(register int j=0;j<n;j++)
{
if(!(j%499&&j%491))
{
continue;
}
kd=(li)yy[i]*yy[i]*yy[i]%(MOD-1)*j*j%(MOD-1);
res=(res+(li)x[j]*y[yy[i]]*pw[kd])%MOD;
}
}
printf("%d\n",res);
}