数学算法介绍
主要介绍一些数学相关的算法。
多项式
DFT:离散傅里叶变换，方式有两种， FFT 和 NTT。
IDFT:逆离散傅里叶变换。
FFT 多项式乘法
最开始学这玩意的时候感觉非常迷，后面数学水平上去了其实也不难。
多项式有两种表示方法，一种是系数法，另一种是点值法，总所周知 \(n\) 个不同点唯一确定一个 \(n-1\) 次多项式。
原理：多项式的两个系数表达式相乘是 \(O(n^2)\) 的，但是其点值表达式相乘却是 \(O(n)\)
的，所以考虑将系数表达式转成点值表达式然后相乘。
事实上，点值表达式和系数表达式的互相转化，如果点值取特殊点，可以做到
\(O(n\log
n)\) ，即使是任意点，即多项式多点求值和多项式多点插值也可以做到
\(O(n \log^2n)\)
设最终多项式次数为 \(n-1\) ，我们进行多项式乘法时选择的点值叫单位根，即
\(x^n=1\) 在复数域上的所有根。
这玩意有一些性质，不过我们得先把次数变为 \(n=2^k\) 形式。
无法想象发明这个东西的人是怎么想到的，可能这就是被记在历史书上的人的水平。
以下内容如果将坐标系视为极坐标系会更好理解 $$ \[\begin{align}
W_n^i &= -W_n^{i+\frac{n}{2}}\\
W_{\frac{n}{2}}^{i} &= W_{n}^{2i}\\
\end{align}\] $$ 考虑这样一个问题，对于一个多项式 \(a_0+a_1x+a_2x^2 \cdots a_{n-1}x^{n-1}\)
，我们需要同时求出它在 \(W_n^{0},W_n^{1}\cdots
W_n^{n-1}\) 处的取值。发现由于第一个性质，貌似可以偷个懒，因为后
\(\frac{n}{2}\) 个 \(W_i\) 就是前 \(\frac{n}{2}\) 个数的相反数。
相反数的性质，奇数次幂的符号要改变。考虑对系数按奇偶性分类，式子变成了这个样子。
\((a_0+a_2x^2+\cdots
+a_{n-2}x^{n-2})+x(a_1+a_3x^2+\cdots + a_{n-1}x^{n-2})\)
然后考虑前后两个部分，都是一个多项式，需要对他们各自求 \(x=W_\frac{2}{n}^0,W_\frac{2}{n}^1\cdots
W_\frac{2}{n}^{\frac{n}{2}-1}\) 处的取值，本质上是求 \(x^2=W_{\frac{2}{n}}^0,W_{\frac{2}{n}}^2\cdots
W_{\frac{2}{n}}^{n-2}\)
处的取值，结合第二个性质，woc，就是两个子问题，解决之后就可以 \(O(n)\) 得到原问题的解，边界显然是 \(n=1\) 。
复杂度 \(T(n)=2T(\frac{n}{2}) + O(n) =
O(n\log n)\)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 void FFT (int now,com *a,int op) if (now==0 ) return ; com a1[1 <<now],a2[1 <<now]; for (i=0 ;i<1 <<now;i+=2 ) { a1[i/2 ]=a[i]; a2[i/2 ]=a[i+1 ]; } FFT (now-1 ,a1,op); FFT (now-1 ,a2,op); com w0=(com{cos (2.0 *Pi/(1 <<now)),op*sin (2.0 *Pi/(1 <<now))}),w=(com){1 ,0 }; for (i=0 ;i<1 <<now-1 ;i++,w=w*w0) { a[i]=a1[i]+w*a2[i]; a[i+(1 <<(now-1 ))]=a1[i]-w*a2[i]; } }
但是如果递归的话，常数会比较拉跨。因为递归必然需要复制数组重新弄成一个下标
\(1-n\)
的问题，无论用什么办法解决，你的高速缓存都会表示意见很大，所以考虑迭代写法。
本质上递归是一层一层合并了两个数组，那么能不能直接模拟这个合并的过程呢，答案是可以的。
观察发现本质上是将下标二进制 reverse
之后逐层合并的，我们也这么做就行。
求 reverse 可以 \(O(n)\) ，如下（如果你不了解运算顺序，请老老实实打括号）
这个原理很简单，不看最后一位，其它位先
reverse，然后处理一下最后一位就行。
1 for (int i=1 ;i<1 <<lim;i++)res[i]=res[i>>1 ]>>1 |(i&1 )<<1 -lim;
下面是迭代写法代码，本质是模拟了递归合并的过程。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 void FFT (com a[],com b[],int op) for (i=0 ;i<1 <<maxn;i++) if (r[i]>i)swap (a[r[i]],a[i]),swap (b[r[i]],b[i]); for (i=1 ;i<=maxn;i++){ for (j=0 ;j<1 <<maxn;j+=1 <<i){ com w0=(com){cos (2.0 *Pi/(1 <<i)),op*sin (2.0 *Pi/(1 <<i))},w=(com){1 ,0 }; for (k=0 ;k<1 <<i-1 ;k++,w=w*w0){ com x=a[j+k],y=a[j+k+(1 <<i-1 )]*w; a[j+k]=x+y; a[j+k+(1 <<i-1 )]=x-y; x=b[j+k],y=b[j+k+(1 <<i-1 )]*w; b[j+k]=x+y; b[j+k+(1 <<i-1 )]=x-y; } } } }
搞完了系数转点值，接下来是点值转系数。
前人告诉我们只需要将单位根改为 \(W_n^{0},W_n^{-1}\cdots
W_n^{-n+1}\) ，再做一遍系数转点值的过程就可以得到系数，但是系数会变成原来的
\(n\) 倍，除掉就行。
给出简要证明，\(i\) 次项的系数为
\(a_i\) ，转一次点值
之后变为 \(b_i\) ，再做一次变成 \(c_i\) \[
\begin{eqnarray}
c_x &=& \sum\limits_{i=0}^{n-1}b_iW_n^{-ix} \\
&=& \sum\limits_{i=0}^{n-1}
W_n^{-ix}\sum\limits_{j=0}^{n-1}a_jW_n^{ij}\\
&=& \sum\limits_{i=0}^{n-1}
\sum\limits_{j=0}^{n-1}a_jW_n^{i(j-x)}\\
\end{eqnarray}
\] 对于 \(j=x\) ，贡献显然为
\(\sum\limits_{i=0}^{n-1}a_xW_n^{i\times0} =
na_x\)
对于 \(j\neq x\) 贡献为 \(a_j
\sum\limits_{i=0}^{n-1}(W_n^{j-x})^{i}\)
对这个式子的求和用等比数列求和公式有贡献为 \(\dfrac{W_n^{n(j-x)}-1}{W_n^{j-x}-1}\)
显然分子为 \(0\) ，分母不为 \(0\) ，所以贡献是 \(0\) ，所以结果就是 \(na_x\)
搞定。
NTT 多项式乘法
FFT 多项式乘法是由缺陷的，由于浮点数精度和运算速度问题，FFT
可能并不能很好的解决一些问题，所以引入了 NTT，NTT
从有限整数域中找到了这样一组具有同样优秀性质的 \(W_n\) ，即 \(g\) ，也就是原根。
原根的内容可以参考数学证明总结中的介绍。
注意，和 FFT 一样，NTT 也需要严格的按照 \(2^k\) 取次数，因为我们利用了 \(W_n^{2i} = W_{\frac{n}{2}}^{i}\)
这一重要性质
所以能取出较大的 \(2^k\)
作为阶的质数才可以作为 NTT 的模数，常见的 NTT 模数是 \(998244353=2^{23}\times 7\times 17 +1\)
，我们可以取它的原根 \(g=3\)
作为基本单位根带入，实际上如果要找到一个应用于 \(n\) 的单位根 \(W_n\) ，需要取 \(W_n=g^{\frac{p-1}{n}}\) 。这样它就满足了我们在
FFT 证明中用到的一切性质。
然后照着 FFT 打一遍就行，只是基本运算这些换为模 \(p\) 意义下的运算就行。
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 void NTT (int *a,int type) for (i=0 ;i<1 <<s;i++) if (rk[i]>i)swap (a[rk[i]],a[i]); for (int len=1 ;len<=s;len++) { int w=1 ,wn=quick (g,mod-1 >>len); if (type==-1 )wn=quick (wn,mod-2 ); for (j=0 ;j+(1 <<len)<=1 <<s;j+=1 <<len,w=1 ) { for (k=j;k<j+(1 <<len-1 );k++,w=1ll *w*wn%mod) { int x=a[k],y=a[k+(1 <<len-1 )]; a[k]=(x+1ll *w*y%mod)%mod,a[k+(1 <<len-1 )]=(x-1ll *w*y%mod)%mod; } } } }
其它多项式乘法和一些优化
FFT 三次变两次
FFT 三次变两次，把 \(b\) 扔到 \(a\) 的虚部去，变成了 \(A(x) = (a_0+b_0i) + (a_1+b_1i)x+ \cdots +
(a_{n-1}+b_{n-1}i)x^{n-1}\)
然后求 \(A^2(x)\) ，得到的系数表达式的虚部就是 \(2ab\) 。
证明显然，会比 NTT 略快。
任意模数 NTT
NTT
解决的是模意义下的乘法问题，当然，如果确保值域不超出模数 ，那么模意义下的结果就是正确结果。
如果模数不是 NTT 模数，就需要写任意模数 NTT。
实现方式有两种，其一是拆系数 FFT，其二是多模数 NTT 后用 CRT
合并。
多模数 NTT
一般选取 998244353 1004535809 469726049 作为 NTT
模数，它们的原根都是 3。分别进行 NTT
后，可以得到在模它们意义下的结果，用中国剩余定理合并，得到一个模约 \(5\times 10^{26}\)
的数的结果，一般来说任意的模数不会超过 \(10^9\) ，所以乘起来的结果不超过 \(10^{24}\) ，直接用合并的结果就行。
提一下实现细节。
一种方式是写个 Num 类，里面放三个 int，重写一下加减乘除，正常做。
另一种方式是写个 Poly 类，封装个乘法，用 for 做三次，都比较阳间。
用 CRT 时记得 __int128。或者用一个科技
一共三个方程，考虑这样合并。
其中 \(mod_1^{-1}\) 表示 \(mod_1\) 在 \(mod_2\) 意义下的逆元。 \[
x = x_1 \pmod {mod1} \\
x = x_2 \pmod {mod2} \\
k_1 mod_1 + x_1 = x_2 \pmod{mod2}\\
x = (x_2-x_1)mod_1^{-1}mod_1 + x_1 \pmod {mod_2*mod_1}
\]
发现这样第二次合并的时候可以直接在模需要的模数 \(p\) 下进行计算。
注意，务必考虑好取模的顺序，计算第二个时不能先对 p
取模，因为可能结果还需要先模上 \(mod_1*mod_2*mod_3\) ，
怎样避免 __int128 ?
把 \((x_2-x_1)mod_1^{-1}\) 对 \(mod_2\)
取模即可，注意处理正负号，保证任意时刻 结果在
[0,对应模数) 内。
拆系数 FFT
FFT 精度不高，一般来说 double 跑个 \(10^{14}\) 问题不大（最后除掉 \(n\)
以后的结果，也就是乘出来的多项式的系数，举个例子就是如果跑 \(n=10^4\) ，那需要保证 \(a_i,b_i\leq 10^5\) ），但是跑 \(10^{24}\) 次方， FFT 表示我做不到，用 long
double 也不行。
所以考虑拆系数，把每个数写成 \(a_1\times2^{15} + a_2\) 的形式。
然后乘法就变成了 \(a_1b_1*2^{30}+(a_1b_2+a_2b_1)*2^{15}+a_2b_2\) 。
设新的四个多项式为 \(A_1,A_2,B_1,B_2\)
一共 8 次 FFT，但可以优化，考虑复多项式 \[
P=B_1+B_2i\\
T_1=P\times A_1 \\
T_2=P\times A_2
\]
发现需要的系数在 \(T_1,T_2\)
的实部和虚部，对这俩做 IDFT，收工。
一共 5 次 DFT。
也可以这样:
先考虑 \[
P=A_1+A_2i\\
P'=A_1-A_2i\\
Q=B_1+B_2i
\]
再考虑： \[
T_1=PQ=(A_1B_1-A_2B_2)+(A_1B_2+A_2B_1)i\\
T_2=P'Q=(A_1B_1+A_2B_2)+(A_1B_2-A_2B_1)i
\] 最后： \[
T_1+T_2=2(A_1B_1+A_1B_2i)\\
T_1-T_2=-2(A_2B_2-A_2B_1i)
\]
弄出 \(T_1,T_2\)
的系数表示法就行。其实复系数多项式乘法和普通系数多项式乘法没有什么区别。
这种 5 次的思路还有其它方式，不一一介绍了。
听说有种方法可以把两次实数域多项式的 DFT 和 IDFT
整成一次复数域的，但我觉得没啥必要会。
DFT 思路是构造一个多项式 \(P=A+Bi\) ，对其做 DFT，然后线性求 \(Q=A-Bi\) 的 DFT 结果。
IDFT 思路不会。
更推荐第一种写法，不会出现精度问题，如果效率不够再考虑改成第二种