三模数NTT样例不过求助
查看原帖
三模数NTT样例不过求助
914961
liaoz123楼主2023/6/3 08:44

检查了一下NTT部分没有写错,CRT部分感觉也没错啊,求助!!!

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=4e5+5;
const int mod1=998244353,mod2=469762049,mod3=1004535809,gg=3;
int limit=1,r[N],l,f[N][4],g[N][4];
int qpow(int a,int b,int mod){
	int ss=1;
	while(b){
		if(b&1)ss=1ll*ss*a%mod;
		a=1ll*a*a%mod;
		b>>=1;
	}
	return ss;
}
void NTT(int *a,int inv,int mod){
	for(int i=0;i<limit;i++)
	if(i<r[i])swap(a[i],a[r[i]]);
	int ginv=qpow(gg,mod-2,mod);
	for(int mid=1;mid<limit;mid<<=1){
		int w0=qpow(inv==1?gg:ginv,(mod-1)/(mid<<1),mod);
		for(int j=0,R=mid<<1;j<limit;j+=R){
			int wn=1;
			for(int z=0;z<mid;z++,wn=1ll*wn*w0%mod){
				int x=a[j+z],y=a[j+mid+z]*1ll*wn%mod;
				a[j+z]=(x+y)%mod;
				a[j+mid+z]=(x-y+mod)%mod;
			}
		}
	}
	if(inv==-1){
		int linv=qpow(limit,mod-2,mod);
		for(int i=0;i<limit;i++)
		a[i]=1ll*a[i]*linv%mod;
	}
}
int n,m,p;
ll inv1,inv2,inv3;
__int128 mod=1;
int crt(ll a1,ll m1,ll a2,ll m2,ll a3,ll m3,ll p){
	__int128 sum=0;
	sum+=mod/m1*inv1*a1;
	sum+=mod/m2*inv2*a2;
	sum+=mod/m3*inv3*a3;
	return int(sum%mod%p);
}
void work(){
	mod=mod*mod1*mod2*mod3;
	inv1=qpow(mod/mod1%mod1,mod1-2,mod1);
	inv2=qpow(mod/mod2%mod2,mod2-2,mod2);
	inv3=qpow(mod/mod3%mod3,mod3-2,mod3);
	NTT(f[1],1,mod1),NTT(g[1],1,mod1);
	for(int i=0;i<limit;i++)f[1][i]=1ll*f[1][i]*g[1][i]%mod1;
	NTT(f[1],-1,mod1);NTT(f[2],1,mod2),NTT(g[2],1,mod2);
	for(int i=0;i<limit;i++)f[2][i]=1ll*f[2][i]*g[2][i]%mod2; 
	NTT(f[2],-1,mod2);NTT(f[3],1,mod3),NTT(g[3],1,mod3);
	for(int i=0;i<limit;i++)f[3][i]=1ll*f[3][i]*g[3][i]%mod3; 
	NTT(f[3],-1,mod3);
	for(int i=0;i<limit;i++)
	f[0][i]=crt(f[1][i],mod1,f[2][i],mod2,f[3][i],mod3,p);
}
int main(){
	scanf("%d%d%d",&n,&m,&p);
	for(int i=0;i<=n;i++)scanf("%d",&f[1][i]),f[1][i]%=p,f[2][i]=f[3][i]=f[1][i];
	for(int i=0;i<=m;i++)scanf("%d",&g[1][i]),g[1][i]%=p,g[2][i]=g[3][i]=g[1][i];
	while(limit<=(n+m))limit<<=1,l++;
	for(int i=1;i<limit;i++)r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
	work();
	for(int i=0;i<=(n+m);i++)
	printf("%d ",f[0][i]);
	return 0;
} 
2023/6/3 08:44
加载中...