FFT的模数很操蛋的时候怎么办

一个说是简单基础，但是很烦人的问题。我个人认为要想清楚还是不太容易的。。至少在我脑子烧坏的情况下。。

还有这篇文章其实和FFT没有什么关系，嗯嗯。

挺没意思的一个东西。

按照picks博客上的做法，取三个可以用来做NTT的模数，然后中国剩余定理合并一下就是了。可能要用int128或者手写高精。

“当然，FFT之后常数乘上三，估计也没有什么题能跑过去了…… ”

首先是，模数$p\le 10^6$的情况，可以直接取一个大质数，比如说$135435435\times 2^{25}+1$(原根为$7$)，每算一次就模$p$，就可以了。

一个细节是ll乘ll模ll(就当ll是long long的缩写吧)。我们希望做到$O(1)$。。这就需要用到double的力量辣！

求$xy\text{ mod } M$。我们用double可以求出$\lfloor xy/M \rfloor$，把它下取整设为$D$，用ll求出$xy-DM$就可以了。其中$xy$是个ll值，实际上等于$xy$模$2^{63}$的值；$DM$也是如此。但是都爆掉之后它们的差是不会变的。

代码

const ll mod = 135435435ll << 25 | 1;
void MOD(ll &x){if (x >= mod) x -= mod; if (x < 0) x += mod;}
ll Tm(ll x, ll y){
  long double D = 1.0 * x * y / (1.0 * mod);
  ll r = x * y - ll(D) * mod; MOD(r);
  return r;
}

然后是$p\le 10^9$的情况。首先吐槽一下如果有int128的话可以直接取更好的质数啊。。比如说$15\times 2^{112}+1$，原根为$11$。。那么如果没有int128的话就只能强上crt了。。当然不需要手写高精。。

取两个模数。。$p_1,p_2$，保证$p_1\times p_2>np^2$。我们可以得到答案模$p_1,p_2$的值分别是$r_1,r_2$。用扩欧可得整数$a_1,a_2$，满足$a_1p_1+a_2p_2=1$。我们知道答案就是$(r_1a_2p_2+r_2a_1p_1)\text{ mod }(p_1p_2)$。换句话说就是$W=(r_1a_2 \text{ mod } p_1)\times p_2+(r_2a_1\text{ mod }p_2)\times p_1$，答案等于$W\text{ mod }p_1p_2$。$W\text{ mod }p$的值可以用上面的方法算出来。

我们判断$W$和$p_1p_2$的大小关系，如果$W\ge p_1p_2$，那么答案是$W\text{ mod }p-(p_1p_2)\text{ mod }p$模$p$；否则答案就是$W\text{ mod }p$。

怎么判断呢。。又是double登场了：令$a=W,b=p_1p_2$，这两个数都是double的表示范围内。。直接比较。。

不过如果$a,b$太过接近，那么就比较$W\text{ mod }2^{63},p_1p_2\text{ mod }2^{63}$的大小关系就好了。

代码中的p代表模数，Tm(x,y)代表(x*y)%p，r1,r2,p1,p2,m1,m2,r分别代表上面的$r_1,r_2,a_1,a_2,p_1,p_2$和答案。

F1.MOD(tmp = F1.Tm(r1, p2)); F2.MOD(t2 = F2.Tm(r2, p1));
r = Tm(tmp, m2 % p) + Tm(t2 % p, m1);
a = double(tmp) * m2 + double(t2) * m1;
b = double(m1) * m2;
if (a > b || (a > b * (1 - 1e-12) && tmp * m2 + t2 * m1 > m1 * m2))
  MOD(r -= Tm(m1, m2 % p));
else MOD(r);

upd：$p\le 10^6$的时候，实数fft精度是足够的?