任意模数NTT

给定两个多项式(f(x),g(x))，求(h(x)=f(x)*g(x))

对(p)取模，(p)不保证可以分解成(a*2^k+1)

三模NTT

由于模数不满足有原根，我们可以找几个有原根的模数，求出结果再(crt)合并一下

考虑卷积后数字最大能达到(p*p*len)，一般是(10^9*10^9*10^5=10^{23})，所以我们选择模数之积应该大于(10^{23})

一般考虑(998244353,1004535809,469762049)，因为原根都是(3)，在(int)范围内且乘积较大

然后跑三个(ntt)（天体常数）得到

[egin{cases}retequiv a_1 (mod p_1)\retequiv a_2 (mod p_2)\retequiv a_3 (mod p_3)\retequiv x (mod p)end{cases} ]

我们现在就是求(x)

由于前三个式子直接合并爆(long long)，我们用一些技巧：

先合并前两个，定义(inv(x,p))为(x)在模(p)意义下的逆元

[retequiv a_1*p_2*inv(p_2,p_1)+a_2*p_1*inv(p_1,p_2) (mod p_1*p_2) ]

记作

[retequiv d (mod m) ]

则

[ret=x*m+d=y*p_3+a_3 ]

[xequiv (a_3-d)*M^{-1} (mod p_3) ]

设(q=(a_3-d)*M^{-1})，那么

[x=k*p_3+q ]

代入(ret)得到

[ret=k*p_1*p_2*p_3+q*M+d ]

然而(ansin [0,p_1*p_2*p_3))，所以(k=0)，(ret)就得出了

#include<bits/stdc++.h>
using namespace std;
namespace red{
#define int long long
#define eps (1e-8)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=4e5+10;
	int mod[3]={469762049,998244353,1004535809};
	int n,m,p;
	int f[N],g[N],pos[N];
	int b[N],ret[N];
	int limit,len;
	inline int fast(int x,int k,int p)
	{
		int ret=1;
		while(k)
		{
			if(k&1) ret=ret*x%p;
			x=x*x%p;
			k>>=1;
		}
		return ret;
	}
	inline int slow(int x,int k,int p)
	{
		int ret=0;
		while(k)
		{
			if(k&1) ret=(ret+x)%p;
			x=(x+x)%p;
			k>>=1;
		}
		return ret;
	}
	struct poly
	{
		int g=3,p,a[N];
		inline void ntt(int limit,int *a,int inv)
		{
			for(int i=0;i<limit;++i)
				if(i<pos[i]) swap(a[i],a[pos[i]]);
			for(int mid=1;mid<limit;mid<<=1)
			{
				int Wn=fast(inv?g:(p+1)/g,(p-1)/(mid<<1),p);
				for(int r=mid<<1,j=0;j<limit;j+=r)
				{
					int w=1;
					for(int k=0;k<mid;++k,w=w*Wn%p)
					{
						int x=a[j+k],y=w*a[j+k+mid]%p;
						a[j+k]=x+y;
						if(a[j+k]>=p) a[j+k]-=p;
						a[j+k+mid]=x-y;
						if(a[j+k+mid]<0) a[j+k+mid]+=p;
					}
				}
			}
			if(inv) return;
			inv=fast(limit,p-2,p);
			for(int i=0;i<limit;++i) a[i]=a[i]*inv%p;
		}
	}fft[3];
	inline int inv(int x,int p)
	{
		return fast(x%p,p-2,p);
	}
	inline void crt()
	{
		int len=n+m;
		int M=mod[0]*mod[1];
		int inv1=inv(mod[1],mod[0]),inv0=inv(mod[0],mod[1]),inv3=inv(M%mod[2],mod[2]);
		int a,b,c,t,k;
		for(int i=0;i<=len;++i)
		{
			a=fft[0].a[i],b=fft[1].a[i],c=fft[2].a[i];
			t=(slow(a*mod[1]%M,inv1,M)+slow(b*mod[0]%M,inv0,M))%M;
			k=((c-t%mod[2])%mod[2]+mod[2])%mod[2]*inv3%mod[2];
			ret[i]=((k%p)*(M%p)%p+t%p)%p;
		}
	}
	inline void main()
	{
		n=read(),m=read(),p=read();
		for(int i=0;i<=n;++i) f[i]=read();
		for(int i=0;i<=m;++i) g[i]=read();
		for(limit=1;limit<=n+m+2;limit<<=1) ++len;
		for(int i=0;i<limit;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
		for(int k=0;k<=2;++k)
		{
			fft[k].p=mod[k];
			for(int i=0;i<=n;++i) fft[k].a[i]=f[i];
			for(int i=0;i<=m;++i) b[i]=g[i];
			for(int i=m+1;i<limit;++i) b[i]=0;
			fft[k].ntt(limit,fft[k].a,1);
			fft[k].ntt(limit,b,1);
			for(int i=0;i<limit;++i)  fft[k].a[i]=fft[k].a[i]*b[i]%mod[k];
			fft[k].ntt(limit,fft[k].a,0);
		}
		crt();
		for(int i=0;i<=n+m;++i) printf("%lld ",ret[i]);
	}
}
signed main()
{
	red::main();
	return 0;
}

MTT

毛爷爷数论变换

对于多项式(P=sumlimits_{i=0}^{n}P_ix^i,Q=sumlimits_{i=0}^{n}Q_ix^i)，求(P*Q)

考虑直接大力(FFT)，发现(10^{23})精度炸飞

考虑降低精度，我们设：

[A=sumlimits_{i=0}^{n}(P_i>>15)x^i,B=sumlimits_{i=0}^{n}(P_i&32767)x^i ]

[C=sumlimits_{i=0}^{n}(Q_i>>15)x^i,D=sumlimits_{i=0}^{n}(Q_i&32767)x^i ]

那么

[P*Q=AC*2^{30}+(AD+BC)*2^{15}+BD ]

需要(8)次(FFT)

我们设(F=A+iC,G=B+iD)

[T1=F*G=AB-CD+i(AD+BC) ]

再设(F^{'}=A-iC)

[T2=F^{'}*G=AB+CD+i(AD-BC) ]

其中

[T1+T2=2(AB+iAD) ]

[T2-T1=2(CD-iBC) ]

我们需要(3)次(DFT),(2)次(IDFT)，一共(5)次(FFT)

其实还能优化！

令

[P(x)=A(x)+iB(x) ]

[Q(x)=A(x)-iB(x) ]

设(P_t)是(P)的(DFT)，(conj(x))表示(x)的共轭复数，(A_i)为(A(x))的(i)次项系数

[P_t(k)=A(omega _n^k)+i B(omega _n^k) ]

[=sumlimits_{j=0}^{n-1}A_jomega _n^{jk}+iB_jomega _n^{jk} ]

[=sumlimits_{j=0}^{n-1}(A_J+iB_j)omega _n^{jk} ]

[Q_t(x)=A(omega _n^k)-i B(omega _n^k) ]

[=sumlimits_{j=0}^{n-1}A_jomega _n^{jk}-iB_jomega _n^{jk} ]

[=sumlimits_{j=0}^{n-1}(A_J-iB_j)omega _n^{jk} ]

[=sumlimits_{j=0}^{n-1}(A_j-i*B_j)(cos(frac{2pi jk}{n})+i*sin(frac{2pi jk}{n})) ]

[=sumlimits_{j=0}^{n-1}(A_jcos(frac{2pi jk}{n})+B_jsin(frac{2pi jk}{n}))+i(A_jsin(frac{2pi jk}{n})-B_jcos(frac{2pi jk}{n})) ]

[=conj(sumlimits_{j=0}^{n-1}(A_jcos(frac{2pi jk}{n})+B_jsin(frac{2pi jk}{n}))-i(A_jsin(frac{2pi jk}{n})-B_jcos(frac{2pi jk}{n}))) ]

[=conj(sumlimits_{j=0}^{n-1}(A_jcos(frac{-2pi jk}{n})-B_jsin(frac{-2pi jk}{n}))-i(A_jsin(frac{-2pi jk}{n})+B_jcos(frac{-2pi jk}{n}))) ]

[=conj(sumlimits_{j=0}^{n-1}(A_j+iB_j)(cos(frac{-2pi jk}{n})+i*sin(frac{-2pi jk}{n}))) ]

[=conj(sumlimits_{j=0}^{n-1}(A_j+iB_j)omega _n^{-jk}) ]

[=conj(sumlimits_{j=0}^{n-1}(A_j+iB_j)omega _n^{(n-j)k}) ]

[=conj(P_t(n-k)) ]

所以我们可以用(P)的点值得到(Q)，也就是说上面的(F^{'})不需要进行(FFT)

注意(n-k)当(k=0)时需要特殊处理，因为是循环卷积，需要移到第(0)项

得到这个结论我们只要(2)次(DFT)，(2)次(IDFT)，一共(4)次(FFT)

碾爆垃圾(9)次(NTT)

注意单位根会乘爆精度，考虑预处理，或者开(long double)

#include<bits/stdc++.h>
using namespace std;
namespace red{
#define int long long
#define eps (1e-8)
	inline int read()
	{
		int x=0;char ch,f=1;
		for(ch=getchar();(ch<'0'||ch>'9')&&ch!='-';ch=getchar());
		if(ch=='-') f=0,ch=getchar();
		while(ch>='0'&&ch<='9'){x=(x<<1)+(x<<3)+ch-'0';ch=getchar();}
		return f?x:-x;
	}
	const int N=1<<18;
	const double pi=acos(-1.0);
	int n,m,p;
	int pos[N];
	int limit,len;
	struct complex
	{
		double x,y;
		complex(double tx=0,double ty=0){x=tx,y=ty;}
		inline complex operator + (const complex &t) const
		{
			return complex(x+t.x,y+t.y);
		}
		inline complex operator - (const complex &t) const
		{
			return complex(x-t.x,y-t.y);
		}
		inline complex operator * (const complex &t) const
		{
			return complex(x*t.x-y*t.y,x*t.y+y*t.x);
		}
		inline complex operator * (const double &t) const
		{
			return complex(x*t,y*t);
		}
		inline complex operator += (const complex &t) const
		{
			return complex(x+t.x,y+t.y);
		}
		inline complex operator ~ () const
		{
			return complex(x,-y);
		}
	}f[N],g[N],a[N],b[N],w[N];
	inline void fft(int limit,complex *a,int inv)
	{
		for(int i=0;i<limit;++i)
			if(i<pos[i]) swap(a[i],a[pos[i]]);
		for(int mid=1;mid<limit;mid<<=1)
		{
			for(int r=mid<<1,j=0;j<limit;j+=r)
			{
				for(int k=0;k<mid;++k)
				{
					complex x=a[j+k],y=w[mid+k]*a[j+k+mid];
					a[j+k]=x+y;
					a[j+k+mid]=x-y;
				}
			}
		}
	}
	inline void main()
	{
		n=read(),m=read(),p=read();
		for(int x,i=0;i<=n;++i) x=read() , f[i].x=x>>15 , f[i].y=x&32767;
		for(int x,i=0;i<=m;++i) x=read() , g[i].x=x>>15 , g[i].y=x&32767;
		for(limit=1;limit<=n+m;limit<<=1) ++len;
		for(int i=0;i<limit;++i) pos[i]=(pos[i>>1]>>1)|((i&1)<<(len-1));
		for(int i=1;i<limit;i<<=1)
    	{
        	w[i]=(complex){1, 0};
        	for(int j=1;j<i;++j)
        	    w[i+j]=((j&31)==1?(complex){cos(pi*j/i), sin(pi*j/i)}:w[i+j-1]*w[i+1]);
    	}
		fft(limit,f,1);fft(limit,g,1);
		for(int i=0;i<limit;++i)
		{
			static complex q,f0,f1,g0,g1;
			q=~f[i?limit-i:0],f0=(f[i]-q)*(complex){0,-0.5},f1=(f[i]+q)*0.5;
			q=~g[i?limit-i:0],g0=(g[i]-q)*(complex){0,-0.5},g1=(g[i]+q)*0.5;
			a[i]=f1*g1,b[i]=f1*g0+f0*g1+f0*g0*(complex){0,1};
		}
		reverse(a+1,a+limit);reverse(b+1,b+limit);
		fft(limit,a,-1),fft(limit,b,-1);
		double k=1.0/limit;
		for(int i=0;i<=n+m;++i)
		{
			printf("%lld ",(((int)(a[i].x*k+.5)%p<<30)+((int)(b[i].x*k+.5)<<15)+(int)(b[i].y*k+.5))%p);
		}
	}
}
signed main()
{
	red::main();
	return 0;
}