任意模数FFT：MTT

前言

注明：本文litble写给自己以后复习用，只是梳理，不做证明，学习MTT推荐下面这篇博客。

感谢这位dalao->here

MTT

将每个系数拆成$a+bM$的形式，将原来的一个多项式写成两个多项式。

则$(a_1+b_{1}M)(a_2+b_{2}M)=a_1{a}_2+a_1{b}_{2}M+b_1{a}_{2}M+b_1{b}_2{M}^2$

当$M=\sqrt{P}$时，$a_1{a}_2$，$a_1{b}_2$，$b_1{a}_2$，$b_1{b}_2$都是$P$级别的，就可以防止炸掉。

不过这样的话，要DFT$a_1$，$a_2$，$b_1$，$b_2$这四个多项式，要IDFT$a_1{a}_2$，$a_1{b}_2+b_1{a}_2$，$b_1{b}_2$，一共做了7次（虽然还是比任意模数NTT少了2次的说……），不够优秀。

于是就有大神提出了优秀的做法。

定义$P(x)=A(x)+B(x)i$，$Q(x)=A(x)-B(x)i$。

则会有$Q({\omega }^k)=conj(P({\omega }^{-k}))=conj(P({\omega }^{n-k}))$

（conj是取共轭复数的意思，$n$是FFT时我们把多项式的度数延伸到的那个2的次方级别的数）

于是我们只用DFT这个$P(x)$，就能得到$Q(x)$被DFT之后的结果。

DFT完成后，我们还可以还原出$A$和$B$的DFT：

$$A({\omega }^k)=\frac{P({\omega }^k)+Q({\omega }^k)}{2}$$

$$B({\omega }^k)=-i\frac{P({\omega }^k)-Q({\omega }^k)}{2}$$

这样我们就可以处理出$a=a_1{a}_2$，$b=a_1{b}_2$，$c=b_1{a}_2$，$d=b_1{b}_2$的说，注意这里的$a,b,c,d$还都是复数的说。

然后构造$M_1=a+bi$和$M_2=c+di$两个多项式，IDFT回去。则分别取这两个IDFT后的多项式的实部和虚部（第一个的实部为$a$，虚部为$b$以此类推），就是IDFT后的$a,b,c,d$，节约2的常数。

这样只要做四次DFT即可。

注意本方法对精度要求很高，记得开long double

代码

洛谷P4245

#include<bits/stdc++.h>
using namespace std;
#define RI register int
int read() {
	int q=0;char ch=' ';
	while(ch<'0'||ch>'9') ch=getchar();
	while(ch>='0'&&ch<='9') q=q*10+ch-'0',ch=getchar();
	return q;
}
typedef long double db;//对精度要求很高
typedef long long LL;
const int N=262150,M=32767;
const db pi=acos(-1);
int n,m,P,ka[N],kb[N],kc[N],rev[N];
struct com{db r,i;}a[N],b[N],Aa[N],Ab[N],Ba[N],Bb[N];
com operator + (com A,com B) {return (com){A.r+B.r,A.i+B.i};}
com operator - (com A,com B) {return (com){A.r-B.r,A.i-B.i};}
com operator * (com A,com B) {return (com){A.r*B.r-A.i*B.i,A.r*B.i+A.i*B.r};}
com conj(com A) {return (com){A.r,-A.i};}

void FFT(com *a,int n) {
	for(RI i=0;i<n;++i) if(rev[i]>i) swap(a[i],a[rev[i]]);
	for(RI i=1;i<n;i<<=1) {
		com wn=(com){cos(pi/i),sin(pi/i)};
		for(RI j=0;j<n;j+=(i<<1)) {
			com t1,t2,w=(com){1,0};
			for(RI k=0;k<i;++k,w=w*wn)
				t1=a[j+k],t2=w*a[j+i+k],a[j+k]=t1+t2,a[j+i+k]=t1-t2;
		}
	}
}
void mul(int *ka,int *kb,int *kc,int n) {
	for(RI i=0;i<n;++i) ka[i]%=P,kb[i]%=P;
	for(RI i=0;i<n;++i)
		a[i]=(com){ka[i]&M,ka[i]>>15},b[i]=(com){kb[i]&M,kb[i]>>15};//构造出两个要相乘的多项式的P
	FFT(a,n),FFT(b,n);
	for(RI i=0;i<n;++i) {
		int j=(n-i)&(n-1);//将P还原回A,B。别漏了&(n-1)这个取模操作
		com kAa=(a[i]+conj(a[j]))*(com){0.5,0};
		com kAb=(a[i]-conj(a[j]))*(com){0,-0.5};
		com kBa=(b[i]+conj(b[j]))*(com){0.5,0};
		com kBb=(b[i]-conj(b[j]))*(com){0,-0.5};
		Aa[j]=kAa*kBa,Ab[j]=kAa*kBb,Ba[j]=kAb*kBa,Bb[j]=kAb*kBb;
		//提前将多项式翻转完成，IDFT的时候就不用翻了
	}
	for(RI i=0;i<n;++i)//构造M
		a[i]=Aa[i]+Ab[i]*(com){0,1},b[i]=Ba[i]+Bb[i]*(com){0,1};
	FFT(a,n),FFT(b,n);
	for(RI i=0;i<n;++i) {
		int kAa=(LL)(a[i].r/n+0.5)%P;//直接提取实部和虚部
		int kAb=(LL)(a[i].i/n+0.5)%P;
		int kBa=(LL)(b[i].r/n+0.5)%P;
		int kBb=(LL)(b[i].i/n+0.5)%P;
		kc[i]=((LL)kAa+((LL)(kAb+kBa)<<15)+((LL)kBb<<30))%P;
	}
}
int main()
{
	n=read(),m=read(),P=read();
	for(RI i=0;i<=n;++i) ka[i]=read();
	for(RI i=0;i<=m;++i) kb[i]=read();
	int kn=1,len=0;
	while(kn<=n+m) kn<<=1,++len;
	for(RI i=0;i<kn;++i) rev[i]=(rev[i>>1]>>1)|((i&1)<<(len-1));
	mul(ka,kb,kc,kn);
	for(RI i=0;i<=n+m;++i) printf("%d ",(kc[i]+P)%P);
	return 0;
}