洛谷P4245 【模板】任意模数NTT

传送门

这里用的拆系数(FFT)（咱实在是不会三模……）

就是说把(F_i)和(G_i)中的每一个数字都拆成(k imes M+b)的形式，这里(M)可以取(2^{15}=32768)

于是把(F)拆出来的两个数列和(G)拆出来的两个数列分别相乘，最终对应的系数要乘上(1,2^{15},2^{15},2^{30})

//minamoto
#include<bits/stdc++.h>
#define R register
#define ll long long
#define fp(i,a,b) for(R int i=a,I=b+1;i<I;++i)
#define fd(i,a,b) for(R int i=a,I=b-1;i>I;--i)
#define go(u) for(int i=head[u],v=e[i].v;i;i=e[i].nx,v=e[i].v)
using namespace std;
char buf[1<<21],*p1=buf,*p2=buf;
inline char getc(){return p1==p2&&(p2=(p1=buf)+fread(buf,1,1<<21,stdin),p1==p2)?EOF:*p1++;}
int read(){
    R int res,f=1;R char ch;
    while((ch=getc())>'9'||ch<'0')(ch=='-')&&(f=-1);
    for(res=ch-'0';(ch=getc())>='0'&&ch<='9';res=res*10+ch-'0');
    return res*f;
}
char sr[1<<21],z[20];int K=-1,Z=0;
inline void Ot(){fwrite(sr,1,K+1,stdout),K=-1;}
void print(R int x){
    if(K>1<<20)Ot();if(x<0)sr[++K]='-',x=-x;
    while(z[++Z]=x%10+48,x/=10);
    while(sr[++K]=z[Z],--Z);sr[++K]=' ';
}
const int N=5e5+5;const double Pi=acos(-1.0);
struct cp{
	double x,y;
	cp(double xx=0,double yy=0){x=xx,y=yy;}
    inline cp operator +(cp b)const{return cp(x+b.x,y+b.y);}
    inline cp operator -(cp b)const{return cp(x-b.x,y-b.y);}
    inline cp operator *(cp b)const{return cp(x*b.x-y*b.y,x*b.y+y*b.x);}
	inline cp operator *(const double &b)const{return cp(x*b,y*b);}
}A[N],B[N],C[N],D[N],H[N],F[N],G[N],w[N];
int r[N],lim=1,l,n,m,P,x;
void FFT(cp *A,int ty){
	fp(i,0,lim-1)if(i<r[i])swap(A[i],A[r[i]]);
	for(R int mid=1;mid<lim;mid<<=1)
		for(R int j=0;j<lim;j+=(mid<<1))
			for(R int k=0;k<mid;++k){
				cp x=A[j+k],y=w[mid+k]*A[j+k+mid];
				A[j+k]=x+y,A[j+k+mid]=x-y;
			}
	if(ty==-1){
		reverse(A+1,A+lim);
		double k=1.0/lim;fp(i,0,lim-1)A[i]=A[i]*k;
	}
}
int main(){
//	freopen("testdata.in","r",stdin);
	n=read(),m=read(),P=read();
	while(lim<=n+m)lim<<=1,++l;
	fp(i,0,lim-1)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	for(R int i=1;i<lim;i<<=1)fp(k,0,i-1)w[i+k]=cp(cos(Pi*k/i),sin(Pi*k/i));
	fp(i,0,n)x=read(),A[i].x=x>>15,B[i].x=x&32767;
	fp(i,0,m)x=read(),C[i].x=x>>15,D[i].x=x&32767;
	FFT(A,1),FFT(B,1),FFT(C,1),FFT(D,1);
	fp(i,0,lim-1)
		F[i]=A[i]*C[i],G[i]=A[i]*D[i]+B[i]*C[i],H[i]=B[i]*D[i];
	FFT(F,-1),FFT(G,-1),FFT(H,-1);
	fp(i,0,n+m)print((((ll)(F[i].x+0.5)%P<<30)+((ll)(G[i].x+0.5)<<15)+((ll)(H[i].x+0.5)))%P);
	return Ot(),0;
}