给定两个多项式(f(x),g(x))，求(h(x)=f(x)*g(x))
对(p)取模，(p)不保证可以分解成(a*2^k+1)
三模NTT
由于模数不满足有原根，我们可以找几个有原根的模数，求出结果再(crt)合并一下
考虑卷积后数字最大能达到(p*p*len)，一般是(10^9*10^9*10^5=10^{23})，所以我们选择模数之积应该大于(10^{23})
一般考虑(998244353,1004535809,469762049)，因为原根都是(3)，在(int)范围内且乘积较大
然后跑三个(ntt)（天体常数）得到
[egin{cases}retequiv a_1 (mod p_1)\retequiv a_2 (mod p_2)\retequiv a_3 (mod p_3)\retequiv x (mod p)end{cases}
]
我们现在就是求(x)
由于前三个式子直接合并爆(long long)，我们用一些技巧：
先合并前两个，定义(inv(x,p))为(x)在模(p)意义下的逆元
[retequiv a_1*p_2*inv(p_2,p_1)+a_2*p_1*inv(p_1,p_2) (mod p_1*p_2)
]
记作
[retequiv d (mod m)
]
则
[ret=x*m+d=y*p_3+a_3
]
[xequiv (a_3-d)*M^{-1} (mod p_3)
]
设(q=(a_3-d)*M^{-1})，那么
[x=k*p_3+q
]
代入(ret)得到
[ret=k*p_1*p_2*p_3+q*M+d
]
然而(ansin [0,p_1*p_2*p_3))，所以(k=0)，(ret)就得出了
#include<bits/stdc++.h>
using namespace std;
namespace red{
#define int long long
#define eps (1e-8)
inline int read()
{
int x=0;char ch,f=1;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
const int N=4e5+10;
int mod[3]={469762049,998244353,1004535809};
int n,m,p;
int f[N],g[N],pos[N];
int b[N],ret[N];
int limit,len;
inline int fast(int x,int k,int p)
{
int ret=1;
while(k)
{
if(k&1) ret=ret*x%p;
x=x*x%p;
k>>=1;
}
return ret;
}
inline int slow(int x,int k,int p)
{
int ret=0;
while(k)
{
if(k&1) ret=(ret+x)%p;
x=(x+x)%p;
k>>=1;
}
return ret;
}
struct poly
{
int g=3,p,a[N];
inline void ntt(int limit,int *a,int inv)
{
for(int i=0;i<limit;++i)
if(i<pos[i]) swap(a[i],a[pos[i]]);
for(int mid=1;mid<limit;mid<<=1)
{
int Wn=fast(inv?g:(p+1)/g,(p-1)/(mid<<1),p);
for(int r=mid<<1,j=0;j<limit;j+=r)
{
int w=1;
for(int k=0;k<mid;++k,w=w*Wn%p)
{
int x=a[j+k],y=w*a[j+k+mid]%p;
a[j+k]=x+y;
if(a[j+k]>=p) a[j+k]-=p;
a[j+k+mid]=x-y;
if(a[j+k+mid]<0) a[j+k+mid]+=p;
}
}
}
if(inv) return;
inv=fast(limit,p-2,p);
for(int i=0;i<limit;++i) a[i]=a[i]*inv%p;
}
}fft[3];
inline int inv(int x,int p)
{
return fast(x%p,p-2,p);
}
inline void crt()
{
int len=n+m;
int M=mod[0]*mod[1];
int inv1=inv(mod[1],mod[0]),inv0=inv(mod[0],mod[1]),inv3=inv(M%mod[2],mod[2]);
int a,b,c,t,k;
for(int i=0;i<=len;++i)
{
a=fft[0].a[i],b=fft[1].a[i],c=fft[2].a[i];
t=(slow(a*mod[1]%M,inv1,M)+slow(b*mod[0]%M,inv0,M))%M;
k=((c-t%mod[2])%mod[2]+mod[2])%mod[2]*inv3%mod[2];
ret[i]=((k%p)*(M%p)%p+t%p)%p;
}
}
inline void main()
{
n=read(),m=read(),p=read();
for(int i=0;i<=n;++i) f[i]=read();
for(int i=0;i<=m;++i) g[i]=read();
for(limit=1;limit<=n+m+2;limit<<=1) ++len;
for(int i=0;i<limit;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
for(int k=0;k<=2;++k)
{
fft[k].p=mod[k];
for(int i=0;i<=n;++i) fft[k].a[i]=f[i];
for(int i=0;i<=m;++i) b[i]=g[i];
for(int i=m+1;i<limit;++i) b[i]=0;
fft[k].ntt(limit,fft[k].a,1);
fft[k].ntt(limit,b,1);
for(int i=0;i<limit;++i) fft[k].a[i]=fft[k].a[i]*b[i]%mod[k];
fft[k].ntt(limit,fft[k].a,0);
}
crt();
for(int i=0;i<=n+m;++i) printf("%lld ",ret[i]);
}
}
signed main()
{
red::main();
return 0;
}
MTT
毛爷爷数论变换
对于多项式(P=sumlimits_{i=0}^{n}P_ix^i,Q=sumlimits_{i=0}^{n}Q_ix^i)，求(P*Q)
考虑直接大力(FFT)，发现(10^{23})精度炸飞
考虑降低精度，我们设：
[A=sumlimits_{i=0}^{n}(P_i>>15)x^i,B=sumlimits_{i=0}^{n}(P_i&32767)x^i
]
[C=sumlimits_{i=0}^{n}(Q_i>>15)x^i,D=sumlimits_{i=0}^{n}(Q_i&32767)x^i
]
那么
[P*Q=AC*2^{30}+(AD+BC)*2^{15}+BD
]
需要(8)次(FFT)
我们设(F=A+iC,G=B+iD)
[T1=F*G=AB-CD+i(AD+BC)
]
再设(F^{'}=A-iC)
[T2=F^{'}*G=AB+CD+i(AD-BC)
]
其中
[T1+T2=2(AB+iAD)
]
[T2-T1=2(CD-iBC)
]
我们需要(3)次(DFT),(2)次(IDFT)，一共(5)次(FFT)
其实还能优化！
令
[P(x)=A(x)+iB(x)
]
[Q(x)=A(x)-iB(x)
]
设(P_t)是(P)的(DFT)，(conj(x))表示(x)的共轭复数，(A_i)为(A(x))的(i)次项系数
[P_t(k)=A(omega _n^k)+i B(omega _n^k)
]
[=sumlimits_{j=0}^{n-1}A_jomega _n^{jk}+iB_jomega _n^{jk}
]
[=sumlimits_{j=0}^{n-1}(A_J+iB_j)omega _n^{jk}
]
[Q_t(x)=A(omega _n^k)-i B(omega _n^k)
]
[=sumlimits_{j=0}^{n-1}A_jomega _n^{jk}-iB_jomega _n^{jk}
]
[=sumlimits_{j=0}^{n-1}(A_J-iB_j)omega _n^{jk}
]
[=sumlimits_{j=0}^{n-1}(A_j-i*B_j)(cos(frac{2pi jk}{n})+i*sin(frac{2pi jk}{n}))
]
[=sumlimits_{j=0}^{n-1}(A_jcos(frac{2pi jk}{n})+B_jsin(frac{2pi jk}{n}))+i(A_jsin(frac{2pi jk}{n})-B_jcos(frac{2pi jk}{n}))
]
[=conj(sumlimits_{j=0}^{n-1}(A_jcos(frac{2pi jk}{n})+B_jsin(frac{2pi jk}{n}))-i(A_jsin(frac{2pi jk}{n})-B_jcos(frac{2pi jk}{n})))
]
[=conj(sumlimits_{j=0}^{n-1}(A_jcos(frac{-2pi jk}{n})-B_jsin(frac{-2pi jk}{n}))-i(A_jsin(frac{-2pi jk}{n})+B_jcos(frac{-2pi jk}{n})))
]
[=conj(sumlimits_{j=0}^{n-1}(A_j+iB_j)(cos(frac{-2pi jk}{n})+i*sin(frac{-2pi jk}{n})))
]
[=conj(sumlimits_{j=0}^{n-1}(A_j+iB_j)omega _n^{-jk})
]
[=conj(sumlimits_{j=0}^{n-1}(A_j+iB_j)omega _n^{(n-j)k})
]
[=conj(P_t(n-k))
]
所以我们可以用(P)的点值得到(Q)，也就是说上面的(F^{'})不需要进行(FFT)
注意(n-k)当(k=0)时需要特殊处理，因为是循环卷积，需要移到第(0)项
得到这个结论我们只要(2)次(DFT)，(2)次(IDFT)，一共(4)次(FFT)
碾爆垃圾(9)次(NTT)
注意单位根会乘爆精度，考虑预处理，或者开(long double)
#include<bits/stdc++.h>
using namespace std;
namespace red{
#define int long long
#define eps (1e-8)
inline int read()
{
int x=0;char ch,f=1;
for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
if(ch=='-') f=0,ch=getchar();
while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
return f?x:-x;
}
const int N=1<<18;
const double pi=acos(-1.0);
int n,m,p;
int pos[N];
int limit,len;
struct complex
{
double x,y;
complex(double tx=0,double ty=0){x=tx,y=ty;}
inline complex operator + (const complex &t) const
{
return complex(x+t.x,y+t.y);
}
inline complex operator - (const complex &t) const
{
return complex(x-t.x,y-t.y);
}
inline complex operator * (const complex &t) const
{
return complex(x*t.x-y*t.y,x*t.y+y*t.x);
}
inline complex operator * (const double &t) const
{
return complex(x*t,y*t);
}
inline complex operator += (const complex &t) const
{
return complex(x+t.x,y+t.y);
}
inline complex operator ~ () const
{
return complex(x,-y);
}
}f[N],g[N],a[N],b[N],w[N];
inline void fft(int limit,complex *a,int inv)
{
for(int i=0;i<limit;++i)
if(i<pos[i]) swap(a[i],a[pos[i]]);
for(int mid=1;mid<limit;mid<<=1)
{
for(int r=mid<<1,j=0;j<limit;j+=r)
{
for(int k=0;k<mid;++k)
{
complex x=a[j+k],y=w[mid+k]*a[j+k+mid];
a[j+k]=x+y;
a[j+k+mid]=x-y;
}
}
}
}
inline void main()
{
n=read(),m=read(),p=read();
for(int x,i=0;i<=n;++i) x=read() , f[i].x=x>>15 , f[i].y=x&32767;
for(int x,i=0;i<=m;++i) x=read() , g[i].x=x>>15 , g[i].y=x&32767;
for(limit=1;limit<=n+m;limit<<=1) ++len;
for(int i=0;i<limit;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
for(int i=1;i<limit;i<<=1)
{
w[i]=(complex){1, 0};
for(int j=1;j<i;++j)
w[i+j]=((j&31)==1?(complex){cos(pi*j/i), sin(pi*j/i)}:w[i+j-1]*w[i+1]);
}
fft(limit,f,1);fft(limit,g,1);
for(int i=0;i<limit;++i)
{
static complex q,f0,f1,g0,g1;
q=~f[i?limit-i:0],f0=(f[i]-q)*(complex){0,-0.5},f1=(f[i]+q)*0.5;
q=~g[i?limit-i:0],g0=(g[i]-q)*(complex){0,-0.5},g1=(g[i]+q)*0.5;
a[i]=f1*g1,b[i]=f1*g0+f0*g1+f0*g0*(complex){0,1};
}
reverse(a+1,a+limit);reverse(b+1,b+limit);
fft(limit,a,-1),fft(limit,b,-1);
double k=1.0/limit;
for(int i=0;i<=n+m;++i)
{
printf("%lld ",(((int)(a[i].x*k+.5)%p<<30)+((int)(b[i].x*k+.5)<<15)+(int)(b[i].y*k+.5))%p);
}
}
}
signed main()
{
red::main();
return 0;
}