lp4245 【模板】任意模数NTT

我们首先考虑一个很暴力的玩法,直接找一个很大很大很大的模数,然后用int128上一个朴素NTT,再将结果对题目给的模数取模,似乎就可以了。
然而仔细考虑一下发现这个做法实在太暴力了(才不是因为找不到合适的模数呢。)我们考虑另一种方案。
现在假设我们已经用好几种模数求出最终的各种值了,那么原值是多少呢?很显然可以用中国剩余定理来求出。
我们计算了一下,发现只要用三个质数,值域就足以用中国剩余定理求出最终值了。
然而我们知道,模运算没有交换律。故而我们不能用常规的中国剩余定理合并,应该需要使用一种类似于扩展中国剩余定理的递推合并方法来合并它们的解。
这就完成了任意模数的NTT。
另外要注意数组大小…

#include<iostream>
#include<cstdio>

inline void Swap(int &A,int &B){
	A^=B^=A^=B;
}

inline long long mlt(long long A,long long X,long long P){
	long long BS=A,RT=0;
	while(X){
		if(X&1){
			RT+=BS;
			RT%=P;
		}
		BS+=BS;
		BS%=P;
		X>>=1;
	}
	return RT;
}

inline long long pw(long long A,long long X,long long P){
	long long BS=A,RT=1;
	while(X){
		if(X&1){
			RT*=BS;
			RT%=P;
		}
		BS*=BS;
		BS%=P;
		X>>=1;
	}
	return RT;
}

const long long MOD[3]={469762049,998244353,1004535809};
const long long MM=468937312667959297;
const long long g0=3,gi[3]={156587350,332748118,334845270};

int L=1,inv[3],R[1<<21|1];
inline void prpr(int LEN){
	int B=0;
	while(L<=LEN){
		L<<=1;
		++B;
	}
	for(int i=0;i<3;++i){
		inv[i]=MOD[i]-(MOD[i]-1)/L;
	}
	for(int i=0;i<L;++i){
		R[i]=R[i>>1]>>1|(i&1)<<(B-1);
	}
}

inline void NTT(int *A,int typ,int P){
	for(int i=0;i<L;++i){
		if(R[i]<i){
			Swap(A[i],A[R[i]]);
		}
	}
	int bs,nw,X,Y,M;
	for(int i=2;i<=L;i<<=1){
		M=i>>1;
		bs=pw(~typ?g0:gi[P],(MOD[P]-1)/i,MOD[P]);
		for(int j=0;j<L;j+=i){
			nw=1;
			for(int k=0;k<M;++k,nw=1ll*nw*bs%MOD[P]){
				X=A[j+k],Y=1ll*nw*A[j+k+M]%MOD[P];
				A[j+k]=(X+Y)%MOD[P];
				A[j+k+M]=(X-Y+MOD[P])%MOD[P];
			}
		}
	}
}
int C[1<<21|1],D[1<<21|1],ans[3][1<<21|1];
inline void FNTT(int *A,int *B,int P){
	for(int i=0;i<L;++i){
		C[i]=A[i],D[i]=B[i];
	}
	NTT(C,1,P);
	NTT(D,1,P);
	for(int i=0;i<L;++i){
		C[i]=1ll*C[i]*D[i]%MOD[P];
	}
	NTT(C,-1,P);
	for(int i=0;i<L;++i){
		ans[P][i]=(int)(1ll*(C[i]+MOD[P])*inv[P]%MOD[P]);
	}
}
int a[1<<21|1],b[1<<21|1],n,m,p;
long long t[3];
void init(){
	scanf("%d%d%d",&n,&m,&p);
	for(int i=0;i<=n;++i){
		scanf("%d",a+i);
	}
	for(int i=0;i<=m;++i){
		scanf("%d",b+i);
	}
	prpr(n+m);
	FNTT(a,b,0);FNTT(a,b,1);FNTT(a,b,2);
	t[0]=pw(MOD[1]%MOD[0],MOD[0]-2,MOD[0]);t[1]=pw(MOD[0]%MOD[1],MOD[1]-2,MOD[1]);t[2]=pw(MM%MOD[2],MOD[2]-2,MOD[2]);//分别求出要用到的三个逆元。 
	long long T1,T2;
	for(int i=0;i<=n+m;++i){
		T1=(mlt(1ll*ans[0][i]*MOD[1]%MM,t[0],MM)+mlt(1ll*ans[1][i]*MOD[0]%MM,t[1],MM))%MM;//直接用CRT合并一二式。 
		T2=((ans[2][i]-T1)%MOD[2]+MOD[2])%MOD[2]*t[2]%MOD[2];//用EXCRT将第三式合并之。 
		printf("%d ",(int)(((T2%p)*(MM%p)%p+T1%p)%p));//求得最终解。 
	}
}

int main(){
	init();
	return 0;
}

发表评论

邮箱地址不会被公开。