三次变两次 FFT
我们发现：
[(F(x)+iG(x))^2=F(x)^2-G(x)^2+2iF(x)G(x)
]
也就是说，我们把 (F(x)) 作为实部，(G(x)) 作为虚部，那么它的平方的虚部的 (1/2) 就是 (F(x)G(x))
可惜精度比较低。
四次 FFT 求任意模数多项式乘法
假设我们要求 (M(x) imes N(x)pmod{p})，因为如果我们直接 FFT 就会爆 double ，所以我们可以把 (M(x)) 拆成 (kA(x)+B(x))，(N(x)) 拆成 (kC(x)+D(x))，其中 (kapprox sqrt{p})，那么，你的值域就大约变为 (p imes n) 的了，但是你就需要 (7) 次 FFT 了。
我们假设有 (Q(x)=A(x)+iB(x),E(x)=C(x)+iD(x))，设 (q(x_1)) 表示 (Q(x)) 在 (x_1) 的取值，(A_j) 表示 (A(x)) 第 (j) 项系数，那么我们就有：
[q(omega_n^x)=sum_{j=0} A_jomega_n^{xj}+iB_jomega_n^{xj}
]
[=sum_{j=0} A_j(cos{frac{2pi xj}{n}}+isin{frac{2pi xj}{m}})+B_j(icos{frac{2pi xj}{n}}-sin{frac{2pi xj}{n}})
]
[q(omega_n^{-x})=sum_{j=0} A_j(cos{frac{2pi xj}{n}}-isin{frac{2pi xj}{n}})+B_j(icos{frac{2pi xj}{n}}+sin{frac{2pi xj}{m}})
]
[a(omega_n^{x})=sum_{j=0} A_j(cos{frac{2pi xj}{n}}+isin{frac{2pi xj}{n}})
]
[b(omega_n^{x})=sum_{j=0} B_j(cos{frac{2pi xj}{n}}+isin{frac{2pi xj}{n}})
]
我们假设 (q(x).r) 表示它的实部，(q(x).f) 表示它的虚部，那么我们就可以得到：
[a(omega_n^{x}).r=frac{q(omega_n^x).r+q(omega_n^{-x}).r}{2},a(omega_n^{x}).f=frac{q(omega_n^x).f-q(omega_n^{-x}).f}{2}
]
[b(omega_n^{x}).r=frac{q(omega_n^{-x}).f+q(omega_n^x).f}{2},a(omega_n^{x}).f=frac{q(omega_n^{-x}).r-q(omega_n^x).r}{2}
]
然后，我们求出 (A(x)E(x)) 和 (B(x)E(x))，我们就可以得到 (A(x)C(x),A(x)D(x),B(x)C(x),B(x)D(x))，然后算就好了。
暂时没有懂 (3.5) 次 FFT 的做法，所以不写了。
Code
#include <bits/stdc++.h>
using namespace std;
#define double long double
#define Int register int
#define MAXN 270005
template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c == '-') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}
struct Complex{
double x,y;
Complex(){}
Complex (double _x,double _y){x = _x,y = _y;}
Complex operator / (const double &p)const{return Complex{x / p,y / p};}
Complex operator + (const Complex &p)const{return Complex{x + p.x,y + p.y};}
Complex operator - (const Complex &p)const{return Complex{x - p.x,y - p.y};}
Complex operator * (const Complex &p)const{return Complex{x * p.x - y * p.y,x * p.y + p.x * y};}
};
#define pi (double)acos(-1)
int l,lim,rev[MAXN];
void fft (Complex *a,int type){
for (Int i = 0;i < lim;++ i) if (i < rev[i]) swap (a[i],a[rev[i]]);
for (Int i = 1;i < lim;i <<= 1){
Complex Wn(cos(pi / i),type * sin(pi / i));
for (Int j = 0,r = i << 1;j < lim;j += r){
Complex w(1,0);
for (Int k = 0;k < i;++ k,w = w * Wn){
Complex x = a[j + k],y = w * a[i + j + k];
a[j + k] = x + y,a[i + j + k] = x - y;
}
}
}
if (type == -1) for (Int i = 0;i < lim;++ i) a[i] = a[i] / lim;
}
int n,m,mod,a[MAXN],b[MAXN],ans[MAXN];
Complex Q[MAXN],E[MAXN],C[MAXN],D[MAXN];
#define ll long long
signed main(){
read (n,m,mod),lim = 1;
while (lim < n + m) lim <<= 1,++ l;
for (Int i = 0;i < lim;++ i) rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << l - 1);int up = (1 << 15) - 1;
for (Int i = 0;i <= n;++ i) read (a[i]),Q[i] = Complex (a[i] >> 15,a[i] & up);
for (Int i = 0;i <= m;++ i) read (b[i]),E[i] = Complex (b[i] >> 15,b[i] & up);
fft (Q,1),fft (E,1);
for (Int i = 0;i < lim;++ i){
int re = (lim - 1) & (lim - i);
C[i] = Complex((Q[i].x + Q[re].x) / 2,(Q[i].y - Q[re].y) / 2) * E[i];
D[i] = Complex((Q[re].y + Q[i].y) / 2,(Q[re].x - Q[i].x) / 2) * E[i];
}
fft (C,-1),fft (D,-1);
for (Int i = 0;i < lim;++ i){
ll v1 = (ll)(C[i].x + 0.5) % mod,v2 = (ll)(C[i].y + D[i].x + 0.5) % mod,v3 = (ll)(D[i].y + 0.5) % mod;
ans[i] = ((v1 << 30) + (v2 << 15) + v3) % mod;
}
for (Int i = 0;i <= n + m;++ i) write ((ans[i] % mod + mod) % mod),putchar (' ');
putchar ('
');
return 0;
}