题解 P4245 【【模板】任意模数NTT】
command_block
·
·
题解
一种奇怪的拆系数FFT?需要5次FFT但是推导很简单,便于记忆,常数也不大。
卷积后，能产生的最大的数也就是mod*mod*len,实战应用中一般是10^9*10^9*10^5=10^{23}。
我们只要弄一个精度够高的多项式乘法,就能做到不丢精度。
一种想法是:跑九次NTT(三次乘法),把答案用中国剩余定理合并,精度可达10^{26}
这种做法常数实在过大,而且10^{26}使用long long存不下,还需要一些黑科技辅助,所以不推荐。
另一种想法是:拆系数FFT。
把一个数拆成a*2^{15}+b的形式,那么a,b<=2^{15}
把a和b弄成两个多项式,这样相乘的值域是n*mod≈10^{14}
那么c_1*c_2=(a_1*2^{15}+b_1)(a_2*2^{15}+b_2)
=a_1a_2*2^{30}+(a_1b_2+a_2b_1)2^{15}+b_1b_2
那么就是要做4次多项式乘法?12次FFT?常数爆炸!
冷静分析:每个多项式只用插值一次,共4次。
点乘复杂度忽略,最后得到的四个多项式各插值一次,共4次。
这样就是8次FFT,常数还是爆炸……
我们考虑推推式子来优化:
我们有四个多项式A1,A2,B1,B2,求这些多项式的两两乘积。
考虑到(a+bi)*(c+di)=a(c+di)*bi(c+di)=ac-bd+(ad+bc)i
那么我们设复多项式P=A1+iA2;\ Q=B1+iB2;
那么T1=P*Q=A1B1-A2B2+(A1B2+A2B1)i
我们又设P'=A1-iA2
那么T2=P'*Q=A1B1+A2B2+(A1B2-A2B1)i
总的FFT次数是:3次DFT+2次IDFT=5次.
不过,值域并不是美好的$10^{14}$,而是$10^{19}$,因为在IDFT出来之前还得除以$n$……
`long double`信仰跑吧,代码并不难写。
```cpp
#include<algorithm>
#include<cstdio>
#include<cmath>
#define Maxn 100500
#define lim 32000
const long double Pi=acos(-1);
using namespace std;
inline int read()
{
register int X=0;
register char ch=0;
while(ch<48||ch>57)ch=getchar();
while(ch>=48&&ch<=57)X=X*10+(ch^48),ch=getchar();
return X;
}
int n,m,p,tr[Maxn<<2];
struct CP
{
long double x,y;
CP operator + (const CP& B) const
{return (CP){x+B.x,y+B.y};}
CP operator - (const CP& B) const
{return (CP){x-B.x,y-B.y};}
CP operator * (const CP& B) const
{return (CP){x*B.x-y*B.y,x*B.y+y*B.x};}
}P1[Maxn<<2],P2[Maxn<<2],Q[Maxn<<2];
void FFT(CP *f,int op,int n)
{
for (int i=0;i<n;i++)
if (i<tr[i])swap(f[i],f[tr[i]]);
for(int p=2;p<=n;p<<=1){
int len=p>>1;
CP tmp=(CP){std::cos(Pi/len),op*std::sin(Pi/len)};
for(int k=0;k<n;k+=p){
CP buf=(CP){1,0};
for(int l=k;l<k+len;l++){
CP tt=buf*f[len+l];
f[len+l]=f[l]-tt;
f[l]=f[l]+tt;
buf=buf*tmp;
}//F(x)=FL(x^2)+x*FR(x^2)
//F(W^k)=FL(w^k)+W^k*FR(w^k)
//F(W^{k+n/2})=FL(w^k)-W^k*FR(w^k)
}
}
}
int ans[Maxn<<2];
int main()
{
scanf("%d%d%d",&n,&m,&p);n++;m++;
for (int i=0,sav;i<n;i++){
sav=read();
P1[i]=(CP){sav/lim,sav%lim};
P2[i]=(CP){sav/lim,-sav%lim};
}
for (int i=0,sav;i<m;i++){
sav=read();
Q[i]=(CP){sav/lim,sav%lim};
}
int len=1;
for (m=m+n-1;len<m;len<<=1);
for (int i=1;i<len;i++)
tr[i]=tr[i>>1]>>1|((i&1)?len>>1:0);
FFT(P1,1,len);FFT(P2,1,len);FFT(Q,1,len);
for (int i=0;i<len;i++){Q[i].x/=len;Q[i].y/=len;}
for (int i=0;i<len;i++)P1[i]=P1[i]*Q[i];
for (int i=0;i<len;i++)P2[i]=P2[i]*Q[i];
FFT(P1,-1,len);FFT(P2,-1,len);
for (int i=0;i<m;i++){
long long a1b1=0,a1b2=0,a2b1=0,a2b2;
a1b1=(long long)floor((P1[i].x+P2[i].x)/2+0.49)%p;
a1b2=(long long)floor((P1[i].y+P2[i].y)/2+0.49)%p;
a2b1=((long long)floor(P1[i].y+0.49)-a1b2)%p;
a2b2=((long long)floor(P2[i].x+0.49)-a1b1)%p;
ans[i]=((a1b1*lim+(a1b2+a2b1))*lim+a2b2)%p;
ans[i]=(ans[i]+p)%p;
printf("%d ",ans[i]);
}return 0;
}
```
如果被卡,考虑预处理单位根。
目测常数是$NTT$的4倍……