CF1054H Epic Convolution

Link
首先我们应该确定的一件事是不需要MTT，模数很小所以精度问题并不大。
因为模数很小，而且根据Euler定理(c^{i^2j^3}equiv c^{i^2j^3mod490018}pmod{490019})，那么我们可以考虑求出：
(s_k=sumlimits_{i=0}^{n-1}A_isumlimits_{j=0}^{m-1}B_i[i^2j^3equiv kpmod{490018}])
最后答案就是(sumlimits s_kc^k)。
分解质因数发现(490018=2*491*499)，因此我们可以CRT合并。
也就是说现在我们要求的是：
(S_{x,y,z}=sumlimits_{i=0}^{n-1}A_isumlimits_{j=0}^{m-1}B_i[i^2j^3equiv xpmod2][i^2j^3equiv ypmod{491}][i^2j^3equiv zpmod{499}])
(x)这一维暴力就行了，(y,z)这两维可以先把(491,499)的倍数拿出来单独搞搞，剩下的用指标化成加。
打表可以发现(2,7)分别是(491,499)的最小原根。
设：
(P_{x,y,z}=sumlimits_{i=0}^{n-1}[i^2equiv xpmod 2][gamma_{491,2}(i^2)=y][gamma_{499,7}(i^2)=z]A_i,Q_{x,y,z}=sumlimits_{i=0}^{m-1}[i^3equiv xpmod 2][gamma_{491,2}(i^3)=y][gamma_{499,7}(i^3)=z]B_i)。
那么(S)就可以用上面两个式子的卷积表示：
(S_{x,y,z}=sumlimits_{a*b=x}sumlimits_{i+j=y}sumlimits_{p+q=z}P_{a,i,p}Q_{b,j,q})
第一维暴力容斥一下，后两维直接2维FFT就行了。

#include<cmath>
#include<cstdio>
#include<cctype>
#include<vector>
#include<cstring>
#include<algorithm>
namespace IO
{
    char ibuf[(1<<21)+1],*iS,*iT;
    char Get(){return (iS==iT? (iT=(iS=ibuf)+fread(ibuf,1,(1<<21)+1,stdin),(iS==iT? EOF:*iS++)):*iS++);}
    int read(){int x=0,c=Get();while(!isdigit(c))c=Get();while(isdigit(c))x=x*10+c-48,c=Get();return x;}
}
using IO::read;
using ld=double;
const ld pi=2*acos(-1);
const int P=490019,pa=491,pb=499,ga=2,gb=7,N=100007,M=1025;
void inc(int&a,int b){a+=b-P,a+=a>>31&P;}
void dec(int&a,int b){a-=b,a+=a>>31&P;}
int mul(int a,int b){return 1ll*a*b%P;}
int pow(int a,int k){int r=1;for(;k;k>>=1,a=mul(a,a))if(k&1)r=mul(a,r);return r;}
struct complex{ld x,y;complex(ld a=0,ld b=0){x=a,y=b;}}w[N],A[2][M][M],B[2][M][M],a[M],b[M];
complex operator+(complex a,complex b){return {a.x+b.x,a.y+b.y};}
complex operator-(complex a,complex b){return {a.x-b.x,a.y-b.y};}
complex operator/(complex a,ld x){return {a.x/x,a.y/x};}
complex operator*(complex a,complex b){return {a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x};}
void operator-=(complex&a,complex b){a=a-b;}
void operator+=(complex&a,complex b){a=a+b;}
int lim=1024,rev[M],sa[P],sb[P],ia[pa],ib[pb],s[P],p[2][M][M];
std::vector<int>vec;
void init()
{
    static ld x;
    for(int i=1;i<lim;++i) rev[i]=(rev[i>>1]>>1)|(i&1? 512:0);
    for(int i=0;i<512;++i) x=pi*i/lim,w[512+i]={cos(x),sin(x)};
    for(int i=511;i;--i) w[i]=w[i<<1];
}
void FFT(complex*a,int f)
{
    static complex x;
    if(!~f) std::reverse(a+1,a+lim);
    for(int i=0;i<lim;++i) if(i<rev[i]) std::swap(a[i],a[rev[i]]);
    for(int i=1;i<lim;i<<=1) for(int j=0,d=i<<1;j<lim;j+=d) for(int k=0;k<i;++k) x=a[i+j+k]*w[i+k],a[i+j+k]=a[j+k]-x,a[j+k]+=x;
    if(!~f) for(int i=0;i<lim;++i) a[i]=a[i]/lim;
}
int main()
{
    int n=read(),m=read(),c=read(),ans=0;init();
    for(int i=0;i<n;++i) inc(sa[1ll*i*i%(P-1)],read());
    for(int i=0;i<m;++i) inc(sb[1ll*i*i*i%(P-1)],read());
    for(int i=0;i<P-1;++i) if(!(i%pa)||!(i%pb)) vec.push_back(i);
    for(int i=0;i<P-1;++i) if(sb[i]) for(int j:vec) if(sa[j]) inc(s[1ll*i*j%(P-1)],mul(sb[i],sa[j]));
    for(int i=0;i<P-1;++i) if(sa[i]&&i%pa&&i%pb) for(int j:vec) if(sb[j]) inc(s[1ll*i*j%(P-1)],mul(sa[i],sb[j]));
    for(int i=0,x=1;i<pa-1;++i) ia[x]=i,x=x*ga%pa;
    for(int i=0,x=1;i<pb-1;++i) ib[x]=i,x=x*gb%pb;
    for(int i=0;i<P;++i)
	if(i%pa&&i%pb)
	{
	    int x=i%2,y=ia[i%pa],z=ib[i%pb];
	    A[x][y][z]=sa[i],B[x][y][z]=sb[i],p[x][y][z]=i;
	}
    for(int i=0;i<pa;++i) for(int j=0;j<pb;++j) A[0][i][j]+=A[1][i][j],B[0][i][j]+=B[1][i][j];
    for(int i=0;i<2;++i) for(int j=0;j<pa;++j) FFT(A[i][j],1),FFT(B[i][j],1);
    for(int i=0;i<2;++i)
	for(int k=0;k<lim;++k)
	{
	    for(int j=0;j<lim;++j) a[j]=A[i][j][k],b[j]=B[i][j][k];
	    FFT(a,1),FFT(b,1);
	    for(int j=0;j<lim;++j) a[j]=a[j]*b[j];
	    FFT(a,-1);
	    for(int j=0;j<lim;++j) A[i][j][k]=a[j];
	}
    for(int i=0;i<2;++i) for(int j=0;j<lim;++j) FFT(A[i][j],-1);
    for(int i=0;i<lim;++i) for(int j=0;j<lim;++j) A[0][i][j]-=A[1][i][j];
    for(int i=0;i<2;++i) for(int j=0;j<lim;++j) for(int k=0;k<lim;++k) inc(s[p[i][j%(pa-1)][k%(pb-1)]],(long long)(A[i][j][k].x+0.5)%P);
    for(int i=0,x=1;i<P;++i) inc(ans,mul(x,s[i])),x=mul(x,c);
    printf("%d",ans);
}