[CF1054H] Epic Convolution——数论、卷积、任意模数NTT

偶耶XJX

已于 2022-08-08 20:11:40 修改

阅读量244

点赞数 1

分类专栏：信息竞赛解题文章标签：算法数论卷积快速数论变换NTT

于 2022-07-28 21:08:19 首次发布

本文链接：https://blog.csdn.net/weixin_43960287/article/details/126024623

版权

信息竞赛解题专栏收录该内容

195 篇文章 5 订阅

订阅专栏

CF1054H Epic Convolution

题解

这道题叫 Epic Convolution，翻译过来就是“史诗卷积”（也叫“屎屎卷积”）

由于模数很小而且是个质数，所以由欧拉定理，我们可以把问题转化为求
$f_x=\sum_{i=0}^{n-1}\sum_{j=0}^{m-1}A_iB_j[i^2j^3\bmod 490018=x]$ 然后为了更高效地求解，我们需要考虑把上面的乘法卷积转化为加法卷积。

可 $490018$ 显然没有原根啊？这要怎么转换？

这个时候把值域剖成乘法环一类的想法显然是不现实的（这就是为什么我是数学的屎），我们需要更加抽象的想法。考虑把 $490018$ 做一个唯一分解：
$490018=2\times 491 \times 499$ 而当我们得到了 $i^2j^3$ 取模 $2, 491, 499$ 的值之后，可以用 CRT 合并起来，所以不妨改一下式子，变为求
$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]$
剩下就简单了，我们只需要对每一维依次卷积，其中第一维可以暴力卷积，后面两维都可以用原根取对数转化成加法卷积。

注意到可能出现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;
}