快速傅里叶变换 | FFT 初学

FFT

前置

多项式：形如 \(A(x)=\sum\limits_{i=0}^{n-1}a_ix^i\) 的式子，其中 \(n\) 表示项数。

多项式乘法：

\[\begin{aligned}C(x) & = A(x)\cdot B(x) \\ & = \sum\limits_{i=0}^{2n-2}c_ix^i\end{aligned} \]

其中，\(c_i=\sum\limits_{j=0}^ia_jb_{i-j}\)。

多项式表示法：

• 系数表示法：将每一项的系数拿出来写在一起共同组成一个向量，如果缺项就补 \(0\)。一般书写的顺序与数组下标对应。

\[A(x)=(a_0,a_1,a_2,\cdots,a_{n-1}) \]

• 点值表示法：任意的 \(n\) 次多项式都可以由二维平面内任意 \(n+1\) 个点唯一确定，因为一个方程可以解出一个系数。故 \(n\) 次多项式可以用一个包含 \(n+1\) 个二维平面内的点集来表示。

\[A(x)=\{(x_0,y_0),(x_1,y_1),(x_2,y_2),\cdots,(x_{n-1},y_{n-1})\} \]

快速傅里叶变换

对两种表示法分别进行多项式乘法，不难发现，系数表示法的过程是基于乘法分配律的，\(O(n^2)\) 的时间复杂度难以优化；而对于点值表示法来说，如果两个多项式分别是 \(n\) 次和 \(m\) 次的，那么最终的多项式是 \(n+m\) 次的，是需要 \(n + m + 1\) 个点来确定的。那么只需要将两个多项式均用 \(n+m+1\) 个点表示，在两两相乘即可得到目标点值，时间复杂度 \(O(n)\)。

但是，一般题目给定了多项式的系数表示，为了快速实现多项式乘法，我们需要快速实现多项式的系数表示转点值表示，这也就是FFT做的事情。

至此，我们大致了解了多项式乘法的基本过程。

• 求值：将多项式 \(A(x),B(x)\) 通过FFT由系数表示转化成点值表示；
• 点乘：然后将点值表示下的 \(A(x),B(x)\) 通过 \(O(n)\) 的点乘得到结果 \(C(x)\) 的点值表示；
• 插值：将多项式 \(C(x)\) 通过FFT的逆操作还原出 \(C(x)\) 的系数表示。

那么，要如何实现求多项式的一个点值呢？这里我们采用霍纳法则将多项式展开。

\[A(x_0)=a_0+x_0(a_1+x_0(a_2+x_0(a_3+\cdots+x_0a_{n-1})\cdots)) \]

那么我们可以在 \(O(n)\) 的时间内求出单次询问定点 \(x_0\) 的 \(A(x_0)\)。

单位根

定义：在复数域下，我们称满足 \(w^n=1\) 的 \(w\) 为 \(n\) 次单位根，记为 \(w_n\)。

由代数基本定理可知，\(n\) 次单位根有 \(n\) 个。以实轴为横轴，虚轴为纵轴建系，以原点为圆心作单位圆，那么这 \(n\) 个 单位复数根均匀分布在以复平面的原点为圆心的单位半径的圆周上，分别为 \(e^{\tfrac{2\pi ik}{n}}\)，其中 \(k=0\sim n-1\)。

\(e^{ix}=\cos(x)+i\sin(x)\)，借助图像可以更好地理解。

下面给出有关单位根的几则引理。

• 消去引理：若 \(n,k\ge 0,d> 0\)，则有 \(w_{dn}^{dk}=w_n^k\)。

  证明：

  \[w_{dn}^{dk}=(e^{\tfrac{2\pi i}{dn}})^{dk}=(e^{\tfrac{2\pi i}{n}})^k=w_n^k \]

  推论：\(w_n^{\tfrac{n}{2}}=w_2=-1\)
• 折半引理：若 \(n\) 是大于 \(0\) 的偶数，则 \(n\) 个 \(n\) 次单位复数根的平方的集合就是 \(\dfrac{n}{2}\) 次单位复数根的集合，共 \(\dfrac{n}{2}\)。通俗点说，\(n\) 次单位复数根的前后两半平方后是对应相等的。

  证明：

  \[(w_n^{\tfrac{n}{2}+k})^2=w_n^{n+2k}=w_n^n\cdot w_n^{2k}=(w_n^k)^2 \]

  推论：\(w_n^{\tfrac{n}{2}+k}=-w_n^k\)
• 求和引理：对于 \(n\ge 1\) 和满足 \(k\nmid n\) 的非负整数 \(k\)，有
  \[\sum\limits_{j=0}^{n-1}(w_n^k)^j=0 \]

下面讲讲如何在 \(O(n\log n)\) 的时间内实现求值。

我们将 \(n\) 个 \(n\) 次单位复数根代入前面展开后的多项式求点值。对于 \(k\in [0,n)\)，和 \(n-1\) 次的多项式 \(A(x)\)，定义

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

而对于一个多项式 \(A(x)\)，我们把奇数下标和偶数下标分别独立出来，定义成 \(2\) 个子多项式。

\[\begin{cases} A^{[0]}(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\tfrac{n}{2}-1} \\ A^{[1]}(x)=a_1+a_3x+a5x^2+\cdots+a_{n-1}x^{\tfrac{n}{2}-1} \end{cases} \]

这样，将 \(A(x)\) 可以表示成 \(A^{[0]}(x^2)+x\cdot A^{[1]}(x^2)\)。

这显然是可以分治递归求解的，每次次数界折半，时间复杂度 \(O(n\log n)\)。

插值的话，就是反过来做，推到一下式子可得：

\[a_k=\dfrac{1}{n}\sum\limits_{i=0}^{n-1}y_i(w_n^{-k})^i \]

发现求系数的过程与前面求点值的过程类似，也可以在 \(O(n\log n)\) 的时间内求解。

优化

• 蝶形变换

  令 \(n=8\)，那么不难发现递归时 \(a\) 长这样。

  \[[0\quad 1\quad 2\quad 3\quad 4\quad 5\quad 6\quad 7] \]
  \[[0\quad 2\quad 4\quad 6]\ [1\quad 3\quad 5\quad 7] \]
  \[[0\quad 4]\ [2\quad 6]\ [1\quad 5]\ [3\quad 7] \]
  \[[0]\ [4]\ [2]\ [6]\ [1]\ [5]\ [3]\ [7] \]

  记最后的数组为 \(a'\)，则有 \(a'_i=a_{rev(i)}\)，故一开始时可以交换 \(a_i\) 和 \(a_{rev(i)}\) 以减少常数。

• 三次变两次

  发现整个过程先要求两个多项式 \(A(x),B(x)\) 的点值，点乘后又要还原出系数，总共进行了 \(3\) 次FFT。

  把 \(B(x)\) 放到 \(A(x)\) 的虚部上，求出 \(A(x)^2\)，然后把它的虚部 \(\times \tfrac{1}{2}\) 即为答案。

  证明：

  \[(a+bi)^2=(a^2-b^2)+(2abi) \]

---