快速傅里叶变换（FFT）

快速傅里叶变换（FFT）

应用

FFT 在 OI 中多用来解决多项式乘法问题，即：

对于两个多项式：

\[f(x) = \sum\limits_{i=0}^{n-1} a_i x^i \\ g(x) = \sum\limits_{i=0}^{m-1} b_i x^i \]

求一个多项式：

\[h(x) = \sum\limits_{i=0}^{n+m-2} c_i x^i \]

其中：

\[c_i = \sum\limits_{j=0}^{i} a_j b_{i-j} \]

这里超出对应数组下标的部分的值均为 \(0\)。

内容

我们发现如果我们直接相乘的话不太好优化，即系数表示法不是很好优化。我们考虑采用点值表示法。

设有一整数 \(N > n-m+2\)，此时我们取 \(N\) 个互不相同的点 \(x_1,x_2,\cdots,x_N\)，将其分别代入 \(f\) 和 \(g\) 中，分别得到了 \(A_1,A_2,\cdots,A_N\) 和 \(B_1,B_2,\cdots,B_N\)，此时我们很容易的可以知道对于每个 \(1 \leq i \leq N\),\((f \cdot g)(x_i) = h(x_i) = A_i B_i\)，即对应位相乘。在这之后我们可以通过拉格朗日插值得到 \(h\) 的各项系数。

但如果我们随意的代 \(N\) 个整数会发现两部分的时间复杂度都是 \(O(N^2)\)，甚至不如暴力。所以我们考虑优化这个过程。

单位根

对于一个整数 \(n\)，我们在复平面上以原点为圆心做单位圆，以点 \((1,0)\) 为起点将圆 \(n\) 等分，每个等分点都是一个 \(n\) 次单位根。共有 \(n\) 个 \(n\) 次单位根 \(\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}\)，其中 \(\omega_n^0 = 1\)。可以发现对于 \(\omega_n^k\)，\(Re(\omega_n^k) = \cos \frac{2k \pi}{n},Im(\omega_n^k) = \sin \frac{2k \pi}{n}\)，所以 \(\omega_n^k = \cos \frac{2k \pi}{n} + i \sin \frac{2k \pi}{n}\)。

为什么我们要引入单位根？因为单位根有一些很好的性质。

性质1

\[\omega_n^k = \omega_n^{k \bmod n} \]

这个性质比较显然，因为 \(\omega_n^n\) 相当于在圆上转了一圈又回到了始边，所以 \(\omega_n^n = \omega_n^0\)。

性质 2

\[\omega_n^{k_1} \omega_n^{k_2} = \omega_{n}^{k_1 + k_2} \]

这个东西感性理解也比较容易，因为先转 \(k_1\) 次 再转 \(k_2\) 次等价于一次直接转 \(k_1 + k_2\) 次。但是我们也可以理性证明。这里需要用到欧拉公式。

欧拉公式

\[e^{ix} = \cos x + i \sin x \]

证明略（太难了）。

这样我们可以把单位根用另一种方法表示：

\[\omega_n^k = \cos \frac{2k \pi}{n} + i \sin \frac{2k \pi}{n} = e ^{i \frac{2k \pi}{n}} \]

这样单位根做乘法就可以表示为：

\[\omega_n^{k_1} \omega_n^{k_2} = e^{i \frac{2 k_1 \pi}{n}} e^{i \frac{2 k_2 \pi}{n}} = e^{i \frac{2(k_1 + k_2) \pi}{n}} = \omega_n^{k_1 + k_2} \]

性质3（蝶形变换）

\[\omega_n^{k + \frac{n}{2}} = - \omega_n^k \]

易知：

\[\omega_n^{\frac{n}{2}} = \cos \frac{n \pi}{n} + i \sin \frac{n \pi}{n} = \cos \pi + i \sin \pi = -1 \]

基于这样良好的性质，我们选择以单位根来当做每个点去求其对应的值。

分治

为了能够分治到底，数组长度必须是 \(2\) 的整次幂。

我们现在有一多项式：

\[A(x) = \sum\limits_{i=0}^{n-1} a_i x^i \]

我们考虑将其划分为两部分：

\[A(x) = (a_0 x^0 + a_2 x^2 + \cdots + a_{n-2} x^{n-2}) + (a_1 x^1 + a_3 x^3 + \cdots + a_{n-1} x^ {n-1}) \]

我们设：

\[A_0(x) = a_0 x^0 + a_2 x^1 + \cdots + a_{n-2} x^{\frac{n}{2}-1} \\ A_1(x) = a_1 x^0 + a_3 x^1 + \cdots + a_{n-1} x^{\frac{n}{2}-1} \]

则：

\[A(x) = A_0(x^2) + x A_1(x^2) \]

代入 \(n\) 次单位根。对于 \(k < \frac{n}{2}\)：

\[\begin{aligned} A(\omega_n^k) &= A_0(\omega_n^{2k}) + \omega_n^k A_1(\omega_n^{2k}) \\ &= A_0(\omega_{\frac{n}{2}}^k) + \omega_n^k A_1(\omega_{\frac{n}{2}}^k) \end{aligned} \]

同理：

\[\begin{aligned} A(\omega_n^{k + \frac{n}{2}}) &= A_0(\omega_n^{2k + n}) + \omega_n^{k + \frac{n}{2}} A_1(\omega_n^{2k + n}) \\ &= A_0(\omega_n^{2k}) - \omega_n^k A_1(\omega_n^{2k}) \\ &= A_0(\omega_{\frac{n}{2}}^k) - \omega_n^k A_1(\omega_{\frac{n}{2}}^k) \end{aligned} \]

此时我们可以将原问题 \(A\) 分治为两个相同大小的子问题 \(A_0\) 和 \(A_1\) 去解决，这样时间复杂度就问 \(O(n \log n)\)。

快速傅里叶逆变换（IFFT）

刚刚我们只解决了正向求值的问题。我们还需要通过一些手段将点值还原为系数。

我们可以将刚刚的求点值的过程用线性代数的方法来表示为：

\[\begin{bmatrix} 1 &1 &1 & \cdots & 1 \\ 1 &\omega_n^1 &\omega_n^2 & \cdots & \omega_n^{n-1} \\ 1 &\omega_n^2 &\omega_n^4 & \cdots & \omega_n^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{bmatrix} \begin{bmatrix} a_0 \\ a_1\\ a_2\\ \vdots \\ a_{n-1}\\ \end{bmatrix} = \begin{bmatrix} y_0\\ y_1\\ y_2\\ \vdots \\ y_{n-1}\\ \end{bmatrix} \]

我们设 \(X\) 为单位根组成的矩阵。我们在等式两边同时左乘一个 \(X^{-1}\)，可以发现：

\[\begin{bmatrix} a_0 \\ a_1\\ a_2\\ \vdots \\ a_{n-1}\\ \end{bmatrix} = X^{-1} \begin{bmatrix} y_0\\ y_1\\ y_2\\ \vdots \\ y_{n-1}\\ \end{bmatrix} \]

所以我们只要求出 \(X^{-1}\) 就能求出 \(a\)。

通过一些观察可以发现：

\[X^{-1} = \frac{1}{n} \begin{bmatrix} 1 &1 &1 & \cdots & 1 \\ 1 &\omega_n^{-1} &\omega_n^{-2} & \cdots & \omega_n^{-(n-1)} \\ 1 &\omega_n^{-2} &\omega_n^{-4} & \cdots & \omega_n^{-2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \cdots & \omega_n^{-(n-1)^2} \end{bmatrix} \]

证明也比较简单。我们考虑直接表示出 \(X^{-1}X\) 中的每一项：

\[\begin{aligned} (X^{-1}X)_{i,j} &= \sum\limits_{k=0}^{n-1} {X^{-1}}_{i,k} X_{k,j} \\ &= \sum\limits_{k=0}^{n-1} \omega_n^{-ik} \omega_n^{kj} \end{aligned} \]

当 \(i=j\) 时，显然：

\[\sum\limits_{k=0}^{n-1} \omega_n^{-ik} \omega_n^{kj} = \sum\limits_{k=0}^{n-1} 1 = n \]

当 \(i \neq j\) 时：

\[\sum\limits_{k=0}^{n-1} \omega_n^{-ik} \omega_n^{kj} = \sum\limits_{k=0}^{n-1} (\omega_n^{j-i})^k \]

我们设 \(a = j-i\)，则原问题就变为了等比数列求和：

\[\sum\limits_{k=0}^{n-1} \omega_n^{ak} \]

我们设 \(S = \sum\limits_{k=0}^{n-1} \omega_n^{ak}\)，则 \(\omega_n^{a} S = \sum\limits_{k=1}^{n} \omega_n^{ak}\)，错位相减可得：

\[(\omega_n^{a}-1)S = \omega_n^{na} - \omega_n^{0a} = \omega_n^0 - \omega_n^0 = 0 \]

又由于 \(\omega_n^{a}-1 \neq 0\)，所以 \(S = 0\)。

所以当 \(i=j\) 时，\((X^{-1}X)_{i,j} = n\)；当 \(i \neq j\) 时，\((X^{-1}X)_{i,j} = 0\)。

我们发现此时 \(X^{-1}\) 中的每一项都相当于是 \(X\) 中的对应项取倒数，所以我们重新做一次 FFT，只不过在处理单位根的时候给它取个倒数就行了。

易知：

\[\omega_n^{-k} = \cos \frac{-2k \pi}{n} + i \sin \frac{-2k \pi}{n} = \cos \frac{2k \pi}{n} - i \sin \frac{2k \pi}{n} \]

这样我们只需要在 FFT 的时候传个系数 \(1,-1\) 就行了。

所以我们就在 \(O(n \log n)\) 的时间复杂度内求解了 FFT 和 IFFT。

优化

可以发现如果要求多项式乘法，我们会进行 \(2\) 次 FFT 和 \(1\) 次 IFFT。由于运算过程中有大量浮点数和三角函数，再加上递归过程中求单位根的过程不能复用，所以常数会比较爆炸。接下来我们考虑迭代实现。这里分为两步：得到底层序列、向上合并。

得到底层序列

可以发现一个数在进行了递归后到达的位置就是他的二进制翻转（不是取反）。所以这里就可以直接得出每个元素的最终位置。

向上合并

我们设置步长 \(w = 1,2,\cdots,\frac{n}{2}\)，每次我们按照长度为 \(2w\) 来分组，对于每组内的 \(i\)，我们再像递归时一样将 \(i\) 和 \(i+w\) 处的值来处理一下就行了。这样每次迭代单位根的值都是可以复用的，相当于只求了 \(O(n)\) 次单位根的值，效率更高。

这样我们就得到了常数更优秀的 FFT。