CF1054H Epic Convolution
题解
这道题叫 Epic Convolution,翻译过来就是“史诗卷积”(也叫“屎屎卷积”)
由于模数很小而且是个质数,所以由欧拉定理,我们可以把问题转化为求
f x = ∑ i = 0 n − 1 ∑ j = 0 m − 1 A i B j [ i 2 j 3 m o d 490018 = x ] f_x=\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}A_iB_j[i^2j^3\bmod 490018=x] fx=i=0∑n−1j=0∑m−1AiBj[i2j3mod490018=x]然后为了更高效地求解,我们需要考虑把上面的乘法卷积转化为加法卷积。
可 490018 490018 490018 显然没有原根啊?这要怎么转换?
这个时候把值域剖成乘法环一类的想法显然是不现实的(这就是为什么我是数学的屎),我们需要更加抽象的想法。考虑把 490018 490018 490018 做一个唯一分解:
490018 = 2 × 491 × 499 490018=2\times 491 \times 499 490018=2×491×499而当我们得到了 i 2 j 3 i^2j^3 i2j3 取模 2 , 491 , 499 2,491,499 2,491,499 的值之后,可以用 CRT 合并起来,所以不妨改一下式子,变为求
f x , y , z = ∑ i = 0 n − 1 ∑ j = 0 m − 1 A i B j [ i 2 j 3 m o d 2 = x ] [ i 2 j 3 m o d 491 = y ] [ i 2 j 3 m o d 499 = z ] f_{x,y,z}=\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}A_iB_j[i^2j^3\bmod 2=x][i^2j^3\bmod 491=y][i^2j^3\bmod 499=z] fx,y,z=i=0∑n−1j=0∑m−1AiBj[i2j3mod2=x][i2j3mod491=y][i2j3mod499=z]
剩下就简单了,我们只需要对每一维依次卷积,其中第一维可以暴力卷积,后面两维都可以用原根取对数转化成加法卷积。
注意到可能出现0无法取对数的情况,在0处暴力卷积即可。
代码
我写的是暴力分类讨论的卷积,而实际上用一点小容斥可以显著减少运算次数。
#include<bits/stdc++.h>//JZM yyds!!
#define ll long long
#define lll __int128
#define uns unsigned
#define fi first
#define se second
#define IF (it->fi)
#define IS (it->se)
#define lowbit(x) ((x)&-(x))
#define END putchar('\n')
#define inline jzmyyds
using namespace std;
const int MAXN=5e5+5;
const ll INF=1e17;
ll read(){
ll x=0;bool f=1;char s=getchar();
while((s<'0'||s>'9')&&s>0){
if(s=='-')f^=1;s=getchar();}
while(s>='0'&&s<='9')x=(x<<1)+(x<<3)+(s^48),s=getchar();
return f?x:-x;
}
int ptf[50],lpt;
void print(ll x,char c='\n'){
if(x<0)putchar('-'),x=-x;
ptf[lpt=1]=x%10;
while(x>9)x/=10,ptf[++lpt]=x%10;
while(lpt)putchar(ptf[lpt--]^48);
if(c>0)putchar(c);
}
//490018=2*491*499,g_491=2,g_499=7
const ll MOD=490019,NM=1000000000245104641ll;//一遍大模数NTT
const lll E=1;
ll ksml(ll a,ll b,ll mo=NM){
ll res=1;
for(;b;b>>=1,a=E*a*a%mo)if(b&1)res=E*res*a%mo;
return res;
}
ll ksm(ll a,ll b,ll mo=MOD){
ll res=1;
for(;b;b>>=1,a=a*a%mo)if(b&1)res=res*a%mo;
return res;
}
#define g 7ll
int rev[1025];
int initrev(){
for(int i=0;i<(1<<10);i++)rev[i]=(rev[i>>1]>>1)|((i&1)<<9);
return 114514;
}
ll omg[1025],cbddl=initrev(),IVN=ksml(1024,NM-2);
void NTT(ll*a,int inv,int n=1024){
for(int i=0;i<n;i++)if(i<rev[i])swap(a[i],a[rev[i]]);
ll x,y,tmp=233,mi=(NM-1)>>1;omg[0]=1;
for(int m=1;m<n;m<<=1,mi>>=1){
if(m>1)tmp=ksml(g,inv>0?mi:NM-1-mi);
for(int i=1;i<m;i++)omg[i]=E*omg[i-1]*tmp%NM;
for(int i=0,om=0;i<n;i+=(m<<1),om=0)
for(int j=i;j<i+m;j++,om++){
x=a[j],y=E*a[j+m]*omg[om]%NM,a[j]=x+y,a[j+m]=x-y;
if(a[j]>=NM)a[j]-=NM;
if(a[j+m]<0)a[j+m]+=NM;
}
}if(inv<0)for(int i=0;i<n;i++)a[i]=E*a[i]*IVN%NM;
}
#undef g
ll f[1025][1025],g[1025][1025];
void polymul(){
for(int i=0;i<500;i++)NTT(f[i],1),NTT(g[i],1);
for(int i=1;i<1024;i++)
for(int j=0;j<i;j++)swap(f[i][j],f[j][i]),swap(g[i][j],g[j][i]);
for(int i=0;i<1024;i++){
NTT(f[i],1),NTT(g[i],1);
for(int j=0;j<1024;j++)f[i][j]=E*f[i][j]*g[i][j]%NM;
NTT(f[i],-1);
}
for(int i=1;i<1024;i++)for(int j=0;j<i;j++)swap(f[i][j],f[j][i]);
for(int i=0;i<1000;i++)NTT(f[i],-1);
}
int ex1[810],lg1[810],ex2[810],lg2[810];
int initexp(){
ex1[0]=ex2[0]=1,lg1[1]=lg2[1]=0;
for(int i=1;i<490;i++)ex1[i]=(ex1[i-1]<<1)%491,lg1[ex1[i]]=i;
for(int i=1;i<498;i++)ex2[i]=ex2[i-1]*7%499,lg2[ex2[i]]=i;
return 1919810;
}
int sh_t=initexp(),n,m;
ll a1[499],b1[499],ans,mi[MAXN],A[MAXN];
void solve(ll(*a)[499],ll(*b)[499],ll(*c)[499],bool CL){
if(CL){
memset(f,0,sizeof(f));
memset(g,0,sizeof(g));
memset(a1,0,sizeof(a1));
memset(b1,0,sizeof(b1));
}
for(int i=1;i<491;i++)
for(int j=1;j<499;j++)
(f[lg1[i]][lg2[j]]+=a[i][j])%=MOD,(g[lg1[i]][lg2[j]]+=b[i][j])%=MOD;
for(int i=1;i<491;i++)
for(int j=1;j<499;j++)(a1[i]+=a[i][j])%=MOD;
for(int i=1;i<491;i++)
for(int j=0;j<499;j++)(b1[i]+=b[i][j])%=MOD;
polymul();
for(int i=0,mi=0;i<1000;i++,mi=i%490)
for(int j=0;j<1000;j++)
(c[ex1[mi]][ex2[j%498]]+=f[i][j])%=MOD;
for(int i=1;i<491;i++)
for(int j=1,to=i;j<491;j++,to=(to+i)%491)
(c[to][0]+=a[i][0]*b1[j]+a1[i]*b[j][0])%=MOD;
memset(a1,0,sizeof(a1));
memset(b1,0,sizeof(b1));
for(int i=1;i<491;i++)
for(int j=0;j<499;j++)(a1[j]+=a[i][j])%=MOD;
for(int i=0;i<491;i++)
for(int j=0;j<499;j++)(b1[j]+=b[i][j])%=MOD;
for(int i=0;i<499;i++)
for(int j=0,to=0;j<499;j++,to=(to+i)%499)
(c[0][to]+=a[0][i]*b1[j]+a1[i]*b[0][j])%=MOD;
}
int crt(int a,int b,int c){
static const ll iv1=492>>1,iv2=ksm(491<<1,499-2,499);
int k1=(b-a+491)*iv1%491,k2=(c-a-k1-k1+4990)*iv2%499;
return (a+(k1<<1)+(k2*491<<1))%(MOD-1);
}
ll a[2][491][499],b[2][491][499],b_[491][499],as[2][491][499],c;
const bool KFC=1;
int main()
{
n=read(),m=read(),c=read(),mi[0]=1;
for(int i=1;i<MOD;i++)mi[i]=mi[i-1]*c%MOD;
if(1ll*n*m<=1000000ll&&KFC){
for(int i=0;i<n;i++)A[i]=read();
for(int j=0;j<m;j++){
ll B=read();
for(int i=0;i<n;i++)
(ans+=A[i]*B*mi[1ll*i*i*j%(MOD-1)*j*j%(MOD-1)])%=MOD;
}return print(ans),0;
}
for(int i=0,k=0;i<n;i++,k=1ll*i*i%(MOD-1))
(a[k&1][k%491][k%499]+=read())%=MOD;
for(int i=0,k=0,d;i<m;i++,k=1ll*i*i*i%(MOD-1))
d=read(),(b[k&1][k%491][k%499]+=d)%=MOD,(b_[k%491][k%499]+=d)%=MOD;
solve(a[1],b[1],as[1],0);
solve(a[0],b_,as[0],1);
solve(a[1],b[0],as[0],1);
for(int i=0;i<2;i++)for(int j=0;j<491;j++)
for(int k=0;k<499;k++)(ans+=as[i][j][k]*mi[crt(i,j,k)])%=MOD;
print(ans);
return 0;
}
文章评论