题目
题目链接：https://www.luogu.com.cn/problem/P4245
给定 (2) 个多项式 (F(x), G(x)) ，请求出 (F(x) * G(x))。
系数对 (p) 取模，且不保证 (p) 可以分解成 (p = a cdot 2^k + 1) 之形式。
(nleq 10^5;a,bleq 10^9;pleq 10^9+9)。
思路
不保证 (p) 是 NTT 模数，所以不能直接用 NTT 做。
一般有两种处理方法，一是用三模数 NTT，最后 CRT 合并，精度可达 (10^{26})。另一种是拆位 FFT。也称 MTT。
拆位 FFT 思路很简单，将多项式每一位系数拆成两部分 (a_1,a_2)，其中 (a_1) 是 (a) 在二进制下前十五位，(a_2) 是 (a) 在二进制下后十五位。
那么
[(F*G)[i]=(2^{15}a_1[i]+a_2[i]) imes (2^{15}b_1[i]+b_2[i])
]
[=2^{30}a_1[i]a_2[i]+2^{15}(a_1[i]b_2[i]+a_2[i]b_1[i])+a_2[i]b_2[i]
]
这样一顿乱搞之后，卷积后每一项系数不超过 (n imes (p^{frac{1}{2}})^2=npleq 10^{14})。可以用 long double 强行存下。
这样的话只需要 (4) 次 DFT 和 (3) 次 IDFT 即可。
但是依然有更优的办法。我们设复多项式 (F'[i]=a_1[i]+a_2[i]i)，(G'[i]=a_1[i]-a_2[i]i)，(H'[i]=b_1[i]+b_2[i]i)，那么
[(F'*H')[i]=a_1[i]b_1[i]-a_2[i]b_2[i]+i(a_1[i]b_2[i]+a_2[i]b_1[i])
]
[(G'*H')[i]=a_1[i]b_1[i]+a_2[i]b_2[i]+i(a_1[i]b_2[i]-a_2[i]b_1[i])
]
那么
[(F'*H')[i]+(G'*H')[i]=2a_1[i]b_1[i]+2a_1[i]b_2[i]i
]
[(F'*H')[i]-(G'*H')[i]=-2a_2[i]b_2[i]+2a_2[i]b_1[i]i
]
我们只需要 (3) 次 DFT 和 (2) 次 IDFT 即可。
注意为了保证精度，最好预处理单位根。
时间复杂度 (O(nlog n))。
事实上我写的 5 次 FFT 还没有 QuantAsk 7 次的做法快 /kk。
代码
#include <bits/stdc++.h>
#define cp complex<long double>
using namespace std;
typedef long long ll;
typedef long double ld;
const int N=300010;
const ld pi=acos(-1);
int n,m,lim,MOD,rev[N];
cp f[N],g[N],h[N],w[N];
void FFT(cp *f,int tag)
{
for (int i=0;i<lim;i++)
if (rev[i]<i) swap(f[rev[i]],f[i]);
w[0]=cp(1,0);
for (int k=1;k<lim;k<<=1)
{
cp tmp(cos(pi/k),tag*sin(pi/k));
for (int i=k-2;i>=0;i-=2)
w[i]=w[i>>1],w[i+1]=w[i]*tmp;
for (int i=0;i<lim;i+=(k<<1))
{
for (int j=0;j<k;j++)
{
cp x=f[i+j],y=w[j]*f[i+j+k];
f[i+j]=x+y; f[i+j+k]=x-y;
}
}
}
}
int main()
{
scanf("%d%d%d",&n,&m,&MOD);
for (int i=0,x;i<=n;i++)
{
scanf("%d",&x);
f[i]=cp(x>>15,x&32767);
g[i]=cp(x>>15,-(x&32767));
}
for (int i=0,x;i<=m;i++)
{
scanf("%d",&x);
h[i]=cp(x>>15,x&32767);
}
lim=1;
while (lim<=n+m) lim<<=1;
for (int i=0;i<lim;i++)
rev[i]=(rev[i>>1]>>1)|((i&1)?(lim>>1):0);
FFT(f,1); FFT(g,1); FFT(h,1);
for (int i=0;i<lim;i++)
f[i]*=h[i],g[i]*=h[i];
FFT(f,-1); FFT(g,-1);
for (int i=0;i<=n+m;i++)
{
ll a=(ll)((f[i]+g[i]).real()/2.0/lim+0.4999)%MOD;
ll b=(ll)((f[i]+g[i]).imag()/2.0/lim+0.4999)%MOD;
ll c=(ll)((f[i]-g[i]).imag()/2.0/lim+0.4999)%MOD;
ll d=(ll)((g[i]-f[i]).real()/2.0/lim+0.4999)%MOD;
printf("%lld ",(((1LL<<30)*a+(1LL<<15)*(b+c)+d)%MOD+MOD)%MOD);
}
return 0;
}