任意模數NTT(MTT)

模板題傳送門

問題的簡單描述為，求解兩個值域為\(\leq 10^9\)的多項式卷積對於\(P\leq 10^9\)取模的結果

\[\ \]

問題不能直接用NTT/FFT求解，因為均超過了值域範圍(double值域承受不了)

\[\ \]

Solution1: 3模數NTT

取幾個互質的模數分別做一次，然後用中國剩餘定理合併

由於值域大，通常需要多次NTT，且中國剩餘定理合併常數也不小

實際程式碼實現也複雜，因此筆者認為不可取

\[\ \]

Solution2: 拆係數FFT

設\(f(x)=\sum a_ix^i\)

核心：將係數\(a_i\)分解成\(a_i=A_i\cdot S+C_i,b_i=B_i\cdot S+D_i\)

(其中\(S\ge \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
}