快速傅里叶变换——加速多项式乘法的神奇魔法

0 为什么要用快速傅里叶变换（FFT）

定义以下两个多项式： $$ f(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_{n-1} x^{n-1} \\ g(x) = b_0 + b_1 x + b_2 x^2 + \dots + b_{m-1} x^{m-1} $$ 想要求 $(f \cdot g)(x)$ 首先我们会想到的是最暴力的做法： $$ (f \cdot g)(x) = \sum^{n-1}_{i=0} \sum^{m-1}_{j=0} a_i b_j x^{i + j} $$ 很显然这种方法非常直观且容易实现，但是缺点也很明显，这个方法的时间复杂度为 $O(nm)$ ，当 $nm \geq {10}^8$ 时，即使是目前最强的计算机也无法快速地得出答案，而这种情况下两个多项式的平均长度也才 ${10}^4$ 的数量级。这显然不能满足很多情况的需求，因为这仅仅只是一次乘法运算，而真正的应用场景不可能只有简单的几个乘法，因此我们需要一种更快的方法来帮助我们解决这个问题。

1 多项式的系数表示法和点值表示法

对于一个多项式 $$ f(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_{n-1} x^{n-1} $$ 我们显然知道它有 $n$ 个系数 ${ a_0, a_1, a_2, \dots, a_{n-1} }$ ，有这 $n$ 个系数就可以唯一确定这个多项式，那么这时我们就来回顾一下如何确定二次函数的解析式：

对于一个二次函数: $f(x) = a x^2 + b x + c$ ，初中我们就学过了可以利用待定系数法求出 $a, b, c$ 的值，即利用三个在函数图像上且 $y$ 不相等的点代入就能得到三元二次方程组： $$ \left\{ \begin{aligned} a x^2_1 + b x_1 + c = y_1 \\ a x^2_2 + b x_2 + c = y_2 \\ a x^2_3 + b x_3 + c = y_3 \end{aligned} \right. ， y_1 \neq y_2 \neq y_3 $$ 显然这个方程组存在唯一的一组解，只要对这个方程组进行求解我们就可以唯一确定这个多项式。

对于 $3$ 项的多项式我们可以这么做，显然对于有 $n$ 项的多项式我们也可以这么做。
找出 $n$ 个不同的点 $(x_i, f(x_i))$ 代入即可唯一确定这个多项式，即： $$ \left\{ \begin{aligned} a_0 + a_1 x_0 + a_2 x^2_0 + \dots + a_{n-1} x^{n-1}_0 & = f(x_0) \\
a_0 + a_1 x_1 + a_2 x^2_1 + \dots + a_{n-1} x^{n-1}_1 & = f(x_1) \\
a_0 + a_1 x_2 + a_2 x^2_2 + \dots + a_{n-1} x^{n-1}_2 & = f(x_2) \\
\vdots \\
a_0 + a_1 x_{n-1} + a_2 x^2_{n-1} + \dots + a_{n-1} x^{n-1}_{n-1} & = f(x_{n-1}) \end{aligned} \right. $$ 同理，我们可以用这个方法求出 $(f \cdot g)(x)$，而显然 $(f \cdot g)(x_i) = f(x_i) \cdot g(x_i)$ ，此时我们就可以用 $O(1)$ 的时间复杂度来解决这个问题，从而达到线性处理问题的效果。

2 多项式的系数表示法和点值表示法之间的转化

知道了存在线性的处理方法，接下来就该考虑如何让多项式的表示在两种形式之间转化，同时考虑要找的点值对 $(x_i, f(x_i))$ ，如何找到合适的 $x_i$ 方便我们计算呢？

复数的单位根

这里先插入一点小知识。
设 $\omega^n_n = 1$， 则这里的 $\omega_n = cos(\frac{2 \pi}{n}) + i sin(\frac{2 \pi}{n})$ 。
由欧拉公式 $$ e^{ix} = cos(x) + isin(x) $$ 所以 $\omega^k_n = cos(\frac{2k \pi}{n}) + i sin(\frac{2k \pi}{n})$ ，特别地：$\omega^0_n = \omega^n_n = 1$
经过简单推导就可以得出下面三个结论： $$ \begin{aligned} \omega^{rk}_{rn} & = cos(\frac{2rk \pi}{rn}) + i sin(\frac{2rk \pi}{rn}) \\ & = cos(\frac{2k \pi}{n}) + i sin(\frac{2k \pi}{n}) \\ & = \omega^k_n \end{aligned} , r \in {\mathbb{N}}^+ , k \in {\mathbb{N}}^+ \\ \begin{aligned} \omega^{k+\frac{n}{2}}_n & = cos((k + \frac{n}{2}) \frac{2 \pi}{n}) + i sin((k + \frac{n}{2}) \frac{2 \pi}{n}) \\ & = cos(\frac{2k \pi}{n} + \pi) + i sin(\frac{2k \pi}{n} + \pi) \\ & = - cos(\frac{2k \pi}{n}) - i sin(\frac{2k \pi}{n}) \\ & = - (cos(\frac{2k \pi}{n}) + i sin(\frac{2k \pi}{n})) \\ & = - \omega^k_n \end{aligned} , k \in {\mathbb{N}}^+ \\ \begin{aligned} \bar{\omega^k_n} & = cos(\frac{2k \pi}{n}) - i sin(\frac{2k \pi}{n})\\ & = cos(2 \pi - \frac{2k \pi}{n}) + i sin(2 \pi - \frac{2k \pi}{n})\\ & = cos((n - k) \frac{2 \pi}{n}) + i sin((n - k) \frac{2 \pi}{n}) \\ & = \omega^{n-k}_n \end{aligned} , k \in {\mathbb{N}}^+ $$

多项式项项的分组

对于多项式（ $n$ 为偶数） $$ f(x) = a_0 + a_1 x + a_2 x^2 + \dots + a_{n - 1} x^{n - 1} $$ 我们按照指数的奇偶性对其分组 $$ \begin{aligned} f(x) & = a_0 + a_2 x^2 \dots + a_{n - 2} x^{n - 2} + a_1 x + a_3 x^3 \dots + a_{n - 1} x^{n - 1} \\ & = a_0 + a_2 x^2 \dots + a_{n - 2} x^{n - 2} + x (a_1 + a_3 x^2 \dots + a_{n - 1} x^{n - 2}) \end{aligned} $$ 设 $$ \begin{aligned} f_1(x) & = a_0 + a_ x \dots + a_{n-2} x^{\frac{n}{2} - 1} \\ f_2(x) & = a_1 + a_3 x \dots + a_{n-1} x^{\frac{n}{2} - 1} \end{aligned} $$ 有 $$ f(x)=f_1(x^2) + x f_2(x^2) $$

系数表示法转化为点值表示法

当 $x = \omega^k_{2n}, k < n, k \in \mathbb{N}$ 时 $$ \begin{aligned} f(\omega^k_{2n}) & = f_1(\omega^{2k}_{2n}) + \omega^k_{2n} f_2(\omega^{2k}_{2n}) \\
& = f_1(\omega^k_n) + \omega^{\frac{k}{2}}{n} f_2(\omega^k_n) \end{aligned} \tag{1} $$ 当 $x = \omega^{k + n}_{2n}, k < n, k \in \mathbb{N}$ 时
$$
\begin{aligned}
f(\omega^{k + n}_{2n}) & = f_1(\omega^{k + n}_n) + \omega^{k + n}_{2n} f_2(\omega^{k + n}_n) \\
& = f_1(\omega^k_n \cdot \omega^n_n) - \omega^k_{2n} f_2(\omega^k_n \cdot \omega^n_n) \\ & = f_1(\omega^k_n) - \omega^k_{2n} f_2(\omega^k_n) \ \end{aligned} \tag{2} $$ 观察 $(1)(2)$ 两式，我们发现知道计算出 $f_1(\omega^k_n), f_2(\omega^k_n)$，就能计算出 $f(\omega^k_{2n}), f(\omega^{k + n}_{2n})$，而此时的 $f_1(\omega^k_n), f_2(\omega^k_n)$ 也可以通过类似的方法计算出来，直到
$$
f_1(\omega^k_1) = f_2(\omega^k_1) =1
$$
至此，我们就可以在 $O(log n)$ 的时间复杂度下计算出 $f(\omega^k_{2n})$ ，从而在 $O(n log n)$ 的时间复杂度下计算出我们要的所有 $f(\omega^i_n)$ 。

3 快速傅里叶变换逆变换

知道了如何进行快速傅里叶变换，接下来我们考虑如何实现它的逆变换。

为了方便表示，不妨设 $y_i = f(x_i)$
稍微整理一下： 我们现在要求解的是下面的方程组 $$ \left\{ \begin{aligned} a_0 + a_1 x_0 + a_2 x^2_0 + \dots + a_{n-1} x^{n-1}_0 & = f(x_0) \\
a_0 + a_1 x_1 + a_2 x^2_1 + \dots + a_{n-1} x^{n-1}_1 & = f(x_1) \\
a_0 + a_1 x_2 + a_2 x^2_2 + \dots + a_{n-1} x^{n-1}_2 & = f(x_2) \\
\vdots \\
a_0 + a_1 x_{n-1} + a_2 x^2_{n-1} + \dots + a_{n-1} x^{n-1}_{n-1} & = f(x_{n-1}) \end{aligned} \right. $$

要使 $y_i = f(x_i)$ 方便计算且满足两两不相等，我们可以令 $x_i = \omega^i_n$ ，即 $$ y_i = a_0 + a_1 \omega^i_n + a_2 (\omega^i_n)^2 + \dots + a_{n-1} (\omega^i_n)^{n-1} $$ 我们不妨把上面的方程组写成矩阵的形式： $$ \left( \begin{matrix} 1 & 1 & 1 & \dots & 1 \\ 1 & (\omega^1_n)^1 & (\omega^1_n)^2 & \dots & (\omega^1_n)^{n-1} \\ 1 & (\omega^2_n)^1 & (\omega^2_n)^2 & \dots & (\omega^2_n)^{n-1} \\ & \vdots & & \ddots & \vdots \\ 1 & (\omega^{n-1}_n)^1 & (\omega^{n-1}_n)^2 & \dots & (\omega^{n-1}n)^{n-1} \
\end{matrix}
\right)
\left(
\begin{matrix}
a_0 \\
a_1 \\
a_2 \\
\vdots \\
a_{n-1}
\end{matrix}
\right)
=
\left(
\begin{matrix}
y_0 \\
y_1 \\
y_2 \\
\vdots \\
y_{n-1} \end{matrix} \right) $$ 设左边的方阵就是对应的 $X$ ，显然 $X$ 可逆，接下来就是要考虑如何求 $X^{-1}$ ，知道了 $X^{-1}$ 我们就可以进行傅里叶逆变换。

这里再引入一个公式 $$ x^n - 1 = (x - 1) (1 + x + x^2 + \dots + x^{n-1}) $$ 当 $x^n - 1 = 0$ 时，显然 $x \neq 0$ ，所以 $1 + x + x^2 + \dots + x^{n-1} = 0$ ；
当 $x = 1$ 时，$1 + x + x^2 + \dots + x^{n-1} = n$ ，且满足等式。

再回到求 $X^{-1}$ 。
对于 $X$ 的第 $i$ 行 $$ \left( \begin{matrix} 1 & (\omega^i_n) & (\omega^i_n)^2 & \dots & (\omega^i_n)^{n-1} \end{matrix} \right) $$ 设 $X^{-1}$ 的第 $j$ 列 $$ \left( \begin{matrix} b_0 \\ b_1 \\ b_2 \\ \vdots \\ b_{n-1} \end{matrix} \right) $$ 当 $i \neq j$ 时，我们要想办法让下面的表达式等于 $0$ $$ b_0 + b_1 (\omega^i_n)^1 + b_2 (\omega^i_n)^2 + \dots + b_{n-1} (\omega^i_n)^{n-1} $$ 这时我们再回顾上面公式的第一种结论，有 $x = \omega^i_n$ ，等式依旧成立。所以我们就可以往这个方向去想，令 $X^{-1}$ 的第 $j$ 行为 $$ \left( \begin{matrix} 1 \\ (\omega^{n-j}_n) \\ (\omega^{n-j}_n)^2 \\
\vdots \\
(\omega^{n-j}_n)^{n-1}
\end{matrix}
\right)
$$
如此， $X^{-1}$ 的第 $v_{i j}$ 个元素可以表示为
$$
\begin{aligned}
v_{i j} & = 1 + \omega^{n-j+i}_n + (\omega^{n-j+i}_n)^2 + \dots + (\omega^{n-j+i}_n)^{n-1} \\ & = \left\{ \begin{aligned} 0 , i \neq j \\ n , i = j \end{aligned} \right. \end{aligned} $$ 利用这个性质我们可以通过构造将方阵 $XX^{-1}$ 的主对角线控制为 $n$ ，而其他元素控制为 $0$ ，例如下面的第二个方阵。 $$ \left( \begin{matrix} 1 & 1 & 1 & \dots & 1 \\ 1 & (\omega^1_n)^1 & (\omega^1_n)^2 & \dots & (\omega^1_n)^{n-1} \\ 1 & (\omega^2_n)^1 & (\omega^2_n)^2 & \dots & (\omega^2_n)^{n-1} \\ & \vdots & & \ddots & \vdots \\ 1 & (\omega^{n-1}_n)^1 & (\omega^{n-1}_n)^2 & \dots & (\omega^{n-1}_n)^{n-1} \end{matrix} \right) \left( \begin{matrix} 1 & 1 & 1 & \dots & 1 \\ 1 & (\omega^1_n)^{n-1} & (\omega^2_n)^{n-1} & \dots & (\omega^{n-1}_n)^{n-1} \\ 1 & (\omega^1_n)^{n-2} & (\omega^2_n)^{n-2} & \dots & (\omega^{n-1}_n)^{n-2} \\ & \vdots & & \ddots & \vdots \\ 1 & (\omega^1_n)^1 & (\omega^2_n)^1 & \dots & (\omega^{n-1}_n)^{n-1} \end{matrix} \right) = nE $$ 右边的矩阵就是我们要求的 $nX^{-1}$ ，再通过复数单位根的第三个结论可以得到 $$ \begin{aligned} X^{-1} & = \frac{1}{n} \left( \begin{matrix} 1 & 1 & 1 & \dots & 1 \\ 1 & (\omega^{n-1}_n)^1 & (\omega^{n-1}_n)^2 & \dots & (\omega^{n-1}_n)^{n-1} \\ 1 & (\omega^{n-2}_n)^1 & (\omega^{n-2}_n)^2 & \dots & (\omega^{n-2}_n)^{n-1} \\ & \vdots & & \ddots & \vdots \\ 1 & (\omega^1_n)^1 & (\omega^1_n)^2 & \dots & (\omega^{n-1}_n)^{n-1} \end{matrix} \right) \\ & = \frac{1}{n} \left( \begin{matrix} 1 & 1 & 1 & \dots & 1 \\ 1 & \bar{(\omega^1_n)}^1 & \bar{(\omega^1_n)}^2 & \dots & \bar{(\omega^1_n)}^{n-1} \\ 1 & \bar{(\omega^2_n)}^1 & \bar{(\omega^2_n)}^2 & \dots & \bar{(\omega^2_n)}^{n-1} \\ & \vdots & & \ddots & \vdots \\ 1 & \bar{(\omega^{n-1}_n)}^1 & \bar{(\omega^{n-1}_n)}^2 & \dots & \bar{(\omega^{n-1}_n)}^{n-1} \end{matrix} \right) \end{aligned} $$ 至此，我们就实现了快速傅里叶变换的逆变换。