Lecture 5: FFT，高斯消元与LU分解

多项式乘法

考虑一个多项式乘法任务：给定两个系数最多为 $n$ 次的多项式 $p(x), q(x)$，求 $r(x) = p(x) \cdot q(x)$

通常，我们已知的是两个多项式的系数，对系数进行乘积计算，朴素的时间复杂度为 $O(n^{2})$

如果，我们得到的是两个多项式在 $n+1$ 个 $x$ 取值对应的点对 $(x_{i}, p(x_i), q(x_i))$，那么利用多项式插值的唯一性，我们可以计算出最多 $2n+1$ 个 $r(x)$ 上的点，从而唯一确定 $r(x)$，时间复杂度为 $O(n)$

可见，多项式的系数表示法和点值表示法对于多项式乘法的时间复杂度不同。我们希望借助点值表示法加速系数表示法下的多项式乘法，但朴素情况下，系数表示与点值表示之间的转换非常耗时（点值表示 → 系数表示需要解线性方程组，时间复杂度达到 $O(n^{3})$）

注意到点值选择是任意的 $n+1$ 个不同点，我们是否可以精心选择一些特殊的点，使得计算过程可以被压缩？

我们的想法是，已知 $p(x),q(x)$ 的系数不超过 $n$

令 $m \geq 2n+1$，选定 $x_{0}, \cdots, x_{m}$，计算 $p(x_{0}),\cdots, p(x_{m})$ 和 $q(x_{0}), \cdots, q(x_{m})$
- 时间复杂度 $O((m+1)n) = O(n^2)$
计算这 $m$ 对 $r(x_{j}) = p(x_{j}) \cdot q(x_{j})$
- 时间复杂度 $O(m) = O(n)$
对 $\lbrace r(x_{j}) \rbrace$ 插值，得到 $r(x) = p(x)q(x)$ 的系数
- 拉格朗日插值时间复杂度 $O(m^{3}) = O(n^{3})$，牛顿法为 $O(m^{2}) = O(n^{2})$

看来时间瓶颈出现在第一步和第三步，有没有什么方法能够将这两步的时间复杂度下降到 $O(n\log n)$ 呢？

复数与单位根

一个复数在代数形式和极坐标形式下分别表示为 $z = a+bi = re^{i\theta}$

指定 $n$，我们定义关于 $z$ 的方程 $z^{n} = 1$ 的解集合为 $n$ 次单位根，比如：

2 次单位根为 $±1$
4 次单位根为 $±1,±i$

$n$ 个单位根在复平面单位圆上呈等距分布，更具体的，可以通过欧拉公式求解： $$ e^{\frac{2\pi l}{n}i} = \cos\frac{2\pi l}{n} + i\sin\frac{2\pi l}{n} $$ 我们需要关注 $n=2^{k}$ 的情况，此时可以得出一个性质：

从上往下为 $2^{k}$ 次单位根的树形表示，注意到 $k>0$ 时：

如果 $x$ 是一个 $2^{k}$ 次单位根，则 $-x$ 也是一个 $2^{k}$ 次单位根（并且 $x, -x$ 是 $2^{k+1}$ 次单位根）

傅里叶变换

离散傅里叶变换

先从典型的离散傅里叶变换开始：

Definition: 离散傅里叶变换（Discrete Fourier Transform, DFT）

给定输入：多项式 $p(x)$ 的系数 $a_{0},\cdots, a_{n}$

输出：选定 $x_{0}, \cdots, x_{m}$ 为 $m+1$ 次单位根，并计算 $p(x_{0}), p(x_{m})$ 的值作为输出：$\lbrace(x_{i}, p(x_{i}))\rbrace$

同时也有逆变换，交换输入和输出，利用插值法求算

（这里给出的定义没有进行归一化，否则会将计算结果 $p(x_{i})$ 除以 $\sqrt{n+1}$ 进行归一化）

我们用一个式子表示 DFT 的计算： $$ p(w^{l}) = \sum_{j=0}^{n} a_{j} \omega^{lj} $$ 以及逆运算： $$ a_{l} = \frac{1}{n+1} \sum_{j=0}^{n} p(\omega^{j}) \omega^{-lj} $$ 注意到，离散傅里叶变换及其逆变换，对应了多项式求值时的两步计算瓶颈过程。如果我们能结合 $2^{k}$ 次单位根的树形性质，写出更快的傅里叶变换，就能降低多项式求值的计算量

快速傅里叶变换

首先，我们改写 $p(x)$。令：

$E(x) = a_{0} + a_{2}x + a_{4}x^{2} + \cdots + a_{n-2}x^{\frac{n}{2}-1}$
$O(x) = a_{1} + a_{3}x + a_{5}x^{2} + \cdots + a_{n-1}x^{\frac{n}{2}-1}$
- 它们的系数都是确定的

得到 $p(x) = E(x^{2}) + xO(x^{2})$，同理，$p(-x) = E(x^{2}) - xO(x^{2})$

可见，只要计算 $E(x^{2}),O(x^{2})$，就可以通过一次加法、减法分别计算出 $p(x), p(-x)$ 的值。我们已经发现，如果 $x$ 是 $2^{k}$ 次单位根，则 $-x$ 也是 $2^{k}$ 次单位根，因此在这里对 $p(-x)$ 的计算也是必要的

现在我们将对 $p(x)$ 与 $p(-x)$ 的计算转换成了两个新的问题：如何计算 $E(x^{2}), O(x^{2})$？比较显然的，如果 $x$ 是 $2^{k}$ 次单位根，则 $x^{2}$ 是 $2^{k-1}$ 次单位根。因此，我们可以继续分解问题，这个计算过程最终成为了一个递归的树结构（或者说，分治操作）

我们直接给出 FFT 的伪代码：

快速傅里叶变换

\[ \begin{array}{l} \textbf{The fast Fourier transform (polynomial formulation)} \newline \hline \texttt{function FFT}(A, \omega) \newline \textbf{Input: } \text{Coefficient representation of a polynomial } A(x) \text{ of degree } \le n-1, \newline \phantom{\textbf{Input:}} \text{where } n \text{ is a power of 2; } \omega \text{, an } n\text{th root of unity} \newline \textbf{Output: } \text{Value representation } A(\omega^0), \dots, A(\omega^{n-1}) \newline \newline \textbf{if } \omega = 1 \textbf{: return } A(1) \newline \text{express } A(x) \text{ in the form } A_e(x^2) + x A_o(x^2) \newline \textbf{call } \text{FFT}(A_e, \omega^2) \text{ to evaluate } A_e \text{ at even powers of } \omega \newline \textbf{call } \text{FFT}(A_o, \omega^2) \text{ to evaluate } A_o \text{ at even powers of } \omega \newline \textbf{for } j = 0 \textbf{ to } n-1 \textbf{:} \newline \quad \text{compute } A(\omega^j) = A_e(\omega^{2j}) + \omega^j A_o(\omega^{2j}) \newline \newline \textbf{return } A(\omega^0), \dots, A(\omega^{n-1}) \end{array} \]

简单来说，我们递归分解多项式计算出各个 $E(x^k), O(x^k)$ 的值，然后借助 $p(x) = E(x^{2}) + xO(x^{2})$ 合并计算出 $p(x_{i})$ 的值。

能够证明，记 $T(n)$ 为 $n$ 阶 FFT 需要的操作次数，则 $T(n) = 2T\left( \dfrac{n}{2} \right) + C_n$（其中 $C_n$ 是与阶数相关的常数），最终可以计算证明：$T(n) = O (n\log n)$

我们记 DFT 的计算结果 $p(w^{l}) = \sum_{j=0}^{n} a_{j} \omega^{lj}$，逆运算 $a_{l} = \frac{1}{n+1} \sum_{j=0}^{n} p(\omega^{j}) \omega^{-lj}$，区别在于把 $\omega \to \omega^{-1}$，且系数多了 $\frac{1}{n+1}$。这对 FFT 也是正确的。

逆 FFT 也只要在 FFT 的基础上，修改上述两处内容即可

矩阵表达形式

挖坑

2026-04-132026-04-13