任意模数NTT(MTT)

Posted on Edited on In

任意模数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，不可接受。

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

构造：

  1. \(f(x)=\sum (A_i,C_i)x^i\)，

  2. \(g(x)=\sum(B_i,D_i)x^i\)，

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

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

  1. \(h(x)=\sum B_ix^i\)，

  2. \(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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
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
}