快速傅里叶变换及其实现

标签：right 变换 text omega 傅里叶 多项式 快速 left

第1章 引言

傅里叶变换（Fourier Transform）是由数学家傅里叶提出的一套对函数进行变换的方法，其主要分为连续傅里叶变换（Continuous Fourier Transform，CFT）和离散傅里叶变换（Discrete Fourier Transform，DFT）两种，在本文中，我们只研究离散傅里叶变换。

离散傅里叶变换虽然在数学层面很有用，但其算法的时间复杂度较高，在算法层面并不实用。继而，后续研究者又提出了快速傅里叶变换（Fast Fourier Transform，FFT）算法，这才彻底解决了问题。

那么，离散傅里叶变换到底有什么用呢？它的用途十分直白：用于计算多项式乘法。

多项式乘法早在中学数学中就已经学过，例如：

\[\left( 1+2x+3x^2 \right) \left( 4+5x+6x^2 \right) =4+13x+28x^2+27x^3+18x^4 \]

相信读者一定有这样的感觉：多项式的加减法很好算，只需要做几次合并同类项就行了，这是因为多项式的加减法是时间复杂度为\(\varTheta(N)\)的算法。但多项式的乘法可就很难算了，上式中，想要计算一次3*3的多项式乘法，则需要进行9次乘法计算，和少量的合并同类项计算。也就是说，如果想要计算一次N*N的多项式乘法，则需要进行\(N^2\)次乘法计算，和少量的合并同类项计算，即：这种计算多项式乘法的算法的时间复杂度为\(\varTheta(N^2)\)。

快速傅里叶变换正是这样一个算法：其能够突破上述算法的\(\varTheta(N^2)\)时间复杂度，将多项式乘法的时间复杂度优化至\(\varTheta(NlogN)\)。

第2章 多项式的系数表达与点值表达

想要研究多项式，就需要先把多项式写出来。在本章中，我们研究多项式的两种表达方式：系数表达与点值表达。

2.1 系数表达

多项式的系数表达我们都非常熟悉：指的就是上面的\(\left( 1+2x+3x^2 \right)\)和\(\left( 4+5x+6x^2 \right)\)这种形式，通过写出多项式中每一项的系数，从而表达出一个多项式是什么样的。

这种写法还可以再省略一些：由于每个系数后面的\(x^n\)写不写出来都一样，所以可以只写出每一项的系数，并构成一个向量：

\[1+2x+3x^2\,\,\Rightarrow \,\,\left( 1, 2, 3 \right) \\ 4+5x+6x^2\,\,\Rightarrow \,\,\left( 4, 5, 6 \right) \]

此时，多项式乘法就有了一种新的表示：

\[\left( 1, 2, 3 \right) \otimes \left( 4, 5, 6 \right) =\left( 4, 13, 28, 27, 18 \right) \]

这是一种全新的向量间的乘法运算，称为：卷积（Convolution），用符号\(\otimes\)表示。

2.2 点值表达

多项式的系数表达很好理解，那么，点值表达又是什么呢？

我们都知道：两点确定一条直线，而直线是一个包含常数项和一次项的多项式；所以，我们也可以说：两点确定一个一次多项式。那么，三点能不能确定一个二次多项式呢？四点又能不能确定一个三次多项式呢？更多的点呢？

答案是肯定的，请看如下引理：

引理1（点值表达的唯一性）：对于任意的n个点构成的点集：\(\left\{ \left( x_0, y_0 \right) , \left( x_1, y_1 \right) , \left( x_2, y_2 \right) , ..., \left( x_{n-1}, y_{n-1} \right) \right\}\)，如果\(x_0\ne x_1\ne x_2\ne ...\ne x_{n-1}\)，则该点集能够唯一确定一个\(n-1\)次多项式。

证明：

\[\text{设}n-1\text{次多项式：}A\left( x \right) =a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1} \\ \text{将}\left( x_0, y_0 \right) \text{，}\left( x_1, y_1 \right) \text{，}\left( x_2, y_2 \right) \text{，}...\text{，}\left( x_{n-1}, y_{n-1} \right) \text{依次带入，得方程组：} \\ \begin{cases} a_0+a_1x_0+a_2{x_0}^2+...+a_{n-1}{x_0}^{n-1}=y_0\\ a_0+a_1x_1+a_2{x_1}^2+...+a_{n-1}{x_1}^{n-1}=y_1\\ a_0+a_1x_2+a_2{x_2}^2+...+a_{n-1}{x_2}^{n-1}=y_2\\ ...\\ a_0+a_1x_{n-1}+a_2{x_{n-1}}^2+...+a_{n-1}{x_{n-1}}^{n-1}=y_{n-1}\\ \end{cases} \\ \text{将上式写为矩阵方程形式：} \\ \left[ \begin{matrix} 1& x_0& x_{0}^{2}& \cdots& x_{0}^{n-1}\\ 1& x_1& x_{1}^{2}& \cdots& x_{1}^{n-1}\\ 1& x_2& x_{2}^{2}& \cdots& x_{2}^{n-1}\\ \vdots& \vdots& \vdots& \ddots& \vdots\\ 1& x_{n-1}& x_{n-1}^{2}& \cdots& x_{n-1}^{n-1}\\ \end{matrix} \right] \left[ \begin{array}{c} a_0\\ a_1\\ a_2\\ \cdots\\ a_{n-1}\\ \end{array} \right] =\left[ \begin{array}{c} y_0\\ y_1\\ y_2\\ \cdots\\ y_{n-1}\\ \end{array} \right] \\ \text{最左边的矩阵为范德蒙德矩阵，且由于任意的}x_i\ne x_j\text{，所以其行列式：} \\ \prod_{0\leqslant i<j\leqslant n-1}{\left( x_j-x_i \right)}\ne 0 \\ \text{所以，矩阵方程存在唯一解，得证。} \]

也就是说，我们不仅可以用n个系数唯一的表达一个n-1次多项式，还可以用这个多项式上的n个点唯一表达。

现在，需要考虑三个问题：

1. 已知一个多项式的系数表达，怎么得到其点值表达？
2. 已知一个多项式的点值表达，怎么得到其系数表达？
3. 点值表达存在的意义是什么？

第一个问题很简单，我们已经知道，想要得到一个多项式的点值表达，需要满足两个条件：

1. 点的数量要够，需要找n个点
2. 每个点的x互不相同

这两个条件可太好满足了，专门挑几个最好算的x带入计算，点不就来了吗？例如：对于多项式\(1+2x+3x^2\)，需要找3个点，那我们就选3个最好算的点：0、1、-1，带进去，点就有了：\(\left\{ \left( 0, 1 \right) , \left( 1, 6 \right) , \left( -1, 2 \right) \right\}\)，这就是多项式\(1+2x+3x^2\)的点值表达。当然，还可以挑别的3个点，算出来的点集都能唯一确定\(1+2x+3x^2\)这个多项式。

第二个问题也不难解决，不是有很多点吗，那就用这些点解上面证明过程中的那个线性方程组，解出来的\(a_0, a_1, a_2, ..., a_{n-1}\)就是系数表达了。

第三个问题很重要。我们之所以要研究点值表达，显然是为了解决多项式的乘法问题的。说到这，点值表达存在的意义就出现了：点值表达下的多项式四则运算，全都是时间复杂度为\(\varTheta(N)\)的算法。这句话可能需要想一下才能明白，这是因为：两个多项式做四则运算，其实就是在这两个多项式的每一对y上都做一遍四则运算。那么，如果我们选择相同的一组x带入两个多项式，得到两组点值表达，多项式乘法就可以变成这两组点之间做一遍乘法了。举个例子：比如我们对\(\left( 1+2x+3x^2 \right)\)和\(\left( 4+5x+6x^2 \right)\)这两个多项式都选0、1、-1带入，就能得到两组x相同的点值表达：\(\left\{ \left( 0, 1 \right) , \left( 1, 6 \right) , \left( -1, 2 \right) \right\}\)和\(\left\{ \left( 0, 4 \right) , \left( 1, 15 \right) , \left( -1, 5 \right) \right\}\)；此时，多项式乘法就可以用这三个点的y进行计算，结果是：\(\left\{ \left( 0, 4 \right) , \left\{ 1, 90 \right\} , \left( -1, 10 \right) \right\}\)。

事实上，我们在这里犯了一个错误：两个二次多项式做乘法，结果应当是一个四次多项式，而四次多项式需要用五个点才能唯一表示，而上面只有三个点，显然是不够的。这个问题提醒了我们：在计算多项式乘法的时候，不能只看现在这个多项式是几次的，还应当看乘法的结果多项式是几次的，后者的次数才能决定一开始要取几个点。

三个问题都讨论完了，但是除了第三个问题的结论很有吸引力以外，前两个问题的结论实在是平平无奇。

先看从系数表达转点值表达的过程，想做这件事，就要选n个点，每个点依次带入多项式算一遍，而多项式里面全是各种高次幂，时间复杂度早已不可接受了。如果使用霍纳法则，将计算一个多项式的时间复杂度优化到\(\varTheta(N)\)，那最终也是一个时间复杂度为\(\varTheta(N^2)\)的算法。

而从点值表达转系数表达就更难算了，由于需要解矩阵方程，考虑线性代数中的LU分解算法，其时间复杂度是\(O(N^3)\)。拉格朗日插值法则更快一些，其基于以下公式：

\[A\left( x \right) =\sum_{i=0}^{n-1}{\left( y_i\frac{\prod_{j\ne i}{\left( x-x_j \right)}}{\prod_{j\ne i}{\left( x_i-x_j \right)}} \right)} \]

不难看出，这个公式的时间复杂度为\(\varTheta(N^2)\)。

看来，我们遇到了木桶效应：虽然在点值表达下，多项式的乘法变得非常好算，但两种表达方式间的来回转换，全都是时间复杂度至少也要\(\varTheta(N^2)\)的算法。还不如用最普通的那个算法呢。

所以，现在的目标是：找到一对算法，能够更快的在多项式的系数表达和点值表达之间进行转换。

第3章 欧拉公式

由于傅里叶变换需要使用复数的相关知识，所以这一章中，我们研究欧拉公式及其推论。

欧拉公式是一个非常有名的公式，其将复数域上的指数函数和三角函数联系在了一起。让我们从\(e^{i\theta}\)的麦克劳林级数展开开始：

\[e^{i\theta}=\sum_{n=0}^{\infty}{\frac{\left( i\theta \right) ^n}{n!}}=1+i\theta -\frac{\theta ^2}{2!}-\frac{i\theta ^3}{3!}+\frac{\theta ^4}{4!}+\frac{i\theta ^5}{5!}-... \\ \text{将上式按奇偶顺序重新排列并整理得到：} \\ =\left( 1-\frac{\theta ^2}{2!}+\frac{\theta ^4}{4!}-... \right) +i\left( \theta -\frac{\theta ^3}{3!}+\frac{\theta ^5}{5!}-... \right) \\ =\sum_{n=0}^{\infty}{\frac{\left( -1 \right) ^n}{\left( 2n \right) !}\theta ^{2n}+i\sum_{n=0}^{\infty}{\frac{\left( -1 \right) ^n}{\left( 2n+1 \right) !}\theta ^{2n+1}}} \\ \text{上式的左右两边分别是}cos\theta \text{和}sin\theta \text{的麦克劳林级数：} \\ =\cos \theta +i\sin \theta \]

这就是著名的欧拉公式：

引理2（欧拉公式）：

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

将\(\theta=\pi\)带入欧拉公式，就能得到著名的欧拉恒等式：

引理3（欧拉恒等式）：

\[e^{\pi i}=-1 \]

欧拉公式和欧拉恒等式在后续的研究中十分重要，读者应熟练掌握。

第4章 n次单位复数根

这一章中，我们研究n次单位复数根。什么是n次单位复数根呢？其指的是以下n次方程的根：

\[\omega ^n=1 \]

在实数域中，不管n是多少，这个方程都最多只有两个根：1和-1。但是，由代数学基本定理：任何一个n次方程都有且仅有n个根。那么，对于\(n>2\)的这些高次方程来说，剩下的根去哪了呢？显然，这些根都是虚数。

根据欧拉公式，我们可以凑出n次单位复数根的一般形式：

引理4（n次单位复数根）：

\[\text{方程：}\omega ^n=1\text{的所有根依次为：}\left( e^{\frac{2\pi i}{n}} \right) ^k, k=0, 1, 2, ..., n-1\text{；记作：}\omega _{n}^{k} \]

证明：

\[\text{将}\left( e^{\frac{2\pi i}{n}} \right) ^k\text{带入原方程得：} \\ \left( e^{\frac{2\pi i}{n}} \right) ^{kn}=e^{2\pi ik} \\ \text{由欧拉公式：} \\ e^{2\pi ik}=\cos 2\pi k+i\sin 2\pi k=1 \\ 得证。 \]

由复数在复平面上的极坐标表示可知：各个n次单位复数根是复平面上单位圆的各个n等分点，故其具有一些特殊的性质。请看：

引理5（折半引理）：\(\omega _{n}^{k}=\omega _{n/2}^{k/2}\)

证明：

\[\omega _{n}^{k}=e^{\frac{2\pi ik}{n}}=e^{\frac{2\pi i\frac{k}{2}}{\frac{n}{2}}}=\omega _{n/2}^{k/2} \]

引理6：\(\omega _{n}^{n}=1\)

证明：

\[\omega _{n}^{n}=e^{\frac{2\pi in}{n}}=e^{2\pi i}=\left( -1 \right) ^2=1 \]

引理7：\(\omega _{n}^{n/2}=-1\)

证明：

\[\omega _{n}^{n/2}=e^{\frac{2\pi i\frac{n}{2}}{n}}=e^{\pi i}=-1 \]

引理8：\(\omega _{n}^{-k}=\overline{\omega _{n}^{k}}\)

证明：

\[\because \omega _{n}^{k}=e^{\frac{2\pi ik}{n}}=\cos \left( \frac{2\pi k}{n} \right) +i\sin \left( \frac{2\pi k}{n} \right) \\ \text{又}\because \omega _{n}^{-k}=e^{\frac{-2\pi ik}{n}}=\cos \left( -\frac{2\pi k}{n} \right) +i\sin \left( -\frac{2\pi k}{n} \right) =\cos \left( \frac{2\pi k}{n} \right) -i\sin \left( \frac{2\pi k}{n} \right) \\ \therefore \omega _{n}^{-k}=\overline{\omega _{n}^{k}} \]

n次单位复数根及其若干引理的用途将在后续章节展开。

第5章 离散傅里叶变换与离散傅里叶逆变换

上一章中，我们研究了n次单位复数根及其性质。那么，这些根有什么用呢？傅里叶变换现在正式登场：

使用所有的n次单位复数根将一个n-1次多项式从系数表达转为点值表达的过程，就称为离散傅里叶变换（Discrete Fourier Transform，DFT）。

原来，上一章中研究的这些n次单位复数根是用来带入的。如果选用这一组根作为x带入多项式，从而得到其点值表达，这个过程就称为离散傅里叶变换。

那么，离散傅里叶逆变换又是什么呢？顾名思义，其是一个变换回来的过程：

将离散傅里叶变换得到的点值表达转为系数表达的过程，就称为离散傅里叶逆变换（Inverse Discrete Fourier Transform，IDFT）。

我们已经知道，想要将一个多项式从系数表达转为点值表达，随便选一组点都是可以的。那为什么选n次单位复数根这组点进行转换，就有了傅里叶变换和傅里叶逆变换这两个专业术语呢？我们不禁猜想：n次单位复数根可能是一组非常特殊的点，将这组点带入多项式时，是可以简化计算的；此外，从这些点得到的点值表达向系数表达转换时，也是可以简化计算的。这样一来，双向的转换过程就都能得到优化了。

事实上，这两个猜想都是正确的。其算法分别被称为：快速傅里叶变换（Fast Fourier Transform，FFT）和快速傅里叶逆变换（Inverse Fast Fourier Transform，IFFT）。在下面的两章中，我们分别研究这两种算法。

第6章 快速傅里叶变换

6.1 快速傅里叶变换的数学原理

上一章中，我们已经知道，如果将n次单位复数根带入多项式，进行系数表达到点值表达的转换，是可以简化计算的，这样的算法被称为快速傅里叶变换。这一章中，我们具体研究这一算法。

既然是一种优化算法，那就不能一上来就带入点进行计算，我们需要先做一些准备：

1. 因为快速傅里叶变换是一个严格二分的算法，所以需要将多项式的项数补齐至2的整数次幂。什么意思呢？比如想对一个二次多项式进行快速傅里叶变换，由于二次多项式有只有三项，所以需要用系数0再补一项，将其补至四项；例如：\(1+2x+3x^2\)就需要补成\(1+2x+3x^2+0x^3\)。同理，具有5、6、7项的多项式都需要用系数0补齐至8项；以此类推
2. 将多项式按奇偶顺序重排并整理：

\[\text{设}n-1\text{次多项式：}A\left( x \right) =a_0+a_1x+a_2x^2+a_3x^3...a_{n-1}x^{n-1} \\ \text{按奇偶顺序重排此多项式，得：} \\ A\left( x \right) =\left( a_0+a_2x^2+...+a_{n-2}x^{n-2} \right) +x\left( a_1+a_3x^2+...+a_{n-1}x^{n-2} \right) \]

3. 将多项式分解为两个子多项式：

\[\text{对于：}A\left( x \right) =\left( a_0+a_2x^2+...+a_{n-2}x^{n-2} \right) +x\left( a_1+a_3x^2+...+a_{n-1}x^{n-2} \right) \\ \text{设：}A^{\left[ 0 \right]}\left( x \right) =a_0+a_2x+a_4x^2+...+a_{n-2}x^{\frac{n}{2}-1} \\ \text{且：}A^{\left[ 1 \right]}\left( x \right) =a_1+a_3x+a_5x^2+...+a_{n-1}x^{\frac{n}{2}-1} \\ \text{则：}A\left( x \right) =A^{\left[ 0 \right]}\left( x^2 \right) +xA^{\left[ 1 \right]}\left( x^2 \right) \]

至此，准备工作就完成了。

现在，我们将n次单位复数根，即\(\omega _{n}^{k}\)带入上式：

\[\text{对于：}A\left( x \right) =A^{\left[ 0 \right]}\left( x^2 \right) +xA^{\left[ 1 \right]}\left( x^2 \right) \\ \text{将}x=\omega _{n}^{k}\text{带入得：} \\ A\left( \omega _{n}^{k} \right) =A^{\left[ 0 \right]}\left( \omega _{n}^{2k} \right) +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n}^{2k} \right) \]

此时，\(\omega _{n}^{2k}\)可以使用折半引理，变成\(\omega _{n/2}^{k}\)。但需要小心：有一半的k是满足\(k\geqslant n/2\)的，这些比较大的k如果直接使用折半引理的话，会变得不好处理，所以，需要分两种情况讨论：

\[\text{对于：}A\left( \omega _{n}^{k} \right) =A^{\left[ 0 \right]}\left( \omega _{n}^{2k} \right) +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n}^{2k} \right) \\ 1. \text{当}k<\frac{n}{2}\text{时：} \\ A\left( \omega _{n}^{k} \right) =A^{\left[ 0 \right]}\left( \omega _{n}^{2k} \right) +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n}^{2k} \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right) +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right) \left( \text{由折半引理} \right) \\ 2. \text{当}k\geqslant \frac{n}{2}\text{时：} \\ \text{设：}k=k^{\prime}+n/2 \\ \text{则：}A\left( \omega _{n}^{k} \right) =A\left( \omega _{n}^{k^{\prime}+n/2} \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{2k^{\prime}+n} \right) +\omega _{n}^{k^{\prime}+n/2}A^{\left[ 1 \right]}\left( \omega _{n}^{2k^{\prime}+n} \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{2k^{\prime}} \right) +\omega _{n}^{k^{\prime}+n/2}A^{\left[ 1 \right]}\left( \omega _{n}^{2k^{\prime}} \right) \left( \text{由引理}6 \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{2k^{\prime}} \right) -\omega _{n}^{k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n}^{2k^{\prime}} \right) \left( \text{由引理}7 \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k^{\prime}} \right) -\omega _{n}^{k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k^{\prime}} \right) \left( \text{由折半引理} \right) \]

对比这两种情况得到的两个结果，可以发现：这两个结果唯一的差别就是一个正负号。也就是说，每当在\(k<n/2\)这半边算出一个点，就可以立即在\(k+n/2\)处再算出一个点。这样一来，计算量直接减少一半。这还没完，在计算\(A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right)\)和\(A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right)\)的时候，这两个多项式不仅长度减半，而且可以继续用这个性质，计算量继续减少一半。这个过程什么时候停下来呢？不难想象：当一个只有两项的多项式被拆成两个只有一项的多项式时，这个过程就停了，因为只有一项的多项式就是一个常数，其值不仅不需要计算，而且与自变量的取值无关。

上述算法的时间复杂度可由主定理证明：为\(\varTheta(NlogN)\)。这就是快速傅里叶变换的数学原理。

6.2 一个手工计算快速傅里叶变换的实例

上一节中，我们研究了快速傅里叶变换的数学原理，但如果只看公式的话，读者可能仍然不理解这个算法是怎么进行的。所以，这一节中，我们就通过一个例子，来真正计算一次快速傅里叶变换。

例：计算多项式\(A\left( x \right) =1+2x+3x^2+4x^3\)的快速傅里叶变换。

按照快速傅里叶变换的算法流程，我们首先要做的是：将这个多项式的项数补齐至2的整数次幂。由于这个多项式的项数是4，其已经是2的整数次幂了，所以，不需要补齐。

接下来，需要按奇偶顺序将多项式重排并整理：

\[A\left( x \right) =1+2x+3x^2+4x^3 \\ \text{设：}A^{\left[ 0 \right]}\left( x \right) =1+3x \\ \text{且：}A^{\left[ 1 \right]}\left( x \right) =2+4x \\ \text{则：}A\left( x \right) =A^{\left[ 0 \right]}\left( x^2 \right) +xA^{\left[ 1 \right]}\left( x^2 \right) \]

此时，我们已经将一个四项的多项式分解为两个两项的多项式了。但这还不够，我们还需要将这两个多项式继续分解为四个只有一项的多项式：

\[\text{对于：}A^{\left[ 0 \right]}\left( x \right) =1+3x \\ \text{设：}A^{\left[ 0 \right] \left[ 0 \right]}\left( x \right) =1 \\ \text{且：}A^{\left[ 0 \right] \left[ 1 \right]}\left( x \right) =3 \\ \text{则：}A^{\left[ 0 \right]}\left( x \right) =A^{\left[ 0 \right] \left[ 0 \right]}\left( x^2 \right) +xA^{\left[ 0 \right] \left[ 1 \right]}\left( x^2 \right) \\ \text{对于：}A^{\left[ 1 \right]}\left( x \right) =2+4x \\ \text{设：}A^{\left[ 1 \right] \left[ 0 \right]}\left( x \right) =2 \\ \text{且：}A^{\left[ 1 \right] \left[ 1 \right]}\left( x \right) =4 \\ \text{则：}A^{\left[ 1 \right]}\left( x \right) =A^{\left[ 1 \right] \left[ 0 \right]}\left( x^2 \right) +xA^{\left[ 1 \right] \left[ 1 \right]}\left( x^2 \right) \]

此时，\(A^{\left[ 0 \right] \left[ 0 \right]}\left( x^2 \right) \text{，}A^{\left[ 0 \right] \left[ 1 \right]}\left( x^2 \right) \text{，}A^{\left[ 1 \right] \left[ 0 \right]}\left( x^2 \right) \text{，}A^{\left[ 1 \right] \left[ 1 \right]}\left( x^2 \right)\)这四个多项式都是只包含常数项的多项式，其值均与输入无关，而分别恒等于其常数项的值。

至此，多项式已经完成分解。我们开始将n次单位复数根带入，并向前倒推多项式的值：

\[\text{对于：}A^{\left[ 0 \right]}\left( x \right) =A^{\left[ 0 \right] \left[ 0 \right]}\left( x^2 \right) +xA^{\left[ 0 \right] \left[ 1 \right]}\left( x^2 \right) \\ \text{要得到}A^{\left[ 0 \right]}\left( x \right) \text{的}FFT\text{，我们就需要将}\omega _{2}^{0}\text{和}\omega _{2}^{1}\text{带入}A^{\left[ 0 \right]}\left( x \right) \text{；} \\ \text{由快速傅里叶变换的性质，我们实际上要做的是：} \\ 1. \text{计算}\alpha =A^{\left[ 0 \right] \left[ 0 \right]}\left( \omega _{1}^{0} \right) \\ 2. \text{计算}\beta =\omega _{2}^{0}A^{\left[ 0 \right] \left[ 1 \right]}\left( \omega _{1}^{0} \right) \\ 3. A^{\left[ 0 \right]}\left( \omega _{2}^{0} \right) =\alpha +\beta \text{；}A^{\left[ 0 \right]}\left( \omega _{2}^{1} \right) =\alpha -\beta \\ \text{由于}A^{\left[ 0 \right] \left[ 0 \right]}\left( \omega _{1}^{0} \right) \text{，}A^{\left[ 0 \right] \left[ 1 \right]}\left( \omega _{1}^{0} \right) \text{的值与输入无关，就是}1\text{和}3\text{，且}\omega _{2}^{0}=1\text{，所以：} \\ \alpha =1\text{；}\beta =3\text{；}A^{\left[ 0 \right]}\left( \omega _{2}^{0} \right) =4\text{；}A^{\left[ 0 \right]}\left( \omega _{2}^{1} \right) =-2 \]

同理：

\[\text{对于：}A^{\left[ 1 \right]}\left( x \right) =A^{\left[ 1 \right] \left[ 0 \right]}\left( x^2 \right) +xA^{\left[ 1 \right] \left[ 1 \right]}\left( x^2 \right) \\ \text{要得到}A^{\left[ 1 \right]}\left( x \right) \text{的}FFT\text{，我们就需要将}\omega _{2}^{0}\text{和}\omega _{2}^{1}\text{带入}A^{\left[ 1 \right]}\left( x \right) \text{；} \\ \text{由快速傅里叶变换的性质，我们实际上要做的是：} \\ 1. \text{计算}\alpha =A^{\left[ 1 \right] \left[ 0 \right]}\left( \omega _{1}^{0} \right) \\ 2. \text{计算}\beta =\omega _{2}^{0}A^{\left[ 1 \right] \left[ 1 \right]}\left( \omega _{1}^{0} \right) \\ 3. A^{\left[ 1 \right]}\left( \omega _{2}^{0} \right) =\alpha +\beta \text{；}A^{\left[ 1 \right]}\left( \omega _{2}^{1} \right) =\alpha -\beta \\ \text{由于}A^{\left[ 1 \right] \left[ 0 \right]}\left( \omega _{1}^{0} \right) \text{，}A^{\left[ 1 \right] \left[ 1 \right]}\left( \omega _{1}^{0} \right) \text{的值与输入无关，就是}2\text{和}4\text{，且}\omega _{2}^{0}=1\text{，所以：} \\ \alpha =2\text{；}\beta =4\text{；}A^{\left[ 1 \right]}\left( \omega _{2}^{0} \right) =6\text{；}A^{\left[ 1 \right]}\left( \omega _{2}^{1} \right) =-2 \]

接下来，让我们回到最开始的目标：

\[\text{对于：}A\left( x \right) =A^{\left[ 0 \right]}\left( x^2 \right) +xA^{\left[ 1 \right]}\left( x^2 \right) \\ \text{要得到}A\left( x \right) \text{的}FFT\text{，我们就需要将}\omega _{4}^{0}\text{，}\omega _{4}^{1}\text{，}\omega _{4}^{2}\text{，}\omega _{4}^{3}\text{带入}A\left( x \right) \text{；} \\ \text{由快速傅里叶变换的性质，我们实际上要做的是：} \\ 1.\text{计算}\alpha _0=A^{\left[ 0 \right]}\left( \omega _{2}^{0} \right) \\ 2.\text{计算}\beta _0=\omega _{4}^{0}A^{\left[ 1 \right]}\left( \omega _{2}^{0} \right) \\ 3.\text{计算}\alpha _1=A^{\left[ 0 \right]}\left( \omega _{2}^{1} \right) \\ 4.\text{计算}\beta _1=\omega _{4}^{1}A^{\left[ 1 \right]}\left( \omega _{2}^{1} \right) \\ 5.A\left( \omega _{4}^{0} \right) =\alpha _0+\beta _0\text{；}A\left( \omega _{4}^{2} \right) =\alpha _0-\beta _0 \\ 6.A\left( \omega _{4}^{1} \right) =\alpha _1+\beta _1\text{；}A\left( \omega _{4}^{3} \right) =\alpha _1-\beta _1 \\ \text{而：}A^{\left[ 0 \right]}\left( \omega _{2}^{0} \right) =4\text{，}A^{\left[ 1 \right]}\left( \omega _{2}^{0} \right) =6\text{，}A^{\left[ 0 \right]}\left( \omega _{2}^{1} \right) =-2\text{，}A^{\left[ 1 \right]}\left( \omega _{2}^{1} \right) =-2\text{，都是已经计算好的} \\ \text{且}\omega _{4}^{0}=1\text{，}\omega _{4}^{1}=e^{\frac{2\pi i}{4}}=\cos \frac{\pi}{2}+i\sin \frac{\pi}{2}=i \\ \text{所以：}\alpha _0=4\text{；}\beta _0=6\text{；}\alpha _1=-2\text{；}\beta _1=-2i \\ \text{所以：}A\left( \omega _{4}^{0} \right) =10\text{；}A\left( \omega _{4}^{2} \right) =-2\text{；}A\left( \omega _{4}^{1} \right) =-2-2i\text{；}A\left( \omega _{4}^{3} \right) =-2+2i \\ \text{又：}\omega _{4}^{0}=1\text{；}\omega _{4}^{1}=i\text{；}\omega _{4}^{2}=\left( \omega _{4}^{1} \right) ^2=-1\text{；}\omega _{4}^{3}=\left( \omega _{4}^{1} \right) ^3=-i \\ \text{所以，}A\left( x \right) \text{的}FFT\text{为：} \\ \left\{ \left( 1,10 \right) ,\left( i,-2-2i \right) ,\left( -1,-2 \right) ,\left( -i,-2+2i \right) \right\} \]

至此，我们就完成了一次快速傅里叶变换的计算。读者可以将4个自变量依次带入多项式，来验证结果的正确性。

6.3 过程更简略的手工计算实例

上一节中，虽然我们花了大量的篇幅来演示一次快速傅里叶变换是怎么计算的，但实际上读者可以发现：其中的大多数过程都只是为了便于读者理解而写的，当熟练掌握后，完全可以只保留真正需要计算的部分，而这部分的计算量是非常少的。从这里就能看出，快速傅里叶变换对多项式求值带来的优化。

这一节中，我们再重新算一次上面这个多项式的快速傅里叶变换；但这一次，我们省去一切不必要的过程，只保留真正需要计算的部分，看看是什么体验。

首先，我们将多项式里面的x全部省去，只留下系数：

1 2 3 4

这里的1 2 3 4表示的就是上面的\(A\left( x \right) =1+2x+3x^2+4x^3\)

接下来，考虑对项数补齐。由于这个多项式的项数是4，刚好是2的整数次幂，所以不需要补齐。

接下来，将系数按奇偶顺序分组：

1   2   3   4
1   3 | 2   4
1 | 3 | 2 | 4

这里的1 | 3 | 2 | 4依次对应着上面的\(A^{\left[ 0 \right] \left[ 0 \right]}\left( x^2 \right) \text{，}A^{\left[ 0 \right] \left[ 1 \right]}\left( x^2 \right) \text{，}A^{\left[ 1 \right] \left[ 0 \right]}\left( x^2 \right) \text{，}A^{\left[ 1 \right] \left[ 1 \right]}\left( x^2 \right)\)，由于这四个多项式的值与输入无关，所以其值分别就是1、3、2、4。

接下来，我们从第二行的系数开始，向上倒推多项式的值。第一次倒推以连续的两个系数为一组，每一组中，相邻的两个系数之间做一对计算，需要用到系数：\(\omega _{2}^{0}=1\)：

1 3 | 2 4

1 + 1 * 3 = 4
1 - 1 * 3 = -2
2 + 1 * 4 = 6
2 - 1 * 4 = -2

4 -2 6 -2

第二次倒推以连续的四个系数为一组，每一组中，第n个系数和第n+2个系数之间做一对计算，需要用到系数：\(\omega _{4}^{0}=1\text{；}\omega _{4}^{1}=i\)：

4 -2 6 -2

4 + 1 * 6 = 10
4 - 1 * 6 = -2
-2 + i * -2 = -2-2i
-2 - i * -2 = -2+2i

10 -2-2i -2 -2+2i

至此，快速傅里叶变换就计算完成了（n次单位复数根的值不需要算出）。由此可见，一旦读者熟练掌握了快速傅里叶变换的计算原理，就可以使用这种非常简洁的计算过程进行快速傅里叶变换的计算了。

第7章 快速傅里叶逆变换

上一章中，我们研究了如何在\(\varTheta(NlogN)\)的时间复杂度下进行快速傅里叶变换，而快速傅里叶变换一旦完成，就可以进行一次时间复杂度为\(\varTheta(N)\)的多项式乘法计算，再之后，我们就需要将多项式的点值表达转为系数表达了，即进行快速傅里叶逆变换。

7.1 快速傅里叶逆变换的数学原理

对于多项式：

\[A\left( x \right) =\sum_{j=0}^{n-1}{a_jx^j} \]

我们现在已经知道其全部的快速傅里叶变换结果：

\[y_i=A\left( \omega _{n}^{i} \right) =\sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{i} \right) ^j} \]

此时，我们需要构造一个新的多项式，这个多项式的系数由\(y_i\)给出，并将单位复数根全部取倒数后带入：

\[\text{设：}B\left( x \right) =\sum_{i=0}^{n-1}{y_ix^i} \\ z_k=B\left( \omega _{n}^{-k} \right) =\sum_{i=0}^{n-1}{y_i\left( \omega _{n}^{-k} \right) ^i} \]

这样做有什么用呢？接下来，将\(y_i\)展开：

\[z_k=B\left( \omega _{n}^{-k} \right) =\sum_{i=0}^{n-1}{y_i\left( \omega _{n}^{-k} \right) ^i} \\ =\sum_{i=0}^{n-1}{\left( \left( \sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{i} \right) ^j} \right) \left( \omega _{n}^{-k} \right) ^i \right)} \\ =\sum_{i=0}^{n-1}{\left( \left( \sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{j} \right) ^i} \right) \left( \omega _{n}^{-k} \right) ^i \right)} \\ =\sum_{i=0}^{n-1}{\sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{j} \right) ^i\left( \omega _{n}^{-k} \right) ^i}} \\ =\sum_{i=0}^{n-1}{\sum_{j=0}^{n-1}{a_j\left( \omega _{n}^{j-k} \right) ^i}} \\ =\sum_{j=0}^{n-1}{a_j\sum_{i=0}^{n-1}{\left( \omega _{n}^{j-k} \right) ^i}} \]

此时，需要分两种情况讨论：

\[1.\text{当}j=k\text{时：} \\ \text{此时，}\left( \omega _{n}^{j-k} \right) ^i=1^i=1 \\ \text{则：}a_k\sum_{i=0}^{n-1}{1}=na_k \\ 2.\text{当}j\ne k\text{时：} \\ \text{此时，}\sum_{i=0}^{n-1}{\left( \omega _{n}^{j-k} \right) ^i}\text{是一个首项为}\left( \omega _{n}^{j-k} \right) ^0=1\text{，公比为}\omega _{n}^{j-k}\text{的等比数列，由等比数列求和公式得：} \\ \sum_{j=0}^{n-1}{a_j\sum_{i=0}^{n-1}{\left( \omega _{n}^{j-k} \right) ^i}} \\ =\sum_{j=0}^{n-1}{a_j\cdot \frac{1\left( 1-\left( \omega _{n}^{j-k} \right) ^n \right)}{1-\omega _{n}^{j-k}}} \\ =\sum_{j=0}^{n-1}{a_j\cdot \frac{1-1^{j-k}}{1-\omega _{n}^{j-k}}}\left( \text{由引理}6 \right) \\ =\sum_{j=0}^{n-1}{a_j\cdot 0} \\ =0 \\ \text{综上：}z_k=na_k\text{，即：}a_k=\frac{z_k}{n} \]

通过这一系列操作，我们找到了一种计算\(a_k\)的算法，基于这个算法，就能从一个多项式的点值表达转为系数表达了。所以，现在只剩下最后一个问题：\(z_k\)表示的是将单位复数根的倒数（而不是单位复数根）带入一个多项式，这并不是快速傅里叶变换的标准做法。那么此时，快速傅里叶变换还能使用吗？

为了研究这个问题，让我们回到快速傅里叶变换推导过程的起点：

\[A\left( x \right) =A^{\left[ 0 \right]}\left( x^2 \right) +xA^{\left[ 1 \right]}\left( x^2 \right) \]

回顾一下：在快速傅里叶变换的推导过程中，我们是将一半的\(\omega _{n}^{k}\)以及另一半的\(\omega _{n}^{k^{\prime}+n/2}\)带入到\(A\left( x \right)\)中的，得到的结论是：

\[1.\text{当}k<\frac{n}{2}\text{时：} \\ A\left( \omega _{n}^{k} \right) =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right) +\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right) \\ 2.\text{当}k\geqslant \frac{n}{2}\text{时：} \\ A\left( \omega _{n}^{k} \right) =A\left( \omega _{n}^{k^{\prime}+n/2} \right) =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k^{\prime}} \right) -\omega _{n}^{k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k^{\prime}} \right) \]

即：只需要计算\(A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right)\)和\(\omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right)\)就行了。

那么，如果将单位复数根的倒数带入，又会怎么样呢？还是分两种情况讨论：离0比较近的一半，和离0比较远的一半：

\[\text{对于：}A\left( \omega _{n}^{-k} \right) =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k} \right) +\omega _{n}^{-k}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k} \right) \\ 1.\text{当}k<\frac{n}{2}\text{时：} \\ A\left( \omega _{n}^{-k} \right) =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k} \right) +\omega _{n}^{-k}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k} \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{-k} \right) +\omega _{n}^{-k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{-k} \right) \left( \text{由折半引理} \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{-k} \right) +\overline{\omega _{n}^{k}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{-k} \right) \left( \text{由引理}8 \right) \\ 2.\text{当}k\geqslant \frac{n}{2}\text{时：} \\ \text{设}k=k^{\prime}+n/2 \\ \text{则：}A\left( \omega _{n}^{-k} \right) =A\left( \omega _{n}^{-k^{\prime}-n/2} \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k^{\prime}-n} \right) +\omega _{n}^{-k^{\prime}-n/2}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k^{\prime}-n} \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k^{\prime}} \right) +\omega _{n}^{-k^{\prime}+n/2}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k^{\prime}} \right) \left( \text{由引理}6 \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n}^{-2k^{\prime}} \right) -\omega _{n}^{-k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n}^{-2k^{\prime}} \right) \left( \text{由引理}7 \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{-k^{\prime}} \right) -\omega _{n}^{-k^{\prime}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{-k^{\prime}} \right) \left( \text{由折半引理} \right) \\ =A^{\left[ 0 \right]}\left( \omega _{n/2}^{-k^{\prime}} \right) -\overline{\omega _{n}^{k^{\prime}}}A^{\left[ 1 \right]}\left( \omega _{n/2}^{-k^{\prime}} \right) \left( \text{由引理}8 \right) \]

可以发现：如果将单位复数根的倒数带入，那么在进行快速傅里叶变换的过程中，只有一个系数的差别（\(A^{\left[ 0 \right]}\)和\(A^{\left[ 1 \right]}\)括号里面的差别可以忽略，因为分解后的多项式最终将与输入无关）。也就是说，只需要对快速傅里叶变换的计算过程做两处微小的改动，我们就可以再一次利用快速傅里叶变换去计算所有\(z_k\)的值，从而计算出所有\(a_k\)的值了。这两处改动分别为：

1. 将\(\omega _{n}^{k}\)换成\(\overline{\omega _{n}^{k}}\)
2. 将每个快速傅里叶变换的结果除以n

这就是快速傅里叶逆变换的数学原理。其时间复杂度与快速傅里叶变换一致，也是\(\varTheta(NlogN)\)。

7.2 一个手工计算快速傅里叶逆变换的实例

上一节中，我们研究了快速傅里叶逆变换的数学原理。这一节中，我们就使用上一章得到的快速傅里叶变换的计算结果10 -2-2i -2 -2+2i，进行一次快速傅里叶逆变换的手工计算，如果计算结果能够回到1 2 3 4，就说明我们的计算是正确的。

首先，考虑是否需要对项数补齐。由于项数是4，所以不需要补齐。

接下来，对系数进行分组：

10   -2-2i   -2      -2+2i
10   -2    | -2-2i   -2+2i

接下来，进行第一次倒推，需要用到系数：\(\overline{\omega _{2}^{0}}=1\)：

10 -2 | -2-2i -2+2i

10 + 1 * -2 = 8
10 - 1 * -2 = 12
-2-2i + 1 * -2+2i = -4
-2-2i - 1 * -2+2i = -4i

8 12 -4 -4i

接下来，进行第二次倒推，需要用到系数：\(\overline{\omega _{4}^{0}}=1\text{；}\overline{\omega _{4}^{1}}=-i\)：

8 12 -4 -4i

8 + 1 * -4 = 4
8 - 1 * -4 = 12
12 + -i * -4i = 8
12 - -i * -4i = 16

4 8 12 16

上面得到的4 8 12 16，分别就是\(z_0, z_1, z_2, z_3\)。最后，我们将其都除以\(n=4\)，就能得到\(a_0, a_1, a_2, a_3\)了：

4  / 4 = 1
8  / 4 = 2
12 / 4 = 3
16 / 4 = 4

可见，计算结果完全符合预期。

第8章 快速傅里叶变换的实现

这一章中，我们研究快速傅里叶变换以及快速傅里叶逆变换的实现，有了这两种变换，就能实现出时间复杂度为\(\varTheta(NlogN)\)的卷积算法。

8.1 对齐至2的整数次幂的算法

快速傅里叶变换的第一步是将多项式的项数对齐至2的整数次幂，所以，我们需要根据输入的项数，来找到需要对齐到的项数。这一需求的朴素算法是使用一个循环，并使用一个从1开始，不断自乘2的数字和输入项数作比较，直至这个数字已经大于等于输入项数时，算法终止。这个算法很简单，读者可以自行尝试。

这里给出一种更为高效的算法：

unsigned __nextPow2(unsigned N)
{
    N--;

    N |= N >> 1;
    N |= N >> 2;
    N |= N >> 4;
    N |= N >> 8;
    N |= N >> 16;

    return N + 1;
}

这个算法不难理解，其要点在于：

1. 如果N不是2的整数次幂，那么，就将N从最高位的1开始，到最低位之间的所有位都变成1；然后，将这个全是1的数字再加1，这些1就都会变成0，并且一个新的1将出现在原最高位的更高一位上，这正是我们需要的数字。例如0b101，我们希望将其变成0b111，再加1，就得到了0b1000，这就是我们需要的数字。那么，具体要怎么操作，才能将0b101变成0b111呢？我们可以从N的最高位的那个1开始，将其右移1位后与N位或，此时，N的最高两位就一定都是1了；接下来，将N右移两位后与N位或，使N的最高4位都变成1；以此类推：接下来使N的最高8、16、32位都变成1。读者在理解这段话时要清楚：这里所说的最高n位，都是在N足够大，确实有这么多位的前提下才成立，否则，N就会因为过多的右移而位或到一个0
2. 如果N已经是2的整数次幂，那么，算法直接返回N就行了。但是，判断一个数字是不是2的整数次幂需要额外的代价，且很明显，这个判定的失败率是很高的，因为绝大多数的整数都不是2的整数次幂。所以，干脆就不要判定这件事了，而是将N减去1。如果N是2的整数次幂，减去1后就会丢失其最高位的1，并变成一个全是1的数字，在经过多次（无用的）位或运算后，又被加上1，回到了原值。例如0b1000，减去1后会变成0b111，其丢失了原数字最高位的1，最终，0b111又会因为加1而回到0b1000并返回。另一方面，如果N并不是2的整数次幂，那就说明N除了最高位的1以外，在低位还有1，这样一来，减去1就不会使N丢失最高位的1，所以，其结果不受影响

8.2 位逆序算法

快速傅里叶变换的第二步是对系数进行分组，分组操作的朴素实现和前面的手工计算过程是一致的，读者可以自行尝试。这里给出的是一种更为高效的算法。

仔细观察分组前后，各个系数的索引值的二进制表示，这里以8个系数为例：

000 001 010 011 100 101 110 111  // 分组前
000 100 010 110 001 101 011 111  // 分组后

不难发现：分组前后的每一对索引值都是位逆序的。

这就意味着，对于输入的每一个系数，我们都可以立即知道这个系数在分组后被放在哪里了：只需要将系数的索引值进行位逆序即可。

那么，怎么实现位逆序呢？朴素的算法是：通过一个循环，将待转换的数字不断右移1位，同时将转换后的数字不断左移1位，并将两个数字的最低位对接即可。读者可以自行尝试。

这里给出的是一种更为高效的算法：

unsigned __bitReverse(unsigned N, unsigned bitWidth)
{
    N = ((0xaaaaaaaa & N) >> 1) | ((0x55555555 & N) << 1);
    N = ((0xcccccccc & N) >> 2) | ((0x33333333 & N) << 2);
    N = ((0xf0f0f0f0 & N) >> 4) | ((0x0f0f0f0f & N) << 4);
    N = ((0xff00ff00 & N) >> 8) | ((0x00ff00ff & N) << 8);

    N = ((N >> 16) | (N << 16)) >> (32 - bitWidth);

    return N;
}

这个算法不难理解，其要点在于：

1. 0xaaaaaaaa是形如0b1010...的位掩码，这个位掩码会保留N的所有奇数位；而0x55555555是形如0b0101的位掩码，这个位掩码会保留N的所有偶数位；将二者的掩码结果一个左移，一个右移，最后再位或到一起，就能使N的所有相邻位发生交换
2. 类似的，0xcccccccc是形如0b1100...的位掩码；而0x33333333是形如0b0011的位掩码；在这两个位掩码，以及后续的左右移位和位或的作用下，N会以每两位为一组发生交换。以此类推，N又会以每4、8位为一组进行交换
3. 最后，N需要以每16位为一组，完成最后一次交换，此时就不需要位掩码了，直接交换即可。至此，N的所有位完成了逆序
4. 上述算法完成的是32位无符号整数的位逆序，而实际输入的数字很可能并没有这么多位。例如：0b001的位逆序应该是0b100，而按照上述算法，最终的结果是：0b100...(后面还有29个0)，这不是我们需要的。所以，最后还需要做一次右移，将多余的0去掉

8.3 快速傅里叶变换与快速傅里叶逆变换的实现

在前面的章节中我们已经知道，快速傅里叶变换和快速傅里叶逆变换的计算过程只有两处微小的不同：

1. 快速傅里叶变换使用的系数是\(\pm \omega _{n}^{k}\)，而快速傅里叶逆变换使用的系数是\(\pm \overline{\omega _{n}^{k}}\)
2. 快速傅里叶逆变换需要在最后对所有的结果除以n

我们可以使用一个要么是1，要么是-1的数字来同时区别这两处不同。请看：

vector<complex<double>> __FFT(const vector<complex<double>> &coefList, double conjNum)
{
    vector<complex<double>> FFTList(coefList.size());

    unsigned bitWidth = __builtin_ctz(FFTList.size());

    for (unsigned idx = 0; idx < FFTList.size(); idx++)
    {
        FFTList[idx] = coefList[__bitReverse(idx, bitWidth)];
    }

    for (unsigned N = 2; N <= FFTList.size(); N *= 2)
    {
        for (unsigned startIdx = 0; startIdx < FFTList.size(); startIdx += N)
        {
            complex<double> curOmega(1.);
            complex<double> mulOmega(cos(2 * M_PI / N), conjNum * sin(2 * M_PI / N));

            for (unsigned leftIdx = startIdx, rightIdx = startIdx + N / 2; leftIdx < startIdx + N / 2; leftIdx++, rightIdx++)
            {
                auto leftNum  = FFTList[leftIdx] + curOmega * FFTList[rightIdx];
                auto rightNum = FFTList[leftIdx] - curOmega * FFTList[rightIdx];

                FFTList[leftIdx]  = leftNum;
                FFTList[rightIdx] = rightNum;

                curOmega *= mulOmega;
            }
        }
    }

    if (conjNum == -1.)
    {
        for (auto &FFTNum: FFTList)
        {
            FFTNum /= FFTList.size();
        }
    }

    return FFTList;
}

__FFT函数用于计算快速傅里叶变换以及快速傅里叶逆变换。当conjNum = 1.时，其处于快速傅里叶变换模式；而当conjNum = -1.时，其处于快速傅里叶逆变换模式（由于这个函数不作为对外接口，所以没有对conjNum使用布尔值或枚举变量等编程手段限制其他错误的值，读者如果对此感到介意，可以自行实现一个更严谨的接口）。

形参方面，coefList为系数列表，所有的系数已经由主调函数从double类型转为了complex<double>类型，并已经进行了对齐处理；conjNum已在上文中说明，其只会传入1.或-1.。

函数中，首先进行的是系数的分组操作。在进行这一操作之前，我们需要知道位逆序所需要的位宽，这是由GCC内置函数__builtin_ctz完成的，其返回输入数字从最低位到第一个1之间的0的数量。分组操作通过循环进行，其将coefList中的系数重排至FFTList列表中。

接下来的代码是一个三重循环。

第一重循环用于遍历N的取值，N从2开始，以不断自乘2的方式递增，直至与多项式的项数一致时终止。

第二重循环用于遍历分组，startIdx存放的是当前分组的起始索引值；而分组的长度（决定了startIdx的循环增量）是恒等于N的。比如，第一次倒推时，以两个数字为一组；第二次倒推时，以四个数字为一组；以此类推。

当确定了N以及当前分组后，就可以对分组内的每一对系数进行计算了。在计算过程中，\(\omega _{n}^{k}\)或\(\overline{\omega _{n}^{k}}\)需要伴随循环而变化，具体来说，每计算一对系数，当前的\(\omega _{n}^{k}\)或\(\overline{\omega _{n}^{k}}\)就需要再乘一次\(\omega _{n}^{1}\)或\(\overline{\omega _{n}^{1}}\)；而\(\omega _{n}^{k}\)或\(\overline{\omega _{n}^{k}}\)的初始值为\(\omega _{n}^{0}\)或\(\overline{\omega _{n}^{0}}\)，均为1。代码方面，curOmega变量用于存放\(\omega _{n}^{k}\)或\(\overline{\omega _{n}^{k}}\)的当前值，其被初始化为1；而mulOmega变量用于存放\(\omega _{n}^{1}\)或\(\overline{\omega _{n}^{1}}\)，其用于在每一对系数计算完成后，将curOmega自乘一次mulOmega。

mulOmega的值基于欧拉公式：

\[\omega _{n}^{1}=e^{\frac{2\pi i}{n}}=\cos \frac{2\pi}{n}+i\sin \frac{2\pi}{n} \\ \overline{\omega _{n}^{1}}=e^{\frac{2\pi i}{n}}=\cos \frac{2\pi}{n}-i\sin \frac{2\pi}{n} \]

conjNum用于控制上式中的正负号。

第三重循环用于计算当前分组内的每一对系数。分组的长度为N，将其分为左右两半，左半边的索引值由leftIdx维护，初始化为startIdx，即当前分组的起始索引值；右半边的索引值由rightIdx维护，初始化为startIdx + N / 2，即当前分组右半边的第一个索引值；这两个索引值同步向前递增，从而访问到分组内的每一对系数。在循环体中，我们同时计算并更新一对系数，然后更新curOmega。

在这个函数的最后，实现的是快速傅里叶逆变换所需的额外操作：将每个系数都除以n。

8.4 卷积的实现

卷积的实现是前面所有准备工作的汇总。让我们先梳理一下实现思路：

1. 形参方面，卷积的输入是两个不保证等长的vector<double>
2. 在进行快速傅里叶变换之前，我们需要先准备好足够多的点来表示结果多项式，即：需要将输入的两个系数列表都用0扩充到足够的长度。那么，需要多少个系数呢？这里需要做一个简单的计算：一个长度为\(n\)的系数列表，表示的是一个\(n-1\)次多项式；而另一个长度为\(m\)的系数列表，表示的是一个\(m-1\)次多项式；这两个多项式相乘的结果是一个\(n+m-2\)次多项式；而这样的多项式一共有\(n+m-1\)项。此外，根据快速傅里叶变换的要求，系数列表的长度必须是2的整数次幂，所以，我们还需要将\(n+m-1\)这个数字对齐到2的整数次幂，作为两个系数列表扩充后的长度
3. 快速傅里叶变换需要的系数列表是vector<complex<double>>类型的，而输入的系数列表是vector<double>类型的，需要进行转换
4. 当两个系数列表都准备好后，进行两次快速傅里叶变换，将两个多项式从系数表达转为点值表达；然后，通过一个循环进行点值表达下的多项式乘法；最后，再进行一次快速傅里叶逆变换，将点值表达转为系数表达
5. 快速傅里叶逆变换的输出是vector<complex<double>>类型的系数列表，而我们最终需要的是vector<double>类型的系数列表，需要进行转换。此外，还需要舍去由于对齐到2产生的系数扩充

下面请看实现：

vector<double> calcVectorConvolution(const vector<double> &leftCoefList, const vector<double> &rightCoefList)
{
    unsigned coefSize  = leftCoefList.size() + rightCoefList.size() - 1;
    unsigned alignSize = __nextPow2(coefSize);

    vector<complex<double>> leftAlignCoefList(alignSize);
    vector<complex<double>> rightAlignCoefList(alignSize);

    for (unsigned idx = 0; idx < leftCoefList.size(); idx++)
    {
        leftAlignCoefList[idx].real(leftCoefList[idx]);
    }

    for (unsigned idx = 0; idx < rightCoefList.size(); idx++)
    {
        rightAlignCoefList[idx].real(rightCoefList[idx]);
    }

    auto leftFFTList  = __FFT(leftAlignCoefList, 1.);
    auto rightFFTList = __FFT(rightAlignCoefList, 1.);

    for (unsigned idx = 0; idx < leftFFTList.size(); idx++)
    {
        leftFFTList[idx] *= rightFFTList[idx];
    }

    auto resFFTList  = __FFT(leftFFTList, -1.);

    vector<double> resCoefList(coefSize);

    for (unsigned idx = 0; idx < resCoefList.size(); idx++)
    {
        resCoefList[idx] = resFFTList[idx].real();
    }

    return resCoefList;
}

形参方面，leftCoefList和rightCoefList是两个多项式的系数表达。

coefSize用于存放结果多项式的项数，其计算公式已经由上文讨论过；alignSize用于存放将coefSize对齐到2的整数次幂后的多项式的项数，其决定了快速傅里叶变换需要的列表长度。

接下来，使用alignSize作为长度生成leftAlignCoefList和rightAlignCoefList，这两个系数列表用于快速傅里叶变换；并将leftCoefList和rightCoefList中的系数分别放入这两个列表的前面部分。

接下来，进行两次快速傅里叶变换，将leftAlignCoefList和rightAlignCoefList从系数表达转为点值表达，存放在leftFFTList和rightFFTList中；然后，使用一个循环进行点值表达下的多项式乘法，将rightFFTList乘入leftFFTList中；最后，再进行一次快速傅里叶逆变换，将leftFFTList从点值表达转为系数表达，存放在resFFTList中。

现在，resFFTList中存放的是alignSize个complex<double>类型的系数，而我们需要的是这个列表中前coefSize个double类型的系数，所以，函数的最后一段用于提取这部分系数并返回。读者可以自行验证：在resFFTList中，除了我们提取出的部分外，其余部分无论是实部还是虚部，都是0。

第9章 讨论

本文中，对时间复杂度的描述多次使用了\(\varTheta\)符号而非\(O\)符号，读者应予以关注。

第三章中研究的欧拉公式的推导过程不能作为其证明过程。这是因为，麦克劳林级数展开需要求出函数的高阶导数，而复数域下的指数函数以及三角函数的导数均依赖于欧拉公式，从而造成循环论证。欧拉公式的严格证明超出了本文的范围。

快速傅里叶变换的朴素实现基于：\(A\left( \omega _{n}^{k} \right) =A^{\left[ 0 \right]}\left( \omega _{n/2}^{k} \right) \pm \omega _{n}^{k}A^{\left[ 1 \right]}\left( \omega _{n/2}^{k} \right)\)这一结论；这是一个递归版本的算法，读者可以自行尝试。

已经存在不要求n为2的整数次幂的快速傅里叶变换算法，但其超出了本文的范围。

本文中使用的"每次计算一对系数"的操作，在相关书籍和文献中被称为蝴蝶操作（Butterfly Operation；《算法导论》中将此术语拼写为Bufferfly Operation，似为勘误）；且\(\pm \omega _{n}^{k}\)和\(\pm \overline{\omega _{n}^{k}}\)被称为旋转因子（Twiddle Factor）。本文作者认为这两个术语不够生动形象，故未在正文中引入。

快速傅里叶变换在自然科学，计算机科学等诸多领域都有着广泛的应用，希望本文能够为读者提供帮助。

樱雨楼

2023.2

标签：right,变换,text,omega,傅里叶,多项式,快速,left
From： https://www.cnblogs.com/yingyulou/p/17138800.html