FFT (快速傅里叶变换)

· · 个人记录

\text{简介}

FFT 在 OI 中主要是用来加速计算 循环卷积 的一种算法。

即这样一个问题，求序列 c。其中的方括号内条件成立时，取值为 1，否则为 0：

c_k=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}[(i+j)\bmod n=k]a_ib_j

（a,b 的长度都是 n）
如果要计算 普通卷积，实际上就是找一个足够大的 n，使 a,b 非零项的长度之和不大于 n（此时后面的项都是零）；也就是让 (i+j)\bmod n =k 当且仅当 i+j=k。

\text{前置知识}

先来说一下单位根。

然后是单位根的几个性质，下面会用到： $$(1) \ \ \omega_n^k=\omega_{nd}^{kd}$$ 根据第一个定义可直接证明。 $$(2) \ \ \omega_n^{k+n/2}=-\omega_n^k \ (n \bmod 2 = 0)$$ 代入第二个定义式，有 $$\omega_n^{k+n/2}=\cos \left( \frac{2\pi k}{n} +\pi\right)+\text i \sin \left( \frac{2\pi k}{n} +\pi\right)$$ $$=-\cos \frac{2\pi k}{n}- \text i \sin \frac{2\pi k}{n}=- \omega_n^k$$ 然后是一个重点 **单位根反演**。它说的是 $n$ 是否整除 $k$ 这个条件，可以表示为 $n$ 个单位根之和除以 $n$： $$(3) \ \ [n\mid k]= \frac 1n\sum_{i=0}^{n-1}\omega_n^{ki}$$ 证明也很简单，在 $n\mid k$ 不成立时利用等比数列求和有 $$\sum_{i=0}^{n-1} \omega_n^{ki}=\frac{1-\omega_n^{nk}}{1-\omega_n^k}$$ 分子为零，而分母不为零，所以原式为零；而 $n\mid k$ 时，直接按原式算，显然为 $1$，于是得证。 **** $$\text{定义}$$ 为了方便下面的讲述，先定义两个式子。 离散傅里叶变换（DFT）为： $$a'_k=\sum_{i=0}^{n-1}\omega_n^{ki}a_i$$ 其逆变换（IDFT）为： $$a_k=\frac 1n \sum_{i=0}^{n-1}\omega_n^{-ki}a_i'$$ 至于为什么是这样，在后面就有讲解。 **** $$\text{问题的转化}$$ 把**单位根反演**套入循环卷积的式子中： $$c_k=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}[n\mid (i+j-k)]a_ib_j$$ $$c_k=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}\left(\frac 1n \sum_{r=0}^{n-1} \omega_n^{ri}\omega_n^{rj}\omega_n^{-rk}\right)a_ib_j$$ 调整一下求和顺序，就有 $$c_k= \frac 1n\sum_{r=0}^{n-1}\omega_n^{-rk} \left(\sum_{i=0}^{n-1}a_i\omega_n^{ri} \right)\left( \sum_{j=0}^{n-1}b_j\omega_n^{rj}\right)$$ 也就是将序列 $a,b$ 分别进行 DFT 后，对应项乘起来，再做一遍 IDFT。 那我们自然会想，通过一些方法来加速 DFT 和 IDFT 的计算。 **** $$\text{加速 DFT/IDFT}$$ DFT/IDFT 在 $n$ 为 $2$ 的整数幂时，有一种优化算法。 由原序列 $a$ 构造一个多项式： $$A(x)=\sum_{i=0}^{n-1}a_ix^i$$ 那么只要求出 $A(x)$ 在 $\omega_n^0,\omega_n^1,\dots ,\omega_n^{n-1}$ 的点值，实际上就完成了 DFT。 要做 IDFT 也很简单，直接带入共轭复数（或将 $a_1,a_2,\dots,a_{n-1}$ 翻转，对比一下系数可证）后除以 $n$ 即可。 利用分治的思想，对 $a$ 按奇偶分类，得到两个多项式： $$A_1(x)=\sum_{i=0}^{\frac n2-1}a_{2i}x^i$$ $$A_2(x)=\sum_{i=0}^{\frac n2-1}a_{2i+1}x^i$$ 不难发现如下性质 $$A(x)=A_1(x^2)+xA_2(x^2)$$ 那么代入单位根就有 $$A(\omega_n^k)=A_1(\omega_n^{2k})+\omega_n^kA_2(\omega_n^{2k})$$ $$=A_1(\omega_{n/2}^k)+\omega_n^kA_2(\omega_{n/2}^k) $$ $$A(\omega_n^{k + n/2})= A_1(\omega_n^{2k+n})+\omega_n^{k+n/2}A_2(\omega_n^{2k+n})$$ $$=A_1(\omega_{n/2}^k)-\omega_n^kA_2(\omega_{n/2}^k) $$ 推到这一步，很明显就可以分治处理了。先用 FFT 算出 $A_1(x)$ 和 $A_2(x)$ 在 $n/2$ 个单位根处的点值，再用这些值线性求出 $A(x)$ 在 $n$ 个单位根处的点值。代码如下： ```cpp void dft(complex *a,int n){ if(n==1) return; // 边界条件 complex a1[n/2],a2[n/2]; for(int i=0;i!=n/2;++i){ // 按奇偶分类 a1[i] = a[i*2]; a2[i] = a[i*2+1]; } dft(a1,n/2); dft(a2,n/2); for(int k=0;k!=n/2;++k){ complex omega = complex(cos(2*pi*k/n),sin(2*pi*k/n)); // 欧拉公式计算单位根 a[k] = a1[k]+omega*a2[k]; a[k+n/2] = a1[k]-omega*a2[k]; } } void idft(complex *a,int n){ reverse(a+1,a+n); dft(a,n); for(int i=0;i!=n;++i) a[i] /= n; } ``` 根据主定理，可以算出 FFT 的时间复杂度为 $\Theta(n \log n)$。 那么，找一个足够大的、是 $2$ 的整数幂的 $n$，用上面的分治算法，就可以计算 **普通卷积** 了。 那么就有一种比较简单的 **循环卷积** 方式：先做一遍普通的，再把第 $i$ 项加到 $i\bmod n$ 上，再将多出来的部分清零。 那么任意长度的 DFT 有没有快速计算的方法呢？ 答案是肯定的，请看下节。 **** $$\text{Bluestein Algorithm}$$ 还是回到 DFT 的那个式子 $$b_k=\sum_{i=0}^{n-1}\omega_n^{ki}a_i$$ 可以把幂拆开 $$ki= -\frac{(k+i)^2-k^2-i^2}{2}$$ $$b_k= \omega_n^{-\frac {k^2}{2}}\sum_{i=0}^{n-1}\omega_n^{\frac{(i+k)^2}{2}}\omega_n^{-\frac{i^2}{2}}a_i$$ 现在的问题就是如何计算这个和式，设 $$f_i=\omega_n^{\frac{i^2}{2}} \quad g_i=\omega_n^{-\frac{i^2}{2}}a_i$$ $$b_k= \omega_n^{-\frac {k^2}{2}}\sum_{i=0}^{n-1}f_{k+i}g_i$$ 把 $f$ 翻转一下，记为 $f'$，那么有 $$b_k=\omega_n^{-\frac {k^2}{2}}\sum_{i=0}^{n-1}f'_{2(n-1)-k-i}g_i$$ 由于 $g$ 在 $i \geq n$ 的项都为 $0$，所以 $$b_k=\omega_n^{-\frac {k^2}{2}}\sum_{i=0}^{2(n-1)-k}f'_{(2(n-1)-k)-i}g_i$$ 这样就可以用 FFT 来算一遍 **普通卷积**，来实现快速计算任意长度的 DFT。 （注意卷积后取第 $2(n-1)-k$ 项对应 $b_k$） 有了这样一个算法，要计算长度为 $n$ 的序列的 **循环卷积** 快速幂时，只需要 DFT 后将每一项快速幂，再 IDFT 回来即可，时间复杂度 $\Theta(n \log n)$；这比普通卷积倍增要优秀很多。当然如果要在模意义下计算的话，需要保证 $\omega_n$ 能被表示。 **** 可以注意到，上面讲的这些关于 FFT 的知识，虽然都是在**复数域**内进行的，但实际上远不止于此。 只要一个数域内存在 $n$ 阶单位根，也就是存在一个元素 $a$，其 $n$ 次方为乘法的单位元，则在此数域内就能进行长度为 $n$ 的 DFT。比如说常见的模质数 $p$ 意义下运算的数域 $\mathbb F_p$，设 $g$ 为某个原根，且 $n \mid (p-1)$，则 $g^{(p-1)/n}$ 就是一个 $n$ 次单位根。