任意模数NTT(MTT)

模板题传送门

问题的简单描述为，求解两个值域为(leq 10^9)的多项式卷积对于(Pleq 10^9)取模的结果

[ ]

问题不能直接用NTT/FFT求解，因为均超过了值域范围(double值域承受不了)

[ ]

Solution1: 3模数NTT

取几个互质的模数分别做一次，然后用中国剩余定理合并

由于值域大，通常需要多次NTT，且中国剩余定理合并常数也不小

实际代码实现也复杂，因此笔者认为不可取

[ ]

Solution2: 拆系数FFT

设(f(x)=sum a_ix^i)

核心：将系数(a_i)分解成(a_i=A_icdot S+C_i,b_i=B_icdot S+D_i)

(其中(Sge sqrt{P})是一个常数，(0 leq A_i,B_i,C_i,D_i<S))

目的是转化后使系数值域变小，double精度可以承受

则最后的答案转化为求解(A_iB_jS^2+(C_iB_j+A_iD_j)S+C_iD_j)

即求解(A_iB_j,C_iB_j,A_iD_j,C_iD_j)，此时值域已经大大缩小

如果直接求解，可以看出要求解4次卷积，需要进行(12)次FFT，不可接受

利用复数的一些性质，有些东西我们可以一起算

构造

(f(x)=sum (A_i,C_i)x^i)

(g(x)=sum(B_i,D_i)x^i)

(f(x)g(x)=sum sum (A_iB_j-C_iD_j, A_iD_j+C_iB_j)x^{i+j})

此时已经得到大部分值了，再构造

(h(x)=sum B_ix^i)

(f(x)h(x)=sum sum (A_iB_j,C_iB_j)x^{i+j})

取一部分即可

最终一共有5次FFT

Tips:

1.这里的负数取整一定要注意，因为C++默认是向0取整，而不是向下取整

2.实际运行表明，这样写用double 很难保证精度，应该要用long double

附:

4次FFT做MTT，但是具体证明比较反人类，而代码非常好看且好写，所以建议直接背板子

Tips: 只要使用了上面提到的最适合FFT的板子，就可以用double，甚至可以开O2

namespace MTT{
	const double PI=acos((double)-1);
	int rev[N];
	struct Cp{
		double x,y;
		Cp(){ ; }
		Cp(double _x,double _y): x(_x),y(_y){ } 
		inline Cp operator + (const Cp &t) const { return (Cp){x+t.x,y+t.y}; }
		inline Cp operator - (const Cp &t) const { return (Cp){x-t.x,y-t.y}; }
		inline Cp operator * (const Cp &t) const { return (Cp){x*t.x-y*t.y,x*t.y+y*t.x}; }
	}A[N],B[N],C[N],w[N/2];

	#define E(x) ll(x+0.5)%P

	void FFT(int n,Cp *a,int f){
		rep(i,0,n-1) if(rev[i]<i) swap(a[i],a[rev[i]]);
		w[0]=Cp(1,0);
		for(reg int i=1;i<n;i<<=1) {
			Cp t=Cp(cos(PI/i),f*sin(PI/i));
			for(reg int j=i-2;j>=0;j-=2) w[j+1]=t*(w[j]=w[j>>1]);
            // 上面提到的最优板子
			for(reg int l=0;l<n;l+=2*i) {
				for(reg int j=l;j<l+i;j++) {
					Cp t=a[j+i]*w[j-l];
					a[j+i]=a[j]-t;
					a[j]=a[j]+t;
				}
			}
		}
		if(f==-1) rep(i,0,n-1) a[i].x/=n,a[i].y/=n;
	}

	void Multiply(int n,int m,int *a,int *b,int *res,int P){
		// [0,n-1]*[0,m-1]->[0,n+m-2]
		int S=(1<<15)-1;

		int R=1,cc=-1;
		while(R<=n+m-1) R<<=1,cc++;
		rep(i,1,R) rev[i]=(rev[i>>1]>>1)|((i&1)<<cc);
		rep(i,0,n-1) A[i]=Cp((a[i]&S),(a[i]>>15));
		rep(i,0,m-1) B[i]=Cp((b[i]&S),(b[i]>>15));
		rep(i,n,R-1) A[i]=Cp(0,0);
		rep(i,m,R-1) B[i]=Cp(0,0);

		FFT(R,A,1),FFT(R,B,1);
		rep(i,0,R-1) {
			int j=(R-i)%R;
			C[i]=Cp((A[i].x+A[j].x)/2,(A[i].y-A[j].y)/2)*B[i];
			B[i]=Cp((A[i].y+A[j].y)/2,(A[j].x-A[i].x)/2)*B[i];
		}
		FFT(R,C,-1),FFT(R,B,-1);

		rep(i,0,n+m-2) {
			ll a=E(C[i].x),b=E(C[i].y),c=E(B[i].x),d=E(B[i].y);
			res[i]=(a+((b+c)<<15)+(d<<30))%P;
		}
	}

	#undef E
}