文章目录

• 1. 多项式和内容概述
• 2. 多项式的表示
• • 2.1 系数表达（求值变为点值表达）
  • 2.2 点值表达（插值变为系数表达）
  • 2.3 基于系数表达的多项式的快速乘法
• 3. DFT与FFT
• 4. 高效FFT实现

两个 $n$ 次多项式相加的最直接方法所需的时间为 $\Theta (n)$ ，但是相乘的最直接方法所需的时间为 $\Theta (n^2)$ 。在这里我们将展示，快速傅里叶变换 the fast Fourier transform, or FFT 如何使多项式相乘的时间复杂度降低为 $\Theta (n\log n)$ 。

傅里叶变换的最常见用途是信号处理，这也是快速傅里叶变换的最常见用途，FFT也有很多日常应用，如压缩技术、编码数字视频和音频信息（包括MP3文件）。信号通常在时域中给出：作为一个把时间映射到振幅的函数 A signal is given in the time domain: as a function mapping time to amplitude（时域图显示了信号振幅随时间的变化情况）。傅里叶分析允许我们将「时域上的信号」表示成「不同频率的相移正弦波的加权叠加」express the signal as a weighted sum of phase-shifted sinusoids of varying frequencies 。和频率相关的权重和相位，在频域中刻画出信号的特征 The weights and phases associated with the frequencies characterize the signal in the frequency domain（频域图表示信号振幅与频率的关系，它只和峰值振幅与频率相关，不显示一个周期内的振幅变化，因为是正弦波）。

1. 多项式和内容概述

一个以 $x$ 为变量的多项式 polynomial 定义在一个代数域 algebraic field $F$ 上，将函数 $A(x)$ 表示为形式和：$$A(x)=\sum _{j=0}^{n-1}{a}_j{x}^j$$ 我们称 $a_0,a_1,\dots ,a_{n-1}$ 为如上多项式的系数 coefficients ，所有系数都属于域 $F$ ，典型的情形是复数集合 $\mathbb{C}$ 。如果一个多项式 $A(x)$ 的最高次的非零系数是 $a_k$ ，则称 $A(x)$ 的次数 degree 是 $k$ ，记为 $\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}(\mathrm{A})=\mathrm{k}$ 。任何严格大于一个多项式次数的整数，都是该多项式的次数界 degree-bound ，因此对次数界为 $n$ 的多项式，其次数可以是 $0\sim n-1$ 之间的任何整数，包括 $0$ 和 $n-1$ 。

我们在多项式上可以定义很多不同的运算。对于多项式加法 polynomial addition ，如果 $A(x)$ 和 $B(x)$ 是次数界为 $n$ 的多项式，则它们的和 $C(x)$ 也是一个次数界为 $n$ 的多项式，对所有属于定义域的 $x$ 都有 $C(x)=A(x)+B(x)$ such that C(x) = A(x) + B(x) for all x in the underlying field 。也就是说，如果 $$A(x)=\sum _{j=0}^{n-1}{a}_j{x}^j$$ 和 $$B(x)=\sum _{j=0}^{n-1}{b}_j{x}^j$$ 则 $$C(x)=\sum _{j=0}^{n-1}{c}_j{x}^j(c_j=a_j+b_j,j=0,1,\dots ,n-1)$$ 例如，如果有多项式 $A(x)=6x^3+7x^2-10x+9$ 和 $B(x)=-2x^3+4x-5$ ，那么 $C(x)=4x^3+7x^2-6x+4$ 。

对于多项式乘法 polynomial multiplication ，如果 $A(x)$ 和 $B(x)$ 都是次数界为 $n$ 的多项式，则它们的乘积 $C(x)$ 是一个次数界为 $2n-1$ 的多项式，对所有属于定义域的 $x$ 都有 $C(x)=A(x)B(x)$ such that C(x) = A(x)B(x) for all x in the underlying field 。多项式乘法的方法是：把 $A(x)$ 中的每一项与 $B(x)$ 中的每一项相乘，然后再合并同类项。例如，对两个多项式 $A(x)=6x^3+7x^2-10x+9$ 和 $B(x)=-2x^3+4x-5$ 进行如下的乘法：

另外一种表示乘积 $C(x)$ 的方法是：$$C(x)=\sum _{j=0}^{2n-2}{c}_j{x}^j\text{\hspace{1em}(30.1)}$$ 这里：$$c_j=\sum _{k=0}^j{a}_k{b}_{j-k}\text{\hspace{1em}(30.2)}$$ 注意，$\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}(\mathrm{C})=\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}(\mathrm{A})+\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}(\mathrm{B})$ ，意味着如果 $A$ 是次数界为 $n_a$ 的多项式，$B$ 是次数界为 $n_b$ 的多项式，那么 $C$ 是次数界为 $n_a+n_b-1$ 的多项式。因为一个次数界为 $k$ 的多项式也是次数界为 $k+1$ 的多项式，所以通常称乘积多项式 $C$ 是一个次数界为 $n_a+n_b$ 的多项式。

（算导30.1节）下面先介绍两种表示多项式的方法：系数表达和点值表达 the coefficient representation and the point-value representation 。当我们用系数表达表示多项式时，The straightforward methods for multiplying polynomials—equations (30.1) and (30.2) 所需运行时间为 $\Theta (n^2)$ ，但采用点值表达表示它们时，运行时间仅为 $\Theta (n)$ 。然而，我们通过在这两种表达法间进行转换、并使用系数表达求多项式乘积，可以使运行时间变为 $\Theta (n\log n)$ 。为了弄清这种方法为什么可行，需要（算导30.2节）首先学习单位复数根 complex roots of unity 。然后，我们运用FFT和它的逆变换 inverse（算导30.2节）来执行上述转换。之后（算导30.3节）展示如何在串行模型和并行模型上快速实现FFT。

2. 多项式的表示

从某种意义上，多项式的系数表达与点值表达是等价的，即用点值形式表示的多项式都对应唯一系数形式的多项式。下面介绍这两种表示方法，并展示如何结合这两种表达法、以使两个次数界为 $n$ 的多项式乘法运算在 $\Theta (n\log n)$ 时间内完成。

2.1 系数表达（求值变为点值表达）

对一个次数界为 $n$ 的多项式 $A(x)=\sum _{j=0}^{n-1}{a}_j{x}^j$ 而言，其系数表达是一个由系数组成的向量 $a=(a_0,a_1,\dots ,a_{n-1})$（矩阵方程中，一般将向量作为列向量看待）。

采用系数表达，对多项式的某些运算是很方便的。例如，对多项式 $A(x)$ 在给定点 $x_0$ 的求值 evaluation 运算（就是计算 $A(x_0)$ 的值），使用霍纳法则 Horner’s rule 我们可以在 $\Theta (n)$ 时间复杂度内完成求值运算：$$A(x_0)=a_0+x_0(a_1+x_0(a_2+\cdots +x_0(a_{n-2}+x_0(a_{n-1}))\dots ))$$ 类似地，对两个分别用系数向量 $a=(a_0,a_1,\dots ,a_{n-1})$ 和 $b=(b_0,b_1,\dots ,b_{n-1})$ 表示的多项式进行相加时，所需的时间是 $\Theta (n)$ ；我们仅输出系数向量 $c=(c_0,c_1,\dots ,c_{n-1})(c_j=a_j+b_j,j=0,1,\dots ,n-1)$ 。

现在考虑两个「用系数形式表示的、次数界为 $n$ 的多项式 $A(x),B(x)$ 」的乘法运算。如果用式30.1和式30.2中描述的方法，完成多项式乘法所需时间就是 $\Theta (n^2)$ ，因为向量 $a$ 中的每个系数必须与向量 $b$ 中的每个系数相乘。当用系数表示时，多项式乘法运算似乎要比求多项式值和多项式加法的运算更困难。由式30.2推导出的系数向量 $c$ ，也称为输入向量 $a$ 和 $b$ 的卷积 convolution ，表示成 $c=a\otimes b$ 。因为多项式乘法与卷积的计算都是最基本的计算问题，在实际应用中非常重要，所以重点讨论有关的高效算法。

2.2 点值表达（插值变为系数表达）

一个次数界为 $n$ 的多项式 $A(x)$ 的点值表达就是一个由 $n$ 个点值对 point-value pairs 所组成的集合： { ( x 0 , y 0 ) , ( x 1 , y 1 ) , … , ( x n − 1 , y n − 1 ) } \{ (x_0, y_0), (x_1, y_1), \dots, (x_{n-1}, y_{n-1} ) \} {(x0​,y0​),(x1​,y1​),…,(xn−1​,yn−1​)} 使得对 $k=0,1,\dots ,n-1$ ，所有 $x_k$ 各不相同，$$y_k=A(x_k)\text{\hspace{1em}(30.3)}$$

一个多项式可以有很多不同的点值表达，因为可以采用 $n$ 个不同的点 $x_0,x_1,\dots ,x_{n-1}$ 构成的集合作为这种表示方法的基。

对一个用系数形式表达的多项式来说，在原则上计算其点值表达是简单易行的，因为我们要做的就是选取 $n$ 个不同 $x_0,x_1,\dots ,x_{n-1}$ ，然后对 $k=0,1,\dots ,n-1$ 求出 $A(x_k)$ 。根据霍纳法则，求出这 $n$ 个点值所需时间复杂度为 $\Theta (n^2)$ 。稍后可以看到，如果巧妙地选取点 $x_k$ ，可以加速这一计算过程，使其运行时间变为 $\Theta (n\log n)$ 。

求值计算的逆（从一个多项式的点值表达确定其系数表达形式）称为插值 interpolation 。下面的定理说明，当插值多项式的次数界等于已知的点值对的数目，插值才是明确的（良定义的） interpolation is well defined when the desired interpolating polynomial must have a degree-bound equal to the given number of point-value pairs 。

定理30.1（插值多项式的唯一性 Uniqueness of an interpolating polynomial ）对于任意 $n$ 个点值对组成的集合 { ( x 0 , y 0 ) , ( x 1 , y 1 ) , … , ( x n − 1 , y n − 1 ) } \{ (x_0, y_0), (x_1, y_1), \dots, (x_{n-1}, y_{n-1})\} {(x0​,y0​),(x1​,y1​),…,(xn−1​,yn−1​)} ，其中所有的 $x_k$ 都不同，那么存在唯一的次数界为 $n$ 的多项式 $A(x)$ ，满足 $y_k=A(x_k)(k=0,1,\dots ,n-1)$ 。
证明：证明主要是根据某个矩阵存在逆矩阵。式30.3等价于矩阵方程：
$$\left[\begin{array}{ccccc}1& x_0& x_0^2& \dots & x_0^{n-1}\\ 1& x_1& x_1^2& \dots & x_1^{n-1}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& x_{n-1}& x_{n-1}^2& \dots & x_{n-1}^{n-1}\end{array}\right]\left[\begin{array}{c}{a}_0\\ a_1\\ ⋮\\ a_{n-1}\end{array}\right]=\left[\begin{array}{c}{y}_0\\ y_1\\ ⋮\\ y_{n-1}\end{array}\right]\text{\hspace{1em}(30.4)}$$ 左边的矩阵表示为 $V(x_0,x_1,\dots ,x_{n-1})$ ，称为范德蒙矩阵 Vandermonde matrix 。（根据算导思考题D-1）该 $n$ 阶范德蒙矩阵 $V_n$ 的行列式值 $D_n$ 为：$$\prod _{0\le j<k\le n-1}(x_k-x_j)$$ 这里的 $\prod _{0\le j<k\le n-1}(x_k-x_j)$ 表示所有同类因子 $(x_k-x_j)$（其中 $j<k$ ）的乘积，即：$$\begin{gathered} \prod _{0\le j<k\le n-1}(x_k-x_j)=(x_{n-1}-x_{n-2})(x_{n-1}-x_{n-3})\dots (x_{n-1}-x_0) \\ (x_{n-2}-x_{n-3})(x_{n-2}-x_{n-4})\dots (x_{n-2}-x_0)\dots (x_2-x_1)(x_2-x_0)(x_1-x_0) \end{gathered}$$

因此，根据（算导D.5）定理 $n\times n$ 矩阵 $A$ 是奇异的/不可逆的、当且仅当 $\mathrm{d}\mathrm{e}\mathrm{t}(A)=0$ 可知，如果 $x_k$ 都不同，则 $(x_k-x_j)$ 不会等于零，该矩阵的行列式非零，即该矩阵 $V$ 是可逆的 invertible（即非奇异的 nonsingular ）。因此，给定点值表达，我们能够唯一确定系数 $a_j$ ：$$a=V(x_0,x_1,\dots ,x_{n-1}{)}^{-1}y\text{\hspace{1em}■}$$

「定理30.1 插值多项式的唯一性」的证明过程，描述了基于求解线性方程组（式30.4）的一种插值算法。

• 利用（算导第28章）$LU$ 分解算法 the LU decomposition algorithms ，我们可以在 $O(n^3)$ 的时间复杂度内求出这些方程的解。
• 一种更快的、基于 $n$ 个点的插值算法，是基于如下的拉格朗日公式 Lagrange’s formula ：$$A(x)=\sum _{k=0}^{n-1}{y}_k\frac{{\prod }_{j\ne k}(x-x_j)}{{\prod }_{j\ne k}(x_k-x_j)}\text{\hspace{1em}(30.5)}$$ 可以自行验证，等式30.5的右端是一个次数界为 $n$ 的多项式，并满足对所有 $k$ 都有 $A(x_k)=y_k$ 。（算导练习30.1-5要求）运用拉格朗日公式，在 $\Theta (n^2)$ 的时间复杂度内计算 $A$ 的所有系数。

因此，（多项式的）$n$ 个点求值运算与插值运算，是良定义的互逆运算，它们将多项式的系数表达与点值表达相互转换 n-point evaluation and interpolation are well-defined inverse operations that transform between the coefficient representation of a polynomial and a point-value representation 。对于这些给出的问题，上面给出算法的运行时间为 $\Theta (n^2)$ 。

对于很多多项式相关的操作，点值表达是很便利的。对于加法，如果 $C(x)=A(x)+B(x)$ ，则对任意点 $x_k$ 满足 $C(x_k)=A(x_k)+B(x_k)$ 。更准确的说，如果已知 $A$ 的点值表达 { ( x 0 , y 0 ) , ( x 1 , y 1 ) , … , ( x n − 1 , y n − 1 ) } \{ (x_0, y_0), (x_1, y_1), \dots, (x_{n-1}, y_{n-1}) \} {(x0​,y0​),(x1​,y1​),…,(xn−1​,yn−1​)} 和 $B$ 的点值表达 （注意，$A$ 和 $B$ 在相同的 $n$ 个位置求值） { ( x 0 , y 0 ′ ) , ( x 1 , y 1 ′ ) , … , ( x n − 1 , y n − 1 ′ ) } \{ (x_0, y_0'), (x_1, y_1'), \dots, (x_{n-1}, y_{n-1}') \} {(x0​,y0′​),(x1​,y1′​),…,(xn−1​,yn−1′​)} 则 $C$ 的点值表达是： { ( x 0 , y 0 + y 0 ′ ) , ( x 1 , y 1 + y 1 ′ ) , … , ( x n − 1 , y n − 1 + y n − 1 ′ ) } \{ (x_0, y_0 + y_0'), (x_1, y_1 + y_1'), \dots, (x_{n-1}, y_{n-1}+y_{n-1}') \} {(x0​,y0​+y0′​),(x1​,y1​+y1′​),…,(xn−1​,yn−1​+yn−1′​)} 因此，对两个点值形式表达的次数界为 $n$ 的多项式相加，所需时间复杂度为 $\Theta (n)$ 。

类似地，对于多项式乘法，点值表达也是方便的。如果 $C(x)=A(x)B(x)$ ，则对于任意点 $x_k$ 满足 $C(x_k)=A(x_k)B(x_k)$ ，并且对 $A$ 的点值表达和 $B$ 的点值表达进行逐点相乘，就可得到 $C$ 的点值表达。不过，我们也必须面对这样一个问题，即 $\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}(\mathrm{C})=\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}(\mathrm{A})+\mathrm{d}\mathrm{e}\mathrm{g}\mathrm{r}\mathrm{e}\mathrm{e}(\mathrm{B})$（如果 $A,B$ 的次数界为 $n$ ，那么 $C$ 的次数界为 $2n$ ） 。也就是说，对于 $A$ 和 $B$ 每个多项式而言，一个标准点值表达是由 $n$ 个点值对所组成。当我们把这些点值对相乘，就得到 $C$ 的 $n$ 个点值对，然而由于 $C$ 的次数界为 $2n$ ，要插值获得唯一的多项式，我们就需要 $2n$ 个点值对（算导练习30.1-4）。因此，必须对 $A$ 和 $B$ 的点值表达进行“扩展”，使每个多项式都包含 $2n$ 个点值对。给定 $A$ 的扩展点值表达： { ( x 0 , y 0 ) , ( x 1 , y 1 ) , … , ( x 2 n − 1 , y 2 n − 1 ) } \{ (x_0, y_0), (x_1, y_1), \dots, (x_{2n-1}, y_{2n-1}) \} {(x0​,y0​),(x1​,y1​),…,(x2n−1​,y2n−1​)} 和 $B$ 的对应扩展点值表达： { ( x 0 , y 0 ′ ) , ( x 1 , y 1 ′ ) , … , ( x 2 n − 1 , y 2 n − 1 ′ ) } \{ (x_0, y_0'), (x_1, y_1'), \dots, (x_{2n-1}, y_{2n-1}')\} {(x0​,y0′​),(x1​,y1′​),…,(x2n−1​,y2n−1′​)} 则 $C$ 的点值表达为： { ( x 0 , y 0 y 0 ′ ) , ( x 1 , y 1 y 1 ′ ) , … , ( x 2 n − 1 , y 2 n − 1 y 2 n − 1 ′ ) } \{ (x_0, y_0y_0'), (x_1, y_1y_1'), \dots, (x_{2n-1}, y_{2n-1}y_{2n-1}') \} {(x0​,y0​y0′​),(x1​,y1​y1′​),…,(x2n−1​,y2n−1​y2n−1′​)}

给定两个点值扩展形式的输入多项式，我们可以看到「使其相乘而得到点值形式的结果」需要 $\Theta (n)$ 时间，比「采用系数形式表示的多项式相乘所需时间」少得多。

最后我们考虑，对一个采用点值表达的多项式、如何求其在某个新点上的值。对这个问题，最简单不过的方法就是，先把该多项式转换为系数形式表达，然后在新点处求值。

2.3 基于系数表达的多项式的快速乘法

我们能否利用「基于点值表达的多项式」的线性时间乘法算法，来加速「基于系数表达的多项式」乘法运算呢？答案的关键，在于能否快速把一个多项式从系数形式转换为点值形式（求值）、以及从点值形式转换为系数形式（插值）。

我们可以采用任何点作为求值点，但通过精心地挑选求值点，可以把两种表示法之间转化所需的时间复杂度变为 $\Theta (n\log n)$ 。（算导30.2节说明）如果选择单位复数根作为求值点，我们可以通过对系数向量进行离散傅里叶变换 discrete Fourier transform, DFT ，得到相应的点值表达；我们也可以通过对点值对执行逆DFT变换 inverse DFT ，而获得相应的系数向量，这样就实现了求值运算的逆运算（插值）。（算导30.2节说明了）FFT如何在 $\Theta (n\log n)$ 的时间复杂度内完成DFT和逆DFT运算。

图30-1说明了这种策略。其中一个细节涉及次数界：两个次数界为 $n$ 的多项式乘积是一个次数界为 $2n$ 的多项式。因此在对输入多项式 $A$ 和 $B$ 进行求值以前，首先给这两个多项式添加 $n$ 个为 $0$ 的高次系数、使其次数界加倍为 $2n$ we first double their degree-bounds to 2n by adding n high-order coefficients of 0 。因为现在这些向量有 $2n$ 个元素，我们可以采用 $2n$ 次单位复数根 complex (2n)th roots of unity ，它在图30-1中被标记为 ${\omega }_{2n}$ 。

给出FFT，我们就有下面的时间复杂度为 $\Theta (n\log n)$ 的方法，该方法把两个次数界为 $n$ 的多项式 $A(x),B(x)$ 进行乘法运算，其中输入与输出均采用系数表达。我们假定 $n$ 是 $2$ 的幂；我们总能通过添加系数为 $0$ 的高阶系数，来满足这个要求 We assume that n is a power of 2; we can always meet this requirement by adding high-order zero coefficients 。

1. 加倍次数界 Double degree-bound ：通过加入 $n$ 个系数为 $0$ 的高阶系数，构造「 $A(x)$ 和 $B(x)$ 为次数界 $2n$ 的多项式」的系数表达 Create coefficient representations of A(x) and B(x) as degree-bound 2n polynomials 。
2. 求值 Evaluate ：通过应用 $2n$ 阶的FFT，计算出 $A(x)$ 和 $B(x)$ 的、长度为 $2n$ 的点值表达。这些点值表达中，包含了两个多项式在 $2n$ 次单位根处的取值。
3. 逐点相乘 Pointwise multiply ：把 $A(x)$ 的值与 $B(x)$ 的值逐点相乘，可以计算出多项式 $C(x)=A(x)B(x)$ 的点值表达，这个表示中包含了 $C(x)$ 在每个 $2n$ 次单位根处的值。
4. 插值 Interpolate ：通过对 $2n$ 个点值对应用FFT来计算其逆FFT applying the FFT on 2n point-value pairs to compute the inverse DFT ，就可以构造出多项式 $C(x)$ 的系数表达。

执行第1步和第3步所需时间为 $\Theta (n)$ ，执行第2步和第4步所需时间为 $\Theta (n\log n)$ 。因此，一旦表明如何运用FFT，我们就已经证明了下面的定理。

定理30.2 当输入和输出多项式均采用系数表达时，我们就能在 $\Theta (n\log n)$ 时间复杂度内，计算出两个次数界为 $n$ 的多项式的乘积。$■$

3. DFT与FFT

4. 高效FFT实现

【算法学习笔记】快速傅里叶变换2