0. 楔子
我们知道，\(\mathtt{FFT}\) 是用复数运算的，当我们需要取模且数据范围较大的时候就没有办法了，比如 这道题。
如果能减小数据范围，我们就可以先算再模，于是，拆系数 \(\mathtt{FFT}\) 就闪亮登场了！
1. 正文
首先计算一下本题的数据范围：\(10^9\times 10^9\times 10^5=10^{23}\)（有 \(a_{x+y}\leftarrow a_xa_y\)，而且这样的 \(x,y\) 的组合是 \(n\) 级别的）。
设 \(\text{base}\) 为一个 \(\sqrt p\) 级别的数，那么类似除法的定义，可以将 \(F(x)\) 分解成 \(\text{base}\cdot A(x)+B(x)\)，这样拆出来的函数的系数是 \(<\text{base}\) 的。
故
这样数据范围就是 \(10^{14}\) 左右，可以用 \(\text{long long}\) 存储。但是这样要做 \(7\) 次 \(\mathtt{FFT}\)（插值 \(4\) 次，转换 \(3\) 次）。为啥不能将三者的点值表达式直接相加呢？显然第一个与第二个点值表达式乘上 \(\rm base\) 就直接爆了，而点值表达式又不能取模，所以只有转回去，取模再相乘、相加。
我们知道，在进行朴素 \(\mathtt{FFT}\) 时虚部并未被使用，我们可以利用这个空间。
考虑计算 \(A,B,C,D\) 的点值表示。令（注意下文将 \(f(\omega_n^k)\) 简写成 \(f(k)\)）
我们发现 \(f(k),g(n-k)\) 是共轭的（\(n\) 是 \(n\) 补全的 \(2\) 的幂）：
需要注意的是，\(f(k),g(n-k)\) 共轭也就代表着 \(f(n)\) 与 \(g(0)\) 配对，所以 \(g(0)\) 需要特殊计算。
综上，计算 \(f,g,h\) 的点值表达式只用 \(2\) 次 \(\mathtt{FFT}\).
令 \(p=f*h,q=g*h\). 所以
将 \(p,q\) 转回系数多项式需要 \(2\) 次 \(\mathtt{FFT}\)，将 \(p,q\) 对应项相加即可解出 \(A*C,A*D\)，从而都解出来。
总共需要 \(4\) 次 \(\mathtt{FFT}\).
2. 代码实现
#include <cstdio>
#define print(x,y) write(x),putchar(y)
template <class T>
inline T read(const T sample) {
T x=0; char s; bool f=0;
while((s=getchar())>'9' || s<'0')
f |= (s=='-');
while(s>='0' && s<='9')
x = (x<<1)+(x<<3)+(s^48),
s = getchar();
return f?-x:x;
}
template <class T>
inline void write(T x) {
static int writ[50],tp=0;
if(x<0) putchar('-'),x=-x;
do writ[++tp] = x-x/10*10, x/=10; while(x);
while(tp) putchar(writ[tp--]^48);
}
#include <cmath>
#include <iostream>
using namespace std;
typedef long long ll;
const int maxn = 4e5+5;
const double PI = acos(-1.0);
const ll bas = (1<<15), bas2 = bas*bas;
struct cp {
double x,y;
cp() {}
cp(const double X,const double Y):x(X),y(Y) {}
cp operator + (const cp& t) { return cp(x+t.x,y+t.y); }
cp operator - (const cp& t) { return cp(x-t.x,y-t.y); }
cp operator * (const cp& t) {
return cp(x*t.x-y*t.y,x*t.y+y*t.x);
}
} f[maxn],g[maxn],h[maxn],w[2][maxn],tmp;
int n,m,mod,lim=1,bit;
int rev[maxn];
void preWork() {
while(lim<n+m-1) lim<<=1,++bit;
for(int i=0;i<lim;++i) {
rev[i] = (rev[i>>1]>>1)|((i&1)<<bit-1);
w[1][i] = cp(cos(PI*2*i/lim),sin(PI*2*i/lim));
w[0][i] = cp(w[1][i].x,-w[1][i].y);
}
}
void FFT(cp* f,bool opt=1) {
for(int i=0;i<lim;++i)
if(i<rev[i]) swap(f[i],f[rev[i]]);
for(int mid=1;mid<lim;mid<<=1)
for(int i=0;i<lim;i+=(mid<<1))
for(int j=0;j<mid;++j) {
tmp = w[opt][lim/mid/2*j]*f[i+j+mid];
f[i+j+mid] = f[i+j]-tmp;
f[i+j] = f[i+j]+tmp;
}
}
int main() {
int val;
n=read(9)+1,m=read(9)+1,mod=read(9);
preWork();
for(int i=0;i<n;++i)
val = read(9), f[i] = cp(val>>15,val&32767);
for(int i=0;i<m;++i)
val = read(9), h[i] = cp(val>>15,val&32767);
FFT(f),FFT(h);
g[0] = cp(f[0].x,-f[0].y);
for(int i=1;i<lim;++i)
g[i] = cp(f[lim-i].x,-f[lim-i].y);
for(int i=0;i<lim;++i) {
h[i].x /= lim, h[i].y /= lim; // 只用除一次~
f[i] = f[i]*h[i];
g[i] = g[i]*h[i];
}
FFT(f,0),FFT(g,0);
for(int i=0;i<n+m-1;++i) {
int AC = ((ll)((f[i].x+g[i].x)/2+0.5))%mod;
int AD = ((ll)((f[i].y+g[i].y)/2+0.5))%mod;
int BD = ((ll)round(g[i].x-AC)%mod+mod)%mod;
int BC = ((ll)round(f[i].y-AD)%mod+mod)%mod;
print((bas2*AC+bas*(AD+BC)+BD)%mod,' ');
}
return 0;
}
