【学习笔记】MTT

用于任意模数多项式乘法。

题目链接

我们所得到的数最大是\(mod^2*Len,\)大概是\(10^{23}\)次方级别。这个时候朴素的\(FFT\)虽然支持取模但是精度会爆炸。

考虑将原来的多项式\(F(x)\)分解成\(A(x)M+B(x).\)这样可以使得这两个多项式的系数不至于过大。

\(A[i]=F[i]/M,B[i]=F[i]\bmod M.\)

这样，\(F(x)*G(x)=(A_0(x)M+B_0(x))(A_1(x)M+B_1(x))=A_0A_1(x)M^2+A_0B_1(x)+A_1B_0(x)+B_1B_0(x)\)

最后将四个多项式加起来即可。一共需要八次\(FFT.\)

#include<bits/stdc++.h>
using namespace std;
typedef long double D;
typedef long long ll;
inline int MAX(int x,int y){return x>y?x:y;}
const D Pi=acos(-1.0);
inline int read() {
	int f=1, r=0; char c=getchar();
	while(!isdigit(c)) { if(c=='-')f=-1; c=getchar(); }
	while(isdigit(c)) { r=r*10+c-'0'; c=getchar(); }
	return f*r;
}
const int MAXN=6e5+10;
struct cp{
	D x,y;
	cp(D xx=0,D yy=0){x=xx,y=yy;}
};
int l,r[MAXN],lt=1,n,m,mod;
inline int add(int x,int y){return (x+y)%mod;}
inline int mul(int x,int y){return 1ll*x*y%mod;}
int F[MAXN],G[MAXN],H[MAXN];
cp operator+(cp a,cp b){return cp(a.x+b.x,a.y+b.y);}
cp operator-(cp a,cp b){return cp(a.x-b.x,a.y-b.y);}
cp operator*(cp a,cp b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
void FFT(cp *P,int tp){
	for(int i=0;i<lt;++i)if(i<r[i])swap(P[i],P[r[i]]);
	for(int i=1;i<lt;i<<=1){
		cp Wn(cos(Pi/i),tp*sin(Pi/i));
		for(int j=0;j<lt;j+=(i<<1)){
			cp w(1,0);
			for(int k=0;k<i;++k,w=w*Wn){
				cp x=P[j+k],y=w*P[i+j+k];
				P[j+k]=x+y;
				P[i+j+k]=x-y;
			}
		}
	}
	if(tp==-1){
		for(int i=0;i<lt;++i)P[i].x/=lt;
	}
}
void MTT(int dg,int *A,int *B,int *C){
	cp a[MAXN],b[MAXN],c[MAXN],d[MAXN];
	cp e[MAXN],f[MAXN],g[MAXN],h[MAXN];
	while(lt<=(dg<<1))lt<<=1,l++;
	for(int i=0;i<lt;++i)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(int i=0;i<lt;++i){
		A[i]%=mod;B[i]%=mod;
		a[i]=A[i]>>15;b[i]=A[i]&0x7fff;
		c[i]=B[i]>>15;d[i]=B[i]&0x7fff;
	}
	FFT(a,1);FFT(b,1);FFT(c,1);FFT(d,1);
	for(int i=0;i<lt;++i){
		e[i]=a[i]*c[i];f[i]=b[i]*c[i];
		g[i]=a[i]*d[i];h[i]=b[i]*d[i];
	}
	FFT(e,-1);FFT(f,-1);FFT(g,-1);FFT(h,-1);
	for(int i=0;i<lt;++i){
		int A1=(((ll(round(e[i].x))%mod)<<30)%mod);
		int A2=((ll(round(f[i].x))%mod)<<15)%mod;
		int A3=((ll(round(g[i].x))%mod)<<15)%mod;
		int A4=ll(round(h[i].x))%mod;
		A1=add(A1,A2);A1=add(A1,A3);A1=add(A1,A4);
		C[i]=A1;
	}
}
int main(){
	n=read();m=read();mod=read();
	for(int i=0;i<=n;++i)F[i]=read();
	for(int i=0;i<=m;++i)G[i]=read();
	MTT(MAX(n,m),F,G,H);
	for(int i=0;i<=n+m;++i)printf("%d ",H[i]);
	return 0;
}