快速傅里叶变换（FFT）的分治实现
本文作者为 JustinRochester。
快速傅里叶变换（FFT）的分治实现
引入
根据 第一篇 的定义，当我们拥有两个度数分别不超过 \(n,m\) （两者同阶）的多项式 \(\displaystyle F(x)=\sum_{i=0}^n f_ix^i, G(x)=\sum_{i=0}^m g_ix^i\) 时，我们可以计算两者的多项式卷积：
\(\displaystyle H(x)=(F\cdot G)(x)=\sum_{i=0}^{n+m} x^i\sum_{p+q=i}f_pg_q\) 。
在这种情况下，我们计算两个多项式的卷积，即是通过长为 \(n+1, m+1\) 的数组 \(f, g\) 求解长为 \(n+m+1\) 的数组 \(h\) ，使得 \(\displaystyle h_k=\sum_{i+j=k}f_ig_j\) 。
在朴素的算法下，我们求解的复杂度是 \(O(n^2)\) 的。
但根据 第二篇 的知识，如果能求出 \(H(x)\) 上的 \(n+m+1\) 个点，我们通过拉格朗日插值可以唯一地求解一个不超过 \(n+m\) 次的多项式 \(H(x)\) 。
虽然暂时这个方法的复杂度是一致的 \(O(n^2)\) 且常数更大，但它提供了一个更好的思路：我们能不能把多项式用一定数量的点来进行表示？
多项式的点值表示
我们规定，我们用 \(n+1\) 个点（即 \(x\) 点的值与对应的 \(y\) 值）来表示一个度数不超过 \(n\) 的多项式 \(\displaystyle F(x)=\sum_{i=0}^n f_ix^i\) 。这称为 \(F(x)\) 的点值表示法。
朴素上，我们可以在 \(O(n^2)\) 的时间内，通过秦九韶算法求解这些点的值。
当我们需要计算 \(H(x)\) 的乘积时，我们需要用 \(n+m+1\) 个点表示这个不超过 \(n+m\) 次的多项式；而恰好，由于 \(F(x), G(x)\) 的次数都不超过 \(n+m\) ，因此也可以用 \(n+m+1\) 个点进行表示。
而当我们求解 \(H(x)\) 的点值表示结果时，假设我们已经通过上述算法求解出了 \(F(x), G(x)\) 的点值表示 \((x_1, F(x_1)), (x_1, G(x_1)), (x_2, F(x_2)), (x_2, G(x_2)), \cdots\) ，我们可以直接在 \(O(n)\) 的时间内求解 \(H(x)\) 的点值表示法：
\((x_i, H(x_i))=(x_i, F(x_i)\cdot G(x_i))\) ，只要对应采样点的 \(y\) 值相乘即可。
于是我们得到了一个不太一样的多项式加减乘运算的方法：
\(F(x), G(x)\stackrel{\text{秦九韶算法采样}}{\longrightarrow}F(x), G(x)\text{ 的点值表示}\stackrel{\text{点值乘积}}{\longrightarrow}H(x)\text{ 的点值表示}\stackrel{\text{拉格朗日插值}}{\longrightarrow}H(x)\)
这个算法的第一步和第三步是 \(O(n^2)\) 的，而第二步是 \(O(n)\) 的。
快速傅里叶变换
很显然，这个算法的第二步并没有什么优化的空间了。于是，我们考虑如何优化第一步和第三步（下一篇的内容）。
对于第一步的采样，如果我们采样的点之间没有相互的关联，我们只能通过秦九韶算法硬采样。那如果我们代入一些有关联的点，比如单位根呢？
我们假设现在对一个不超过 \(n-1\) 次的多项式 \(\displaystyle F(x)=\sum_{i=0}^{n-1} f_ix^i\) 代入 \(n\) 个 \(n\) 次单位根的各幂次。则 \(F(x)\) 的第 \(i\) 个点值表示法的结果为：\((\omega_n^i, F(\omega_n^i))\) 。
为了方便，我们不妨设 \(n\) 为偶数（不足的可以前补 \(0\)）。我们令 \(\displaystyle F_L(x)=\sum_{i=0}^{{n\over 2}-1}f_{2i}x^i, F_R(x)=\sum_{i=0}^{{n\over 2}-1}f_{2i+1}x^i\) 。
这样一来，我们可以很容易的发现，\(F(x)=F_L(x^2)+xF_R(x^2)\) ，实际上 \(F_L(x),F_R(x)\) 就是将数组 \(F(x)\) 进行奇偶分列。
我们再次代入采样点，会发现 \(\displaystyle F(\omega_n^i)=F_L(\omega_n^{2i})+\omega_n^iF_R(\omega_n^{2i})=F_L(\omega_{n\over 2}^i)+\omega_n^iF_R(\omega_{n\over 2}^i)\)。
我们不妨假设 \(i<{n\over 2}\) ，则再检查 \(\omega_n^{i+{n\over 2}}\) 的结果，就会同样有：
\(\displaystyle F(\omega_n^{i+{n\over 2}})=F_L(\omega_{n\over 2}^{i+{n\over 2}})+\omega_n^{i+{n\over 2}}F_R(\omega_{n\over 2}^{i+{n\over 2}})\) 。
这里同样根据 第四篇 的内容，我们知道：\(\displaystyle \omega_{n\over 2}^{i+{n\over 2}}=\omega_{n\over 2}^i, \omega_n^{i+{n\over 2}}=-\omega_n^i\) 。
于是我们代入化简就会得到：\(\displaystyle F(\omega_n^{i+{n\over 2}})=F_L(\omega_{n\over 2}^i)-\omega_n^iF_R(\omega_{n\over 2}^i)\) 。
对比两式，我们归纳得到：当 \(i<{n\over 2}\) 时：
\(\begin{cases} \begin{aligned} F(\omega_n^i)&=&F_L(\omega_{n\over 2}^i)+\omega_n^iF_R(\omega_{n\over 2}^i)\\ F(\omega_n^{i+{n\over 2}})&=&F_L(\omega_{n\over 2}^i)-\omega_n^iF_R(\omega_{n\over 2}^i)\\ \end{aligned} \end{cases}\)
因此我们发现，求解小于 \(n\) 次的多项式采样点 \(F(\omega_n^i)\) 的问题，就简化为了如何求解两个小于 \({n\over 2}\) 次的多项式采样点 \(F_L(\omega_{n\over 2}^i), F_R(\omega_{n\over 2}^i)\) 的问题。
而后者一旦求解完毕，我们可以在 \(O(n)\) 的时间内，合并得到所有的 \(F(\omega_n^i)\) 。
而对于后者的求解，我们会发现，它的求解结构和原问题一模一样！我们只需要继续递归就可以了。
递归的边界在于多项式减小到不超过 \(1\) 次时，我们已经有 \(F_i(x)=f_i\) 了。此时的递归可以直接返回结果，\(F_i(\omega_1^1)=f_i\) 。
于是时间复杂度为 \(\displaystyle T(n)=2T({n\over 2})+\Theta(n)=\Theta(n\log n)\) 。这个复杂度的证明我们放在后续。
这种实现方法一般叫做 DIT-FFT ，即基于时间抽取的 FFT。另一个常见的实现方法我们会在后续介绍，是 DIF-FFT，即基于频率抽取的 FFT。
快速傅里叶变换的分治实现
当然，在递归求解的过程中，也可能出现次数变为奇数的问题。这个问题可以在递归的过程中，像上文那样补充一项的 \(0\) 使得其变成偶数。但我们一般的实现过程中，为了避免这种麻烦，我们一般都是直接将 \(n\) 补充为 \(2^k\) 。
这个方法的另一个好处是：我们可以预处理所有需要用到的单位根。例如 \(\displaystyle \omega_{n\over 8}^i=\omega_{2^{k-3}}^i=\omega_{2^k}^{8i}=\omega_n^{8i}\) ，我们只需要修改步长即可。
若采用每次升为最近偶数的方法，当多项式次数为 \(17\) 时，我们补充到 \(18\) 并获得了 \(\omega_{18}^0,\cdots, \omega_{18}^{17}\) ；
而第一次递归分解为了两个次数为 \(9\) 的多项式，我们又需要补充到 \(10\) ，然后需要重新计算 \(\omega_{10}^0,\cdots, \omega_{10}^9\) ，因为 \(\omega_{10}^i\) 不能全部在 \(\omega_{18}^j\) 中找到；
继续第二次递归后，得到两个次数不超过 \(5\) 的，又需要补充到 \(6\) ，然后重新计算 \(\omega_6^0,\cdots, \omega_6^5\) ......
依次类推，每次都需要重新计算单位根，加大了运算量。
因此我们将快速傅里叶变换的分治实现总结如下：
输入：
一个度数小于 \(n\) 的多项式 \(\displaystyle F(x)=\sum_{i=0}^{n-1}f_ix^i\) 。
预处理：
将 \(n\) 补充为 \(2^k\) ，预处理 \(\omega_n^0, \cdots, \omega_n^{n-1}\)。
过程：
- 若 \(n=1\) 则直接返回；
- 将 \(F(x)\) 奇偶分列为 \(F_L(x), F_R(x)\) ；
- 递归求解 \(F_L(x), F_R(x)\) 的结果；
- 对 \(i<{n\over 2}\) ，计算 \(\begin{cases} \begin{aligned} F(\omega_n^i)&=&F_L(\omega_{n\over 2}^i)+\omega_n^iF_R(\omega_{n\over 2}^i)\\ F(\omega_n^{i+{n\over 2}})&=&F_L(\omega_{n\over 2}^i)-\omega_n^iF_R(\omega_{n\over 2}^i)\\ \end{aligned} \end{cases}\) ；
- 返回计算结果；
快速傅里叶变换的分治实现-参考代码
int N;
vir w[MAXN];//单位根
inline void FFT(vir *f, int n) {
static vir tmp[MAXN];
if(n==1) return ;
//奇偶分列
for(int i=0; i<n; ++i) tmp[i]=f[i];
vir *fl=f, *fr=f+n/2;
for(int i=0, j=0; i<n; i+=2, ++j) fl[j]=tmp[i];
for(int i=1, j=0; i<n; i+=2, ++j) fr[j]=tmp[i];
//递归求解
FFT(fl, n/2); FFT(fr, n/2);
//O(n) 合并
vir *o=w;
int pace = N/n;
for(int i=0; i<n/2; ++i) {
tmp[i] = fl[i] + *o * fr[i];
tmp[i + n/2] = fl[i] - *o * fr[i];
o += pace;
}
for(int i=0; i<n; ++i) f[i]=tmp[i];
}
inline void mul(vir *f, vir *g, int n, int m) {
for(N=1; N<n+m-1; N<<=1);//求最小的 2 的幂次
vir omega = get_omega(N);
w[0]=1;
for(int i=1; i<N; ++i) w[i]=w[i-1]*omega;
FFT(f, N); FFT(g, N);
for(int i=0; i<N; ++i) f[i]*=g[i];
IFFT(f, N);
}
附录：FFT 复杂度的简单求解
根据 FFT 的复杂度计算式：\(\displaystyle T(n)=2T({n\over 2})+\Theta(n)\) 。
我们不妨假设，当 \(n\leq N\) 时，\(T(n)\leq C_1\) ，以及 \(n\geq N\) 时 \(\displaystyle T(n)\leq 2T({n\over 2})+C_2n\)。
当 \(n=N\cdot 2^k\) 时，有：
\(\begin{aligned} &T(n)\\ =&T(N\cdot 2^k)\\ \leq&2T(N\cdot 2^{k-1})+C_2\cdot N\cdot 2^k\\ \leq&4T(N\cdot 2^{k-2})+2\cdot C_2\cdot N\cdot 2^{k-1}+C_2\cdot N\cdot 2^k\\ \leq&4T(N\cdot 2^{k-2})+2\cdot C_2\cdot N\cdot 2^k\\ \leq&\cdots\\ \leq&2^kT(N)+k\cdot C_2\cdot N\cdot 2^k\\ \leq&2^k\cdot C_1+k\cdot C_2\cdot N\cdot 2^k\\ =&{C_1\over N}\cdot n+C_2(n\cdot \log_2{n\over N})\\ =&C_2\cdot n\log_2 n+({C_1\over N}-C_2\cdot \log_2 N)\cdot n\\ \end{aligned}\)
同理，不妨假设 \(n\geq N_1\) 时 \(\displaystyle T(n)\geq 2T({n\over 2})+C_3n\) 可得到 \(\displaystyle T(n)\geq C_3\cdot n\log_2 n - C_3\cdot \log_2 N\cdot n\)
因此，\(\exists n_0, c_d=C_3-1, c_u=C_2+1\text{ s.t. }n\geq n_0\to c_d\cdot n\log_2 n\leq T(n)\leq c_u\cdot n\log_2 n\) 即 \(T(n)=\Theta(n\log n)\)
否则，我们令 \(n_0\leq N\cdot 2^{k-1}<n<N\cdot 2^k\) ，有：
由 \(n<N\cdot 2^k\) 得 \(T(n)\leq T(N\cdot 2^k)\leq c_u\cdot k\cdot 2^k< c_u\cdot {n\over N}\cdot \log_2 {n\over N}={c_u\over N}\cdot n\log_2 n - {c_u\over N}\log_2 N\cdot n\leq {c_u\over N}\cdot n\log_2 n\leq c_u\cdot n\log_2 n\)
另一方向 \(n>N\cdot 2^{k-1}\) 同理可得 \(T(n)\geq c_d\cdot n\log_2 n\)
于是，此时也有 \(T(n)=\Theta(n\log n)\)
综上，该算法的复杂度是 \(\Theta(n\log n)\) 的。
