FFT&NTT（以及扩展）
FFT&NTT（以及扩展）
预备知识：用于NTT
NTT/FFT其实本质相同，用途是快速求解 多项式乘积
前言
FT: 傅里叶变换:
这是一个工程上的概念，可以简述为：一个周期性的信号波段可以用 若干个正弦曲线 的带权和表示
DFT: 离散傅里叶变换，这是傅里叶变换在离散情况下的变种
FFT: 快速傅里叶变换
NTT: 快速数论变换
谈及核心思想
- 单位根:
构造 \(\omega_n\) 为 \(n\) 阶单位根(不知道 \(\omega_n\) 的值域)，满足性质 \(\omega_n^n=\omega_n^0=1\)，即在幂次上呈现 \(n\) 元循环。
对于 \(2|n\) , \(\omega _n^{\frac{n}{2}}=-1\)，
显然 \(\omega_n\) 满足一个非常简单的性质：折半引理 \(\begin{aligned} \forall 2|i\and 2|n , \omega_n^i=\omega_{\frac{n}{2}}^{\frac{i}{2}}\end{aligned}\)。
- 多项式与点值式的转化：
一个 \(n\) 阶多项式一般的表示就是 \(F(x)=\sum_{i=0}^{n-1} a_ix^i\)。
然而，多项式也可以用 \(n\) 个互不相关的点表示，即 \((x_0,y_0),(x_1,y_1),\cdots,(x_{n-1},y_{n-1})\)，两者可以互相转化。
对于点形式的多项式，在两个多项式求乘积时，同一个 \(x_i\) 对应的 \(y_i\) 可以直接相乘。
FFT / NTT 做多项式卷积的核心过程是：
多项式 \(\stackrel{\text{DFT}}{\Longrightarrow}\) 点值式 \[\Longrightarrow\] 点值式对应相乘 \[\stackrel{\text{IDFT}}{\Longrightarrow}\] 多项式。
而用单位根来构造快速的多项式与点值式的转化，即 \(\text{DFT}/\text{IDFT}\)。
- 分治思想
用于降低多项式与点值式转换的复杂度。
FFT的单位根
\((x,y)\) 指复数 \(i=\sqrt{-1},(x,y)=x+yi\)。
基本运算 \((x,y)+(a,b)=(x+a,y+b),(x,y)\cdot (a,b)=(ax-by,ay+bx)\)。
FFT的单位根是: \(\omega_n\) = \((cos(\frac{2\pi}{n}),sin(\frac{2\pi}{n}))\)。
而 \(\omega _n^i=(cos(\frac{2\pi}{n}\cdot i),sin(\frac{2\pi}{n}\cdot i))\) (展开发现就是三角函数求和公式)。
显然满足单位根的性质。
（实际上可以发现，这个说是点值其实就是信号序列的三角函数表示）。
NTT
相信您已经了解了原根的一些性质， \(\text{NTT}\) 的单位根常用原根构造。
\(\text{NTT}\) 的单位根实际有较大的局限性，对于质数 \(P\) 只能构造出 \(n|P-1,\omega_n=g^{\frac{P-1}{n}}\)。
计算在模意义下就能满足单位根的性质。
通常我们 \(P\) 取 \(998244353\) ， \(2^{23}|(P-1)\) ，它的一个原根是 \(3\)。
实际上，为了满足下面分治需要，构造的模数通常满足 \(P-1=s\cdot 2^t\) 的 \(t\) 较大，这类模数我们常称作 \(\text{NTT}\) 模数。
多项式转点值式
接下来我们考虑如何将多项式转化为点值式。
对于点值式，我们构造的点横坐标为 \(x_i=\omega_n^i\)。
具体目标是对于函数 \(F(x)\) ，求出在 \(x_0,x_1,\cdots ,x_{n-1}\) 上的函数值。
即求出 \(F(x_i)=a_0\omega_n^0+a_1\omega_n^{i}+a_2\omega_n^{2i}+\cdots\)。
接下来就是核心的分治思想，注意，这里的分治是子问题严格等大的。
对于当前问题，分成两部分子问题求解(实际是可以分成多部分的，但是这个是特殊情况暂时不予讨论)，即求解。
令 \(m=\frac{n}{2}\)。
\(\displaystyle 2|i,G(x_i)=a_0\omega_{m}^0+a_2\omega_{m}^{\frac{i}{2}}+a_4\omega_{m}^{\frac{i}{2}\cdot 2}+\cdots\)。
\(\displaystyle 2|i,H(x_i)=a_1\omega_{m}^0+a_3\omega_{m}^{\frac{i}{2}}+a_5\omega_{m}^{\frac{i}{2}\cdot 2}+\cdots\)。
更简洁的描述为。
\(i<m,G(x_i')=a_0\omega_{m}^0+a_2\omega_{m}^{i}+a_4\omega_{m}^{2i}+\cdots\)
\(i<m,H(x_i')=a_1\omega_{m}^0+a_3\omega_{m}^{i}+a_5\omega_{m}^{2i}+\cdots\)
由于 \(G(x'_i),H(x'_i)\) 计算的是 \([0,m-1]\) 项，而求 \(F(x_i)\) 时用到的是 \(0,2,4,\cdots\) 项，实际需要访问 \(G(x^2_i),H(x^2_i)\)
和 \(F(x_i)\) 的式子比较，我们得到合并的式子为：
\[ F(x_i)=G(x^2_i)+x_i H(x^2_i) \] 带入折半引理，实际等价于：
\[ F(x_i)=G(x'_i)+x_i H(x'_i) \] 注意 \(x_i=x'_{i\mod m}\)。
为了保证复杂度，尽量使得每次分治的子问题都分为两部分，这样的复杂度为 \(O(n\log n)\)。
附：实际上，分为 \(d\) 个子问题时，每次合并的复杂度为 \(O(n\cdot d)\) ，因此复杂度为 \(nd \log_d n\)。
保证每次分治为两个严格等大的子问题，可以从一开始就把 \(n\) 扩充为 \(2\) 的幂次：
1 | int N=1; |
附: \(d\) 个子问题时，设子问题答案为 \(G_j(x_i)\) ，则合并的式子为：
\[ \begin{aligned} F(x_i)=\sum_{j=0}^{d-1}x_i^jG_j(x_i^d)=\sum_{i=0}^{d-1}x_i^jG_j(x'_{i\mod \frac{n}{d}})\end{aligned} \]
点值式转多项式
核心性质：单位根反演 \(\sum_{j=0}^{n-1}\omega_n^{ij}= \left\{\begin{aligned} \frac{\omega_n^{in}-1}{\omega_n^i-1}=0 && i\ne 0\\ n && i=0\end{aligned} \right.\)。
设点值式对应 \(y_i\) 的序列为 \(b_i\)，则 \(n\cdot a_i=\sum_{j=0}^{n-1}\omega_n^{-ij} b_j\)。
证明如下：
\(\begin{aligned} \sum_{j=0}^{n-1}\omega_n^{-ij}b_j=\sum_{j=0}^{n-1} \omega_n^{-ij}(\sum_{k=0}^{n-1}a_k\omega_n^{jk})\end{aligned}\)
\(\begin{aligned} \sum_{j=0}^{n-1}\omega_n^{-ij}b_j= \sum_{k=0}^{n-1}a_k\sum_{j=0}^{n-1}\omega_n^{j(k-i)} \end{aligned}\)
由上面的式子，发现只有 \(k-i=0\) 时右边的求和式有值，故上式成立。
因此点值式转多项式直接把系数改为 \(\omega_n^{-i}\) 即可。
Tips:
- 由于单位根的循环特性，溢出会直接溢出到本来的式子里
因此，如果乘法过后的多项式产生了超过 \(>n\) 的项 \(x^i\) ，会溢出到 \(x^{i\mod n}\)。
- 点值式并不是不满足除法，只是除法得到的多项式并不一定是一个 \(n\) 元以内的多项式，除了恰好整除的情况，得到的通常是一个无穷级数的式子，如 \(\begin{aligned} \frac{1}{1-x}=\frac{1-x^{\infty}}{1-x}=\sum_{i=0}^{\infty}x^i\end{aligned}\)。
真正要求除法，通常是求前 \(n\) 项的结果，即需要用到多项式乘法逆。
代码实现与优化
然后我们得到一份优美的代码 （FFT）。
（ Complex 是 C++ 库自带的复数，M_PI 是 C++
自带 \(\pi\) 常量）
1 | void FFT(int n,Complex *a,int f) { |
由于 \((\omega_n)^{\frac{n}{2}}=-1\) ，所以还可以简化为：
1 | Complex w(cos(2*M_PI/n),f*sin(2*M_PI/n)),e(1,0); |
由于用了 double，最后输出要取整。
蝴蝶优化
我们加一点优化，取代递归的分治过程。
可以看到，分治时我们按照 \(i \mod 2\) 分成两组，然后继续分。
这个过程中，实际上我们就是将 \(i\) 的二进制位前后翻转。
所以我们可以暴力处理出 \(i\) 分治底层的位置。
1 | rep(i,0,n-1) { |
当然也是有 \(O(n)\) 处理方法的。
1 | int N=1,c=-1; |
（建议自己模拟一下）。
有了这个翻转数组，我们可以直接从分治底层开始解决整个问题，每次合并操作完全相同。
每次分治问题的大小，依次合并每一个子问题区间即可。
为了在一个数组上完成操作，还需要注意合并顺序。
代码解释 \(i\) ：分治子问题大小为 \(2i\) ， \(l\) :合并区间的左端点为 \(l\) ，右端点为 \(l+2i\) ， \(j\) 枚举合并位置。
1 | void FFT(int n,Complex *a){ |
事实上我们还有更快的写法，就是将 \(\omega_n^i\) 预处理出来。
（注意这个 \(\text{FFT}\) 的预处理很考验double精度，不能每次都直接累乘上去，隔几个就要重新调用依次三角函数）。
当然如果自己写复数会更快。
关于点值式转多项式的优化
由于每次求得点值是 \(\omega_n^{-i}=\omega_n^{n-i}\)。
所以可以直接用 多项式转点值式的函数， 最后把 \([1,n-1]\) 这一段翻转，每个数除掉 \(n\) 即可。
对于加减运算取模的优化
三目运算：
1 | a+=b,a=a>=P?a-P:a; |
逻辑运算优化（原理是逻辑预算会在第一个确定表达式值的位置停下）。
1 | a+=b,((a>=P)&&(a-=P)); |
关于系数预处理优化（以NTT为例）
带入上面已经提到的优化，无预处理系数的 \(\text{NTT}\) 大概是这样的：
1 | ll qpow(ll x,ll k=P-2){ |
一种简单的预处理是，每次对于每个分治大小，预处理依次系数：
1 | ll qpow(ll x,ll k=P-2){ |
另一种是在一开始就把所有的系数用一个数组存下来，具体过程可以描述为：
对于每个分治长度 \(n\) ，我们只需要访问 \(\omega_n^{0},\omega_n^{1},\cdots,\omega_n^{\frac{n}{2}-1}\)。
那么对于分治长度 \(n\) ，我们在 \(w\) 数组的第 \(\frac{n}{2}\) ~ \(n-1\) 项依次存储这些值。
优化:我们只需要对于最大的分治长度处理，剩下的部分发现可以直接用折半引理访问得到。
1 | ll qpow(ll x,ll k=P-2){ |
三份代码在 duck.ac 上的评测结果表明，不预处理系数将近慢一倍。
单组数据来看，预处理系数会慢一点；多组来看，预处理系数会快。实际差距不大，都可以使用。
但是在某些层面来说，下面这份板子才是最好的（适用NTT,FFT且精度较高），不需要预处理。
1 | ll qpow(ll x,ll k=P-2){ |
拓展
1. 分治+NTT
常用于处理多个计数背包的快速合并 （实际无权值01背包也是可以的）。
我们可以用 NTT \(n\log n\) 合并两个大小为 \(n\) 的背包。
分治时，每次合并两个分治子问题，总共的时间就是 \(\sum size\log n\)。
每个背包的 \(size\) 会被计算 \(\log n\) 次，所以总共复杂度是 \(n \log ^2 n\)。
2. CDQ+NTT
对于形如 \(dp_i=\sum_{j=0}^{i-1}dp_jg_{i-j}\) 的 \(dp\) 转移（就是dp转移与差值有关）。
由于求 \(dp_i\) 时，需要保证 \(dp_0,dp_1,\cdots,dp_{i-1}\) 才能卷积，这个限制，我们可以用CDQ分治解决。
对于当前分治区间 \([L,R]\)。
依次考虑 \([L,mid]\) 内部转移， \([L,mid]\) 向 \([mid+1,R]\) 的转移(用FFT/NTT解决)， \([mid+1,R]\) 内部转移。
算法流程：
1 | void Solve(l,r){ |
3. MTT(任意模数NTT)
4.\(n\) 元点值式
练习建议：
9.图上 \(dp\) ：
联通图个数：BZOJ-3456 题解
森林数量和带限制森林数量：HDU - 5279 题解
10.点分治+FFT：CodeChef-PRIMEDST 题解
更多应用和优化参见毛啸2016论文
（如：两次FFT做卷积，4次FFT做MTT。。。）。