codeforces1054H. Epic Convolution 2dFFT 原根

codeforces1054H. Epic Convolution

题目传送门

分析

直接处理不好处理，考虑转化。
难搞的是$c^{i^2{j}^3}$
根据费马小定理，可以让$i^2{j}^3\mathrm{m}\mathrm{o}\mathrm{d}490018$
发现模数较小，一个基本的思路是，考虑求${\sum }_{i,j}{a}_i{b}_j[i^2{j}^3\equiv k\mathrm{m}\mathrm{o}\mathrm{d}490018]$。
乘法不好搞，考虑用原根把乘法转化成加法。
我们假设模数是质数$P$，原根是$w$，设$i=w^x,j=w^y,k=w^z$
那么问题转化成
${\sum }_{i,j}{a}_i{b}_j[2x+3y\equiv z\mathrm{m}\mathrm{o}\mathrm{d}P-1]$
只需要将所有的$a_i$映射到$a_{2x}$，所有的$b_j$映射到$b_{3y}$上就可以用$FFT$解决了。
不是质数直接搞肯定不成，于是考虑分解一下模数然后用$CRT$合并
$490018=2\cdot 491\cdot 499$
模数很给力，刚好仨质数。
将一个数$i$分解成$(i\mathrm{m}\mathrm{o}\mathrm{d}2,i\mathrm{m}\mathrm{o}\mathrm{d}491,i\mathrm{m}\mathrm{o}\mathrm{d}499)$
$\mathrm{m}\mathrm{o}\mathrm{d}2$这个玩意儿不用太在意，考虑分奇偶性。奇数的答案一定是奇数卷奇数得到的，偶数的答案用全部的答案减去奇数的答案就好了。
于是每个数$i$按奇数和偶数分类后转化成二元组$(i\mathrm{m}\mathrm{o}\mathrm{d}491,i\mathrm{m}\mathrm{o}\mathrm{d}499)$
然后映射成原根$(g_1^{i_1},g_2^{i_2})$
注意这个时候如果是$499,491$的倍数要单独考虑，直接抓出来暴力就可以了，复杂度过得去。
接下来我们要解决的问题变成了一个二维卷积。
也就是
$A_{a,b}={\sum }_{c,d,e,f}[c+e\equiv a\mathrm{m}\mathrm{o}\mathrm{d}490][d+f\equiv b\mathrm{m}\mathrm{o}\mathrm{d}499]a_{c,d}\cdot b_{e,f}$
看上去很麻烦。不过我们考虑优化一维。
把每一行看成一个多项式，那么相当于$A_a={\sum }_{c,e}[c+e\equiv a\mathrm{m}\mathrm{o}\mathrm{d}490]a_c\ast b_e$
为了方便我们用$\ast$表示卷积，$\cdot$表示乘法，外面的$\sum$是多项式加法。
考虑先对行$DFT$，那么有$[DFT(A_a){]}_k={\sum }_{c,e}[c+e\equiv a\mathrm{m}\mathrm{o}\mathrm{d}490][DFT(a_c){]}_k\cdot [DFT(b_e){]}_k$
这个时候问题转化成了每一列的卷积，只不过卷积的东西是对行多项式$DFT$完之后的点值。
同样地，我们可以再对列进行$FFT$，这样的话就可以转化成
$A_{a,b}=a_{a,b}\cdot b_{a,b}$
$IDFT$是完全一样的一个过程。
实际上上面描述的算法叫做$2d-FFT$，它的公式是：
$$DFT::f(x,y)=\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}F(u,v)e^{2\pi i(\frac{xu}{M}+\frac{yv}{N})}$$
$$IDFT::F(x,y)=\frac{1}{NM}\sum _{u=0}^{M-1}\sum _{v=0}^{N-1}f(u,v)e^{-2\pi i(\frac{xu}{M}+\frac{yv}{N})}$$
然后问题就解决啦！
复杂度$O(K)$

代码

代码略长。。。但是应该会比标解的400+行好一丢丢。

#include<bits/stdc++.h>
int ri() {
	char c = getchar(); int x = 0, f = 1; for(;c < '0' || c > '9'; c = getchar()) if(c == '-') f  = -1;
	for(;c >= '0' && c <= '9'; c = getchar()) x = (x << 1) + (x << 3) - '0' + c; return x * f;
}
const int L = 1024, N = 1e5 + 10, P = 490019, phi = P - 1;
const double pi = acos(-1.0);
struct cp {
    double r, i;
    cp(double _r = 0, double _i = 0) : r(_r), i(_i) {}
    cp operator + (const cp &a) {return cp(r + a.r, i + a.i);}
    cp operator - (const cp &a) {return cp(r - a.r, i - a.i);}
    cp operator * (const cp &a) {return cp(r * a.r - i * a.i, i *a.r + a.i * r);}
    cp operator / (const double &x) {return cp(r / x, i / x);}
    void operator *= (const cp &x) {*this = *this * x;}
    void operator /= (const double &x) {*this = *this / x;}
}A[2][L][L], B[2][L][L], w[L];
int R[L], pwc[phi], st[N], a[N], b[N], pw1[500], Lg1[500], pw2[500], Lg2[500], c, tp, n, m, Ans;
int Pow(int x, int k, int P) {
	int r = 1;
	for(;k; k >>= 1, x = 1LL * x * x % P)
		if(k & 1) r = 1LL * r * x % P;
	return r;
}
int Iv(int x, int P) {return Pow(x, P - 2, P);}
void Chi(int i, int &x, int &y, int &z)  {x = i & 1; y = i % 491; z = i % 499;}
int iChi(int a, int b, int c) {
	int x = 491LL * 499 * Iv(491 * 499, 2) % phi;
	int y = 2LL * 499 * Iv(2 * 499, 491) % phi;
	int z = 2LL * 491 * Iv(2 * 491, 499) % phi;
	return (1LL * a * x + 1LL * b * y + 1LL * c * z) % phi;
}
void Pre() {
    int x = 9;
    for(int i = 1;i < L; ++i)
        R[i] = R[i >> 1] >> 1 | (i & 1) << x;
    for(int i = 0;i < L; ++i)     
        w[i] = cp(cos(2 * pi * i / L), sin(2 * pi * i / L));
	pwc[0] = 1;
	for(int i = 1;i < phi; ++i)
		pwc[i] = 1LL * pwc[i - 1] * c % P;
	pw1[0] = 1;
	for(int i = 1;i < 490; ++i)
		pw1[i] = pw1[i - 1] * 7 % 491,
		Lg1[pw1[i]] = i;
	pw2[0] = 1;
	for(int i = 1;i < 498; ++i)
		pw2[i] = pw2[i - 1] * 7 % 499,
		Lg2[pw2[i]] = i;
}
void DFT(cp *F, int ri) {
    for(int i = 1;i < L; ++i)
        if(i < R[i])
            std::swap(F[i], F[R[i]]);
    for(int i = 1, d = L >> 1;i < L; i <<= 1, d >>= 1)
        for(int j = 0;j < L; j += i << 1) {
            cp *l = F + j, *r = F + j + i, *p = w, tp;
            for(int k = i; k--; ++l, ++r, p += d)
                tp = *r * *p, *r = *l - tp, *l = *l + tp;
        }
    if(!~ri) {
    	for(int i = 1;i < L >> 1; ++i)
    		std::swap(F[i], F[L - i]);
    	for(int i = 0;i < L; ++i)
    		F[i] /= L;
    }
}
void Tw_DFT(cp F[L][L], int ir) {
	for(int i = 0;i < L; ++i)
		DFT(F[i], ir);
	for(int i = 0;i < L; ++i) {
		static cp tp[L];
		for(int j = 0;j < L; ++j)
			tp[j] = F[j][i];
		DFT(tp, ir);
		for(int j = 0;j < L; ++j)	
			F[j][i] = tp[j];
	}
}
void Solve0() {
	tp = 0;
	for(int i = 0;i < std::max(n, m); i += 491)
		st[++tp] = i;
	for(int i = 0;i < std::max(n, m); i += 499)
		st[++tp] = i;
	std::sort(st + 1, st + tp + 1);
	tp = std::unique(st + 1, st + tp + 1) - st - 1;
	for(int i = 0;i < n; ++i) {
		int i2 = 1LL * i * i % phi; long long sum = 0;
		for(int j = 1;j <= tp && st[j] < m; ++j) {
			int j3 = 1LL * st[j] * st[j] * st[j] % phi;
			sum += 1LL * a[i] * b[st[j]] * pwc[1LL * i2 * j3 % phi] % P;
		}
		Ans = (Ans + sum) % P;
	}
	for(int j = 0;j < m; ++j) {
		int j3 = 1LL * j * j * j % phi; long long sum = 0;
		for(int i = 1;i <= tp && st[i] < n; ++i) {
			int i2 = 1LL * st[i] * st[i] % phi; 
			sum += 1LL * a[st[i]] * b[j] * pwc[1LL * i2 * j3 % phi] % P;
		}
		Ans = (Ans + sum) % P;
	}
	for(int i = 1;i <= tp && st[i] < n; ++i) {
		int i2 = 1LL * st[i] * st[i] % phi; long long sum = 0;
		for(int j = 1;j <= tp && st[j] < m; ++j) {
			int j3 = 1LL * st[j] * st[j] * st[j] % phi;
			sum += 1LL * a[st[i]] * b[st[j]] * pwc[1LL * i2 * j3 % phi] % P;
		}
		Ans = (Ans - sum) % P;
	}
}
int main() {
	n = ri(); m = ri(); c = ri();
	for(int i = 0;i < n; ++i)
		a[i] = ri();
	for(int i = 0;i < m; ++i)
		b[i] = ri();
	Pre(); 
	Solve0();
	for(int i = 0;i < n; ++i) {
		int x, y, z; Chi(i, x, y, z);
		if(!y || !z) continue;
		y = (Lg1[y] << 1) % 490;
		z = (Lg2[z] << 1) % 498;
		A[x][y][z].r += a[i];
		if(x)
			A[0][y][z].r += a[i];
	}
	for(int j = 0;j < m; ++j) {
		int x, y, z; Chi(j, x, y, z);
		if(!y || !z) continue;
		y = (Lg1[y] * 3) % 490;
		z = (Lg2[z] * 3) % 498;
		B[x][y][z].r += b[j];
		if(x)
			B[0][y][z].r += b[j];
	}
	Tw_DFT(A[0], 1); Tw_DFT(A[1], 1);
	Tw_DFT(B[0], 1); Tw_DFT(B[1], 1);
	for(int x = 0;x < 2; ++x)
		for(int y = 0;y < L; ++y)
			for(int z = 0;z < L; ++z)
				A[x][y][z] *= B[x][y][z];
	Tw_DFT(A[0], -1); Tw_DFT(A[1], -1);
	for(int x = 0;x < 2; ++x)
		for(int y = 0;y < L; ++y)
			for(int z = 0;z < L; ++z) {
				long long val = (long long)(A[x][y][z].r + 0.5);
				if(!x)
					val -= (long long)(A[1][y][z].r + 0.5);
				val %= P;
				Ans = (Ans + val * pwc[iChi(x, pw1[y % 490], pw2[z % 498])]) % P;
			}
	printf("%d\n", (Ans + P) % P);
	return 0;
}