快速傅里叶变换(FFT)

快速傅里叶变换 (FFT)

• 快速傅里叶变换 (FFT)
  • 1. 前置知识
    • 1.1. 复数
      • 1.1.1. 复数的引入
      • 1.1.2. 复数的几何意义
      • 1.1.3. 复数的运算
        • 1.1.3.1. 复数的加法与减法
        • 1.1.3.2. 复数的乘法与除法
        • 1.1.3.3. 复数运算与复平面上的向量的对应关系
    • 1.2. 欧拉公式
    • 1.3. 多项式
      • 1.3.1. 多项式的表示
      • 1.3.2. 朴素的多项式乘法
  • 2. 算法介绍
  • 3. 原理概述
  • 4. 原理详解
    • 4.1. 傅里叶变换
      • 4.1.1. 离散傅里叶变换/逆变换的计算公式
    • 4.2. 多项式乘法与快速傅里叶变换
    • 4.3. 离散傅里叶变换计算式的证明
      • 4.3.1. 单位复根
      • 4.3.2. 离散傅里叶变换计算式
    • 4.4. 使用离散傅里叶变换计算多项式乘法
    • 4.5. 使用快速傅里叶变换加速对离散傅里叶变换计算式的计算
      • 4.5.1. 快速傅里叶变换
      • 4.5.2. 代码实现
      • 4.5.3. 快速傅里叶变换的优化
        • 4.5.3.1. 蝴蝶变换
        • 4.5.3.2. "三步变两步"优化
  • 5. 代码模板
  • 6. 相关题目
    • 6.1. 使用快速傅里叶变换加速计算大整数乘法

1. 前置知识

1.1. 复数

1.1.1. 复数的引入

从方程的角度看，负实数能不能开平方，就是方程 $x^2+a=0(a>0)$ 有没有解，进而可以归结为方程 ${\mathbf{x}}^2+1=0$ 有没有解。

回顾已有的数集扩充过程，可以看到，每次扩充都与实际需求密切相关。例如，为了解决正方形对角线的度量，以及 $x^2-2=0$ 这样的方程在有理数集中无解的问题，人们把有理数集扩充到了实数集。数集扩充后，在实数集中规定的加法运算、乘法运算，与原来在有理数集中规定的加法运算、乘法运算协调一致，并且加法和乘法都满足交换律和结合律，乘法对加法满足分配律。

依照这种思想，为了解决 $x^2+1=0$ 这样的方程在实数系中无解的问题，我们设想引入一个新数 $i$，使得 $x=i$ 是方程 $x^2+1=0$ 的解，即使得 $i^2=-1$。

我们之前将有理数集扩充到实数集后，原先的运算性质都得到了保持，如实数 $1+\sqrt{3}$ 为一个有理数加上一个无理数，其仍能与实数集中的任意数进行原先的任意运算，运算结果仍为一个实数。
现在我们把新引进的数 $i$ 添加到实数集中，我们希望数 $i$ 和实数之间仍然能像实数那样进行加法和乘法运算，并希望加法和乘法都满足交换律、结合律，以及乘法对加法满足分配律。那么，实数系经过扩充后，得到的新数系由哪些数组成呢？

依照以上设想，把实数 $b$ 与 $i$ 相乘，结果记作 $a+bi$；把实数 $a$ 与 $bi$ 相加，结果记作 $a+bi$。注意到所有实数以及 $i$ 都可以写成 $\mathbf{a}+\mathbf{b}\mathbf{i}(\mathbf{a},\mathbf{b}\in \mathbf{R})$ 的形式，从而这些数都在扩充后的新数集中。

1.1.2. 复数的几何意义

我们知道了 $a+bi$ 这样类似的形式的数被称为复数，并且给出了定义和分类，我们还可以挖掘一下更深层的性质。

我们把所有实数都放在了数轴上，并且发现数轴上的点与实数一一对应。我们考虑对复数也这样处理。

首先我们定义复数相等：两个复数 $z_1=a+bi,z_2=c+di$ 是相等的，当且仅当 $a=c$ 且 $b=d$。

也就是说，我们可以用唯一的有序实数对 $(a,b)$ 表示一个复数 $z=a+bi$。这样，联想到平面直角坐标系，我们可以发现复数集与平面直角坐标系中的点集一一对应。好了，我们找到了复数的一种几何意义。

那么这个平面直角坐标系就不再一般，因为平面直角坐标系中的点具有了特殊意义——表示一个复数，所以我们把这样的平面直角坐标系称为复平面，$x$ 轴称为实轴，$y$ 轴称为虚轴。我们进一步地说：复数集与复平面内所有的点所构成的集合是一一对应的。

我们考虑到学过的平面向量的知识，发现向量的坐标表示也是一个有序实数对 $(a,b)$ ，显然，复数 $z=a+bi$ 对应复平面内的点 $Z(a,b)$，那么它还对应平面向量 $\stackrel{⟶}{\mathrm{OZ}}=(a,b)$，于是我们又找到了复数的另一种几何意义：复数集与复平面内的向量所构成的集合是一一对应的（实数 $0$ 与零向量对应）。

于是，我们由向量的知识迁移到复数上来，定义复数的模就是复数所对应的向量的模，即复数 $z=a+bi$ 的模为 $|z|=\sqrt{a^2+b^2}$。

于是为了方便，我们常把复数 $z=a+bi$ 称为点 $Z$ 或向量 $\stackrel{⟶}{\mathrm{OZ}}$，并规定相等的向量表示同一个复数。

并且由向量的知识我们发现，虚数不可以比较大小（但是实数是可以的）。

1.1.3. 复数的运算

1.1.3.1. 复数的加法与减法

我们规定，复数的加法规则如下：

设 $z_1=a+bi,z_2=c+di$，那么

$$z_1+z_2=(a+c)+(b+d)i$$

很明显，两个复数的和仍为复数。

考虑到向量的加法运算，我们发现复数的加法运算符合向量的加法运算法则，这同样证明了复数的几何意义的正确性。

同样可以验证，复数的加法满足交换律和结合律。

减法作为加法的逆运算，我们可以通过加法法则与复数相等的定义来推导出减法法则：

$$z_1-z_2=(a-c)+(b-d)i$$

这同样符合向量的减法运算。

1.1.3.2. 复数的乘法与除法

我们规定，复数的乘法规则如下：

设 $z_1=a+bi,z_2=c+di$，那么

$$\begin{array}{rl}{z}_1{z}_2& =(a+bi)(c+di)\\ & =ac+bci+adi+bd{i}^2\\ & =(ac-bd)+(bc+ad)i\end{array}$$

可以看出，两个复数相乘类似于两个多项式相乘，只需要把 $i^2$ 换成 $-1$，并将实部与虚部分别合并即可。

复数确实与多项式有关，因为复数域是实系数多项式环模掉 $x^2+1$ 生成的理想。（暂时还不理解）

复数的乘法与向量的向量积形式类似，是由于复数集是数环。

于是容易知道，复数乘法满足交换律，结合律和对加法的分配律。由于满足运算律，我们可以发现实数域中的乘法公式在复数域中同样适用。

除法运算是乘法运算的逆运算，我们可以像计算 $\frac{1}{\sqrt{3}+1}$ 那样推导一下：

$$\begin{array}{rl}\frac{a+bi}{c+di}& =\frac{(a+bi)(c-di)}{(c+di)(c-di)}\\ & =\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}i(c+di\ne 0)\end{array}$$

为了分母实数化，我们乘了一个 $c-di$，这个式子很有意义。

我们定义，当两个虚数实部相等，虚部互为相反数时，这两个复数互为共轭复数。通常记 $z=a+bi$ 的共轭复数为 $\overline{z}=a-bi$。我们可以发现，两个复数互为共轭复数，那么它们关于实轴对称。

1.1.3.3. 复数运算与复平面上的向量的对应关系

前面我们提到了复数集与复平面内的向量所构成的集合有一一对应关系。而复数加减法运算规则和它们所对应的向量的加减法的运算规则相同，因此复数的加减法和向量的加减法也有一一对应关系。

但是复数的乘法却与向量的乘法有所区别，复数的乘法运算规则与实数集的乘法运算规则相似，而向量的乘法运算规则却与它们不同，那么两个复数相乘后得到的复数所对应的向量与原来的两个复数所对应的向量有没有什么联系呢？

我们再回到刚才的乘法公式：

$$\begin{array}{rl}{z}_1{z}_2& =(a+bi)(c+di)\\ & =ac+bci+adi+bd{i}^2\\ & =(ac-bd)+(bc+ad)i\end{array}$$

即复数 $z_1=a+bi,z_2=c+di$ 相乘之后的结果 $z_3=(ac-bd)+(bc+ad)i$，它们对应的向量分别是$\overline{z_1}=(a,b),\overline{z_2}=(c,d),\overline{z_3}=(ac-bd,bc+ad)$，它们的模分别是

$$\begin{array}{rl}|\overline{z_1}|& =\sqrt{a^2+b^2},\\ |\overline{z_2}|& =\sqrt{c^2+d^2},\\ |\overline{z_3}|& =\sqrt{(ac-bd{)}^2+(bc+ad{)}^2}\\ & =\sqrt{(ac{)}^2-2abcd+(bd{)}^2+(bc{)}^2+2abcd+(ad{)}^2}\\ & =\sqrt{(ac{)}^2+(bd{)}^2+(bc{)}^2+(ad{)}^2}\\ & =\sqrt{a^2(c^2+d^2)+b^2(c^2+d^2)}\\ & =\sqrt{a^2(c^2+d^2)+b^2(c^2+d^2)}\\ & =\sqrt{(a^2+b^2)(c^2+d^2)}\\ & =|\overline{z_1}|\cdot |\overline{z_2}|\end{array}$$

因此两个复数相乘的得到的复数的模就等于这两个复数的模相乘，两个复数相乘后得到的复数对应的向量的模即为两个复数对应的向量的模相乘。

我们再设这三个复数对应的向量 $\overline{z_1},\overline{z_2},\overline{z_3}$ 与 x 轴的夹角分别为 ${\theta }_1,{\theta }_2,{\theta }_3$，则 $\sin {\theta }_1=\frac{b}{\sqrt{a^2+b^2}},\cos {\theta }_1=\frac{a}{\sqrt{a^2+b^2}},\sin {\theta }_2=\frac{d}{\sqrt{c^2+d^2}},\cos {\theta }_2=\frac{c}{\sqrt{c^2+d^2}},\sin {\theta }_3=\frac{bc+ad}{\sqrt{(a^2+b^2)(c^2+d^2)}}，\cos {\theta }_3=\frac{ac-bd}{\sqrt{(a^2+b^2)(c^2+d^2)}}$

观察到

$$\begin{gathered} \begin{array}{rl}\sin ({\theta }_1+{\theta }_2)& =\sin {\theta }_1\cos {\theta }_2+cos{\theta }_1\sin {\theta }_2\\ & =\frac{b}{\sqrt{a^2+b^2}}\cdot \frac{c}{\sqrt{c^2+d^2}}+\frac{a}{\sqrt{a^2+b^2}}\cdot \frac{d}{\sqrt{c^2+d^2}}\\ & =\frac{bc+ad}{\sqrt{(a^2+b^2)(c^2+d^2)}}\\ & =\sin {\theta }_3\end{array} \\ \begin{array}{rl}\cos ({\theta }_1+{\theta }_2)& =\cos {\theta }_1\cos {\theta }_2-\sin {\theta }_1\sin {\theta }_2\\ & =\frac{a}{\sqrt{a^2+b^2}}\cdot \frac{c}{\sqrt{c^2+d^2}}-\frac{b}{\sqrt{a^2+b^2}}\cdot \frac{d}{\sqrt{c^2+d^2}}\\ & =\frac{ac-bd}{\sqrt{(a^2+b^2)(c^2+d^2)}}\\ & =\cos {\theta }_3\end{array} \end{gathered}$$

因此可知 ${\theta }_1+{\theta }_2={\theta }_3$，即两个复数相乘之后所得到的复数对应的向量与 x 轴的夹角等于这两个复数对应的向量与 x 轴的夹角之和

于是我们得到了两个复数的相乘时，它们所对应的向量与乘得的结果对应的向量的关系：两个复数相乘，在复平面上表示为两个向量模长相乘，辐角相加

由于向量没有除法，这里不讨论复数除法与向量的关系。

1.2. 欧拉公式

先给出欧拉公式：$e^{i\theta }=\cos \theta +i\sin \theta$
欧拉公式揭示了复数指数幂与复数本身之间的联系

这个公式我们无法在不知道它的情况下有明确目的性地去推得它，它实际上是我们将实数集扩充到复数集并为复数集定义加减乘除以及复数指数幂运算之后自然而然满足的式子，在我们不知道它的情况下，有很多线索可以帮助我们去找到它，欧拉就是在观察泰勒级数时发现了这个公式，从而将复数指数幂运算与实数建立了联系，在我们知道它的情况下，我们也可以用很多方法去证明它。

困难：
我们现在已经将实数集扩充得到了复数集，且为复数规定了一系列实数所具备的加减乘除运算法则，除了这些运算，实数还规定了实数指数幂运算，我们如何应该如何规定复数指数幂运算呢？

我们将有理数集扩充为实数集时，对于无理数我们虽然无法知道其精确值，但是我们能够知道其近似值，其近似值还是有理数，能够求得一个数的无理数指数幂的无限近似值，因此无理数作为指数是很自然的，但是现在对于复数我们不知道其精确值($i=\sqrt{-1}$)，一个数的复数指数幂应该得到什么样的结果我们无法得知，因此要将复数指数幂运算纳入复数集，我们就必须将复数的复数指数幂与复数本身之间建立起联系

线索：

1. 复数的乘法与其对应的复平面向量的关系：我们前面证明了两个复数相乘，在复平面上表示为两个向量模长相乘，辐角相加，于是我们考虑这样一个复数 $\cos \theta +\sin \theta i$，这个复数对应的向量的模长为 1，辐角为 $\theta$，因此考虑其幂次显然可以得到：$(\cos \theta +\sin \theta \cdot i{)}^n=\cos (n\theta )+\sin (n\theta )\cdot i$，我们将其设为一个函数 $f(\theta )=\cos \theta +\sin \theta \cdot i$，于是有 $f^n(\theta )=f(n\theta )$，这恰好符合实数中指数函数的性质：$g(t)=x^t,g^n(t)=(x^t{)}^n=x^{nt}=g(nt)$

2. 还是考虑复数乘法与其对应复平面向量的关系，在复平面上，1 显然就是实数轴 x 轴上从 0 到 1 的一条直线，如果我们将它乘以 $(1+\frac{1}{n}\cdot i)$，它会向上旋转一定的角度并增长一定的长度，而当 n 很大时，$(1+\frac{1}{n}\cdot i)$ 的模长近似于 1，与 1 相乘对其模长基本没有影响，而其辐角却可以不断累加，而 $n\to \infty$ 时，1 乘上 $(1+\frac{1}{n}\cdot i{)}^n$ 基本上就是在围绕原点半径为 1 的圆上旋转，而我们知道 $e=\underset{n\to \infty }{\lim }(1+\frac{1}{n}{)}^n$，若引入复数指数幂，可以推得 $\underset{n\to \infty }{\lim }(1+\frac{i}{n}{)}^n=\underset{n\to \infty }{\lim }(1+\frac{i}{n}{)}^{\frac{n}{i}\cdot i}=e^i$，即 $e^i$ 代表复平面上的向量的逆时针旋转变换。(实际上 $e^{\pi i}+1=0$，即 $e^{\pi i}$ 表示向量旋转了 π 弧度，则一个复平面的向量乘以 $e^i$ 就表示这个向量逆时针旋转了 1 弧度，这也是复数指数幂的几何意义)

3. 观察 $e^x,\sin x,\cos x$ 的泰勒展开式（也是欧拉公式的证明）
   $e^x=1+x+\frac{1}{2!}{x}^2+\frac{1}{3!}{x}^3+...$
   $\sin x=x-\frac{1}{3!}{x}^3+\frac{1}{5!}{x}^5+...$
   $\cos x=1-\frac{1}{2!}{x}^2+\frac{1}{4!}{x}^4+...$
   观察到这三个式子形式的相似之处，将 $x=i\theta$ 带入 $e^x$ 可得：
   $\begin{array}{rl}{e}^{i\theta }& =1+i\theta +\frac{{(i\theta )}^2}{2!}+\frac{{(i\theta )}^3}{3!}+\frac{{(i\theta )}^4}{4!}+\frac{{(i\theta )}^5}{5!}+\frac{{(i\theta )}^6}{6!}+\frac{{(i\theta )}^7}{7!}...\\ & =1+i\theta -\frac{{\theta }^2}{2!}-\frac{i{\theta }^3}{3!}+\frac{{\theta }^4}{4!}+\frac{i{\theta }^5}{5!}-\frac{{\theta }^6}{6!}-\frac{i{\theta }^7}{7!}+...\\ & =(1-\frac{{\theta }^2}{2!}+\frac{{\theta }^4}{4!}-\frac{{\theta }^6}{6!})+(i\theta -\frac{i{\theta }^3}{3!}+\frac{i{\theta }^5}{5!}-\frac{i{\theta }^7}{7!})+...\\ & =(1-\frac{{\theta }^2}{2!}+\frac{{\theta }^4}{4!}-\frac{{\theta }^6}{6!})+i(\theta -\frac{{\theta }^3}{3!}+\frac{{\theta }^5}{5!}-\frac{{\theta }^7}{7!})+...\\ & =\cos \theta +i\sin \theta \end{array}$

外加一个求导证明：
令 $f(x)=\frac{e^{ix}}{\cos x+i\sin x}=\frac{u}{v}$，对该函数求导：
$f^{\prime }(x)=\frac{u^{\prime }v-u{v}^{\prime }}{v^2}=\frac{i{e}^{ix}(\cos x+i\sin x)-e^{ix}(-\sin x+i\cos x)}{(\cos x+i\sin x{)}^2}=0$
因此 $f(x)$ 是一个常值函数，随意带入一个值 $f(0)=1$，因此 $e^{ix}=\cos x+i\sin x$

1.3. 多项式

1.3.1. 多项式的表示

1. 系数表示法

系数表示法就是用一个多项式各个项的系数来表达这个多项式，即使用一个系数序列来表示多项式：
f ( x ) = a 0 + a 1 x + a 2 x 2 + . . . + a n − 1 x n − 1 ⟺ f ( x ) = { a 0 , a 1 , . . . , a n − 1 } f(x)=a_0+a_1x+a_2x^2+...+a_{n-1}x^{n-1} \Longleftrightarrow f(x)={\{a_0,a_1,...,a_{n-1}\}} f(x)=a0​+a1​x+a2​x2+...+an−1​xn−1⟺f(x)={a0​,a1​,...,an−1​}

2. 点值表示法

点值表示法是把这个多项式看成一个函数，从上面选取 $n$ 个点，从而利用这 $n$ 个点来唯一地表示这个函数。($n$ 个点可以唯一确定这个 $n$ 个参数的函数，即 $n$ 个方程可以解 $n$ 个未知数)

1.3.2. 朴素的多项式乘法

令 $f(x)=a_1{x}^2+b_{1}x+c_1，g(x)=a_2{x}^2+b_{2}x+c_2$
则它们的朴素乘法为：
$h(x)=f(x)\ast g(x)=a_1{x}^2\ast a_2{x}^2+a_1{x}^2\ast b_{2}x+a_1{x}^2\ast c_2+b_{1}x\ast a_2{x}^2+b_{1}x\ast b_{2}x+b_{1}x\ast c_2+c_1\ast a_2{x}^2+c_1\ast b_{2}x+c_1\ast c_2$
对于这两个三度多项式相乘，朴素乘法进行了 9 次运算，复杂度是 $O(n^2)$ ，若 n 非常大，这个复杂度会非常高

2. 算法介绍

离散傅里叶变换（Discrete Fourier Transform，缩写为 DFT），是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其离散时间傅里叶变换（DTFT）的频域采样。

FFT 是一种高效实现 DFT 的算法，称为快速傅里叶变换（Fast Fourier Transform，FFT）。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换进行加速实现，它对傅里叶变换的理论并没有新的发现，但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。快速数论变换 (NTT) 是快速傅里叶变换（FFT） 在 数论基础上的实现。

在 1965 年，Cooley 和 Tukey 发表了快速傅里叶变换算法。事实上 FFT 早在这之前就被发现过了，但是在当时现代计算机并未问世，人们没有意识到 FFT 的重要性。一些调查者认为 FFT 是由 Runge 和 König 在 1924 年发现的。但事实上高斯早在 1805 年就发明了这个算法，但一直没有发表。

快速傅里叶变换支持在 $O(n\log n)$ 的时间内计算两个 $n$ 度的多项式的乘法，比朴素的 $O(n^2)$ 算法更高效。由于两个整数的乘法也可以被当作多项式乘法，因此这个算法也可以用来加速大整数的乘法计算。

3. 原理概述

1. 离散傅里叶变换：离散傅里叶变换在频域中所用/所得到的是在离散时间傅里叶变换的频域中的采样，这样做是因为计算机无法处理无限连续域，因此只能将函数定义在有限离散点上，即进行采样来表示函数，因此离散傅里叶变换(DFT)就是为了便于计算机计算而产生的。

2. 离散傅里叶变换与多项式的关系：离散傅里叶变换/逆变换是对时域中的一般函数而言的，可以将一般的函数的采样与三角函数的积分的采样相互转换，而如果要在多项式乘法的计算中利用上离散傅里叶变换，我们就需要将一些特定的三角函数值(即在离散时间傅里叶变换的频域中的一些采样)代入多项式的未知数 x 中，使多项式在形式上满足离散傅里叶变换计算式

3. 多项式的系数表示法转点值表示法的过程与离散傅里叶变换的关系：

   • 多项式的系数表示法转为点值表示法： 我们用特定的复数(三角函数)来代入 x，使我们要计算的多项式在形式上满足离散傅里叶变换计算式，从而可以使用离散傅里叶变换计算得到时域中的若干点，即多项式的点值表示法
   • 把一个多项式的系数表示法转化为点值表示法的过程(频域->时域)，就叫做DFT(离散傅里叶变换)

4. 多项式的点值表示法转系数表示法与离散傅里叶逆变换的关系：

   • 我们已知一个多项式的点值表示法(时域)，于是相对地我们可以利用离散傅里叶逆变换公式得到其在频域中的若干样本点(系数乘以三角函数)，其这若干个样本点的系数即为多项式的系数表示法。
   • 把一个多项式的点值表示法转化为系数表示法的过程(时域->频域)，就是 IDFT(离散傅里叶逆变换)

5. 快速傅里叶变换(FFT)：FFT(快速傅里叶变换) 就是充分利用 DFT/IDFT 计算式中的指数因子所具有的对称性质和周期性质(即奇、偶、虚、实等特性)，进而求出这些短序列相应的 DFT 并进行适当组合，达到删除重复计算，减少乘法运算和简化结构的目的(实际上就是递归分治)，从而加速 IDFT/DFT 的过程。

6. 利用快速傅里叶变换加速多项式乘法：既然多项式可以通过代入特定的三角函数值满足离散傅里叶变换计算式，从系数表达法转换为点值表达法(频域->时域)，并通过离散傅里叶逆变换计算式还原，那我们可以将两个要相乘的多项式均通过此方法转换为点值表示法，然后将它们的点值表达法对应相乘($O(n)$),再将相乘后的结果还原回系数表示法，这样就得到了这两个多项式相乘的结果，其中的两次离散傅里叶变换以及一次离散傅里叶逆变换可以通过递归分治来加速，每次的复杂度为 $O(n\log n)$

4. 原理详解

4.1. 傅里叶变换

首先介绍一下不同种类的信号的傅里叶变换：

1. 连续时间周期信号:处理时间连续并且具有周期性的信号，其频域上离散，非周期。(傅里叶级数)
2. 连续时间非周期信号:处理时间连续但是不具有周期性的信号，其频域上连续，非周期。(连续傅里叶变换)
3. 离散时间非周期信号:处理时间离散，不具有周期性的信号，其频域上连续，有周期性。(离散时间傅里叶变换)
4. 离散时间周期信号:处理时间离散，具有周期性的信号，其频域上离散，有周期性。(近似：离散傅里叶变换将时域信号的采样变换为在离散时间傅里叶变换频域的采样可近似为 此种类型的信号)

结论：时域的周期性导致频域的离散性，时域的连续函数在频域形成非周期频谱；而时域的离散性造成频域的周期延拓，时域的非周期性对应于频域的连续函数形式。

离散傅里叶变换是离散时间信号的频谱在[0,2π]上的 N 点等间隔采样,也就是对序列频谱的离散化,这就是 DFT 的物理意义。

我们通常所说的傅里叶变换一般指的是连续傅里叶变换，即上面所说的第 2 种时域信号到频域信号的变换，能够将时域中平方可积的函数变换为频域中复指数函数的积分形式，即能将满足一定条件(1. 有限个间断点 2.有限个极值点 3.绝对可积)的某个函数表示成三角函数（正弦和/或余弦函数）或者它们的积分的线性组合。(就像是在分析这个函数的组成成分)

而离散傅里叶变换是傅里叶变换在时域和频域上都呈现离散的形式，它能将时域信号的采样变换为在 离散时间傅里叶变换（DTFT） 频域的采样

算法竞赛中我们常用离散傅里叶变换(DFT) 再加上递归分治来快速计算两个多项式的乘法，多项式的系数表示法转换为点值表示法(相当于从频域到时域)，相对地，把一个多项式的点值表示法转化为系数表示法的过程(从时域到频域)，就是离散傅里叶逆变换(IDFT)。(可能在许多其他地方将时域->频域的变换称为变换，频域->时域的变换称为逆变换，原理都一样，只是叫法有区别)

4.1.1. 离散傅里叶变换/逆变换的计算公式

1. 时域->频域：
   • 序列 { X n } n = 0 N − 1 {\{X_n\}}^{N-1}_{n=0} {Xn​}n=0N−1​的离散傅里叶变换为

$$x[k]=\underset{n=0}{\sum ^{N-1}}{X}_n{e}^{-i2\pi kn/N}$$

2. 频域->时域：

$$X_n=\frac{1}{N}\underset{n=0}{\sum ^{N-1}}x[k]e^{i2\pi kn/N}$$

4.2. 多项式乘法与快速傅里叶变换

回到原来的问题，计算两个 n 度多项式的乘法，如果我们能得到这两个 n 度多项式的点值表示法，那么将对应的点值相乘即可 $O(n)$ 地得到它们相乘的结果的点值表示法，然后将其还原回系数表达式即可，但是即使我们得到了点值表示法，想要将其还原回系数表示法的形式也是非常困难的，可以使用高斯消元，但是其复杂度为 $O(n^2)$，和系数表示法朴素计算的复杂度一样。

这就需要引入傅里叶变换了，一般我们讲的傅里叶变换是对连续函数而言，虽然多项式是连续函数，但是我们只能存储其离散点值样本，怎么使用傅里叶变换？(①) 所幸我们有离散时间傅里叶变换，但离散时间傅里叶变换处理时域中离散但非周期性的信号(即我们这里的多项式点值样本)，其在频域中是连续且周期性的，我们这里多项式的系数也是有限且非周期的，要如何表示其频域连续且周期性的特征呢？(②) 所幸我们还有离散傅里叶变换，我们可以将这有限个系数看做原频域的采样，由于原频域是周期性的，因此我们只需要在一个周期内均匀采样即可(后面再说明一下这里的‘均匀采样’在这个问题中的含义)，但就算我们解决了时域与频域中样本个数(点值个数与系数个数)的问题，我们的目的是求多项式的点值表达式，相乘之后，再将其还原回系数表达式，要如何让多项式的系数表达式满足离散傅里叶变换的条件呢？(③) 我们可以将一些特定的复数值代入多项式的未知数 x，使其满足离散傅里叶变换的计算式形式，从而将其和其系数的乘积视作频域内的样本点，变换到时域(点值表示法)；而点值表示法(时域)可以直接根据离散傅里叶逆变换计算式还原为系数表示法(频域)，最后一个问题，虽然我们能够通过离散傅里叶变换和逆变换的计算式来完成多项式的点值表示法和系数表示法的相互转换，但是如果朴素计算，其复杂度还是 $O(n^2)$，那么我们前面那么多的努力不是就白费了-_-||(④)，所幸，离散傅里叶变换的计算式中的复数(三角函数)有着许多良好的性质，我们可以通过递归分治来减少重复计算，从而使得其时间复杂度降低为 $O(n\log n)$，我们将这样的算法叫做快速傅里叶变换(FFT)。

再解释说明一下：

• 时域内的样本点即：时域中的函数(即原函数)的若干样本点的点值(即多项式的点值表示法中的点值)

• 频域内的样本点即：频域中该点的函数(即系数乘上复数(三角函数)，这若干样本点的系数即为多项式的系数表示法)

• 频域中的均匀采样即：

  1. 在系数表示法转点值表示法过程中将原多项式系数表示法中的 n 个系数看做多项式所表示的函数在频域中一个周期内均匀分布的样本点的系数
  2. 在点值表示法转系数表示法过程中得到这 n 个均匀分布的样本点的系数(即多项式的系数表示法中的 n 个系数)

以上，我们就解决了对多项式乘法计算使用傅里叶变换遇到的四个问题：
①、时域样本点离散(多项式的点值表示法只有有限个点，且计算机也只能存储有限个离散点)
解决方案：一般的连续傅里叶变换->离散时间傅里叶变换
②、频域样本点离散(多项式的系数表示法只有有限个系数，且计算机也只能存储有限个)
解决方案：离散时间傅里叶变换->离散傅里叶变换
③、含有未知数的多项式形式上与离散傅里叶变换计算式要求的形式不符合
解决方案：将 n 个特殊的复数代入到多项式的未知数 x 中，使其满足离散傅里叶变换的计算式
④、朴素按照离散傅里叶变换计算式计算，复杂度还是 $O(n^2)$
解决方案：使用快速傅里叶变换(FFT)，即根据离散傅里叶变换计算式中复数(三角函数)的性质，使用递归分治的思想简化计算，将复杂度降低为 $O(n\log n)$

于是现在对于两个 n 度多项式的乘法，我们可以通过使用快速傅里叶变换将它们从系数表示法转换为点值表示法，然后 $O(n)$ 相乘，再通过逆变换还原为系数表示法，一共做了三次 $O(n\log n)$ 的 FFT 外加一次 $O(n)$ 的相乘，总体时间复杂度 $O(n\log n)$

4.3. 离散傅里叶变换计算式的证明

现在我们来证明离散傅里叶变换的计算式。

再给出一次离散傅里叶变换的计算式：

1. 时域->频域：
   • 序列 { X n } n = 0 N − 1 {\{X_n\}}^{N-1}_{n=0} {Xn​}n=0N−1​的离散傅里叶变换为

$$x[k]=\underset{n=0}{\sum ^{N-1}}{X}_n{e}^{-i2\pi kn/N}$$

2. 频域->时域：

$$X_n=\frac{1}{N}\underset{k=0}{\sum ^{N-1}}x[k]e^{i2\pi kn/N}$$

4.3.1. 单位复根

这两个式子虽然看上去令人头皮发麻，但是观察一下我们会发现其复杂之处主要在于那个 $e$ 的幂次 $e^{i2\pi kn/N}$，能够看出这其中就有我们在前面学习欧拉公式时提到的 $e^i$，前面我们提到过，一个复平面的向量乘以 $e^i$ 就代表将这个向量逆时针旋转一弧度，我们将这个式子变换一下形式 $e^{i2\pi kn/N}=e^{\frac{i\cdot 2\pi }{N}\cdot kn}$，可以看出，一个复平面的向量乘以 $e^{\frac{i\cdot 2\pi }{N}}$ 实际上就表示将这个向量逆时针旋转 $\frac{1}{N}$ 圈，即 $\frac{2\pi }{N}$ 弧度，如此一来，$N$ 是常量，而剩下的两个量 $n,k$ 都是与时域频域的样本的编号相关的变量，这样这个式子也就不难理解了。

既然 $e^i$ 对应的是复平面上辐角为 1 弧度的单位向量，$e^{\frac{2\pi i}{N}}$ 对应的是复平面上辐角为 $\frac{2\pi }{N}$ 弧度的单位向量，考虑式中一般的 $e^{k\cdot \frac{2\pi i}{N}},(N\in Z^{\ast },0\le k\le N-1)$，我们显然可以得到关于这个复数的一些性质：

1. $(e^{k\cdot \frac{2\pi i}{N}}{)}^N=e^{2k\pi i}=\cos (2k\pi )+i\sin (2k\pi )=1$（旋转一圈，得到 1）
2. $e^{(k\cdot 2)\cdot \frac{2\pi i}{N\cdot 2}}=e^{\frac{2\pi i\cdot k}{N}}$（将一圈的总份数以及该向量的辐角所占的份数都乘以 2，值不变）
3. $e^{(k+\frac{N}{2})\cdot \frac{2\pi i}{N}}=e^{k\cdot \frac{2\pi i}{N}+\pi i}=e^{k\cdot \frac{2\pi i}{N}}\cdot e^{\pi i}=-e^{k\cdot \frac{2\pi i}{N}}$（旋转 $\frac{1}{2}$ 圈，变为相反数）

快速傅里叶变换就是根据此复数的这些性质对计算式进行递归分治地计算从而降低复杂度的

通过观察性质 1：$(e^{k\cdot \frac{2\pi i}{N}}{)}^N=1$，我们可以发现，这个式子不就跟方程 $x^N=1$ 一样吗，只不过以往我们解这个方程得到的解只有 1(N 为奇数)以及 ±1(N 为偶数)，而这个式子表明了它的解还可以是 $e^{k\cdot \frac{2\pi i}{N}}$，且 k 可以取 0~N-1，因此它的解一下子增多到了 N 个(1 和 -1(N 为偶数时) 会在 k 取 0 和 $\frac{N}{2}$(N 为偶数时) 时被取到)！

由此，我们给出可以如下定义：我们称 $x^n=1$ 在复数意义下的解是 $n$ 次复根。显然，这样的解有 $n$ 个，设 ${\omega }_n=e^{\frac{2\pi i}{n}}$ ，则 $x^n=1$ 的解集表示为 { ω n k ∣ k = 0 , 1 , . . . , n − 1 } \{\omega_n^k|k=0,1,...,n-1\} {ωnk​∣k=0,1,...,n−1}。我们称 ${\omega }_n$ 是 $n$ 次单位复根（the n-th root of unity）。

举个例子，当 $n=4$ 时，${\omega }_n=i$，即 $i$ 就是 $4$ 次单位复根：

当 $n=4$ 的时候，相当于把单位圆等分 $n=4$ 份。将每一份按照极角编号，那么我们只要知道 ${\omega }_4^1$（因为它的角度是相当于单位角度），就能知道 ${\omega }_4^0,{\omega }_4^1,{\omega }_4^2,{\omega }_4^3$。
${\omega }_4^0$ 恒等于 1，${\omega }_4^2$ 的角度是 ${\omega }_4^1$ 的两倍，所以 ${\omega }_4^2=({\omega }_4^1{)}^2=i^2=-1$，依次以此类推。

刚才我们提到 $e^{k\cdot \frac{2\pi i}{N}},(N\in Z^{\ast },0\le k\le N-1)$ 的性质就可以表示为：
对于任意正整数 n 和整数 k：

1. ${\omega }_n^n=1$
2. ${\omega }_n^k={\omega }_{2n}^{2k}$
3. ${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k\to {\omega }_{2n}^{k+n}=-{\omega }_{2n}^k$

4.3.2. 离散傅里叶变换计算式

再回到离散傅里叶变换的计算式，现在我们可以将它这样表示(变换和逆变换正好和之前公式的顺序相反，因为算法竞赛中常将使用离散傅里叶变换将多项式系数表示法转点值表示法的过程叫做变换，将其反过程叫做逆变换)：

1. 频域->时域(系数表示法->点值表示法)：

$$\begin{array}{rl}{X}_n& =\underset{k=0}{\sum ^{N-1}}x[k]e^{i2\pi kn/N}\\ & \Downarrow \\ X_n& =\underset{k=0}{\sum ^{N-1}}x[k]{\omega }_N^{k\cdot n}\\ & \Downarrow \\ X_j& =\underset{i=0}{\sum ^{n-1}}{x}_i{\omega }_n^{i\cdot j}(j=0,1,2,3,...,n-1)(1)\end{array}$$

2. 时域->频域(点值表示法->系数表示法)：
   • 序列 { X n } n = 0 N − 1 {\{X_n\}}^{N-1}_{n=0} {Xn​}n=0N−1​的离散傅里叶变换为

$$\begin{array}{rl}x[k]& =\frac{1}{n}\underset{n=0}{\sum ^{N-1}}{X}_n{e}^{-i2\pi kn/N}\\ & \Downarrow \\ x[k]& =\frac{1}{n}\underset{n=0}{\sum ^{N-1}}{X}_n{\omega }_N^{-k\cdot n}\\ & \Downarrow \\ x_i& =\frac{1}{n}\underset{j=0}{\sum ^{n-1}}{X}_j{\omega }_n^{-i\cdot j}(i=0,1,2,3,...,n-1)(2)\end{array}$$

我们将式(1)(2)分别写成矩阵形式：

$$\begin{array}{rl}{X}_j=& \underset{i=0}{\sum ^{n-1}}{x}_i{\omega }_n^{i\cdot j}(j=0,1,2,3,...,n-1)(1)\\ & \Downarrow \\ \left[\begin{array}{c}{X}_0\\ X_1\\ ⋮\\ X_{n-1}\end{array}\right]& =\left[\begin{array}{cccc}{\omega }_n^0& {\omega }_n^0& \cdots & {\omega }_n^0\\ {\omega }_n^0& {\omega }_n^1& \cdots & {\omega }_n^{n-1}\\ ⋮& ⋮& \ddots & ⋮\\ {\omega }_n^0& {\omega }_n^{n-1}& \cdots & {\omega }_n^{(n-1{)}^2}\end{array}\right]\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{n-1}\end{array}\right]\\ & \Downarrow \\ \left[\begin{array}{c}{X}_0\\ X_1\\ ⋮\\ X_{n-1}\end{array}\right]& =T\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{n-1}\end{array}\right](3)\end{array}$$

其中 $T$ 是 ${\omega }_n^{ij}$ 的矩阵，同理，对于式(2):

$$\begin{array}{rl}{x}_i=& \frac{1}{n}\underset{j=0}{\sum ^{n-1}}{X}_j{\omega }_n^{-i\cdot j}(i=0,1,2,3,...,n-1)(2)\\ & \Downarrow \\ \left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{n-1}\end{array}\right]& =\frac{1}{n}\left[\begin{array}{cccc}{\omega }_n^0& {\omega }_n^0& \cdots & {\omega }_n^0\\ {\omega }_n^0& {\omega }_n^{-1}& \cdots & {\omega }_n^{-(n-1)}\\ ⋮& ⋮& \ddots & ⋮\\ {\omega }_n^0& {\omega }_n^{-(n-1)}& \cdots & {\omega }_n^{-(n-1{)}^2}\end{array}\right]\left[\begin{array}{c}{X}_0\\ X_1\\ ⋮\\ X_{n-1}\end{array}\right]\\ & \Downarrow \\ \left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{n-1}\end{array}\right]& =U\left[\begin{array}{c}{X}_0\\ X_1\\ ⋮\\ X_{n-1}\end{array}\right](4)\end{array}$$

其中 $U$ 是 ${\omega }_n^{-ij}$ 的矩阵乘以 $\frac{1}{n}$

结合 (3)(4)，我们可以得到：

$$\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{n-1}\end{array}\right]=U\left[\begin{array}{c}{X}_0\\ X_1\\ ⋮\\ X_{n-1}\end{array}\right]=UT\left[\begin{array}{c}{x}_0\\ x_1\\ ⋮\\ x_{n-1}\end{array}\right](5)$$

于是我们只需要证明 (5) 式中的 $UT=I_n$ ($I_n$ 为 n 阶单位矩阵)，就可证明 (1)(2) 式成立，下面来证明：

$$\begin{array}{rl}& \because T_{ij}={\omega }_n^{ij}，U_{ij}=\frac{1}{n}{\omega }_n^{-ij}\\ & \therefore (TU{)}_{ik}=\stackrel{n-1}{\sum _{j=0}}{T}_{ij}{U}_{jk}=\frac{1}{n}\stackrel{n-1}{\sum _{j=0}}{\omega }_n^{ij}{\omega }_n^{-jk}=\frac{1}{n}\stackrel{n-1}{\sum _{j=0}}{\omega }_n^{j(i-k)}(6)\end{array}$$

令 $p=i-k$，有

$$\begin{array}{rl}\stackrel{n-1}{\sum _{j=0}}{\omega }_n^{j(i-k)}& =\stackrel{n-1}{\sum _{j=0}}{\omega }_n^{jp}\\ & ={\omega }_n^{0p}+{\omega }_n^{1p}+\cdots +{\omega }_n^{(n-1)p}\end{array}$$

当 $p\ne 0$ (即 $i\ne k$) 时，显然这是一个公比 $q={\omega }_n^p$ 的等比数列，由等比数列求和公式得到：

$$\stackrel{n-1}{\sum _{j=0}}{\omega }_n^{jp}=\frac{({\omega }_n^p{)}^n-1}{{\omega }_n^p-1}=\frac{({\omega }_n^n{)}^p-1}{{\omega }_n^p-1}=\frac{(1{)}^p-1}{{\omega }_n^p-1}=0(7)$$

而当 $p=0$ (即 $i=k$)时，有：

$$\stackrel{n-1}{\sum _{j=0}}{\omega }_n^{jp}=\stackrel{n-1}{\sum _{j=0}}{\omega }_n^0=\stackrel{n-1}{\sum _{j=0}}1=n(8)$$

综合(7)(8)，得到：

$$(TU{)}_{ik}=\frac{1}{n}\stackrel{n-1}{\sum _{j=0}}{\omega }_n^{j(i-k)}=\{\begin{array}{ll}1,& i=k\\ 0,& i\ne k\end{array}$$

$\therefore TU=I_n$ 得证，(1)(2)式成立

证毕。

4.4. 使用离散傅里叶变换计算多项式乘法

在上面，通过单位复根的引入以及单位复根的一些性质，我们完成了对离散傅里叶变换计算式的证明，现在我们再来看看已经被我们证明过的这两个式子：

1. $X_j=\underset{i=0}{\sum ^{n-1}}{x}_i{\omega }_n^{i\cdot j}(j=0,1,2,3,...,n-1)$

2. $x_i=\frac{1}{n}\underset{j=0}{\overset{n-1}{\sum}}X_j \omega_n^{-i\cdot j} :(i=0,1,2,3,…,n-1) $

如何将这两个变换式与多项式的点值表示法和系数表示法的相互转换联系起来呢？

前面我们提到过，将多项式看做函数，将一些特定的复数值带入其未知数 x 中即可使其满足变换式

观察第一个式子我们可以发现：

将多项式表示成函数，即令 $f(x)=\underset{i=0}{\sum ^{n-1}}{a}_i{x}^i$，再将 n 个 n 次单位复根 ${\omega }_n^j(j=0,1,2,...,n-1)$ 分别带入 $f(x)$，可得：

$f({\omega }_n^j)=\underset{i=0}{\sum ^{n-1}}{a}_i{({\omega }_n^j)}^i=\underset{i=0}{\sum ^{n-1}}{a}_i{\omega }_n^{i\cdot j}(j=0,1,2,...,n-1)$

再令 $f({\omega }_n^j)=X_j,a_i=x_i$，即可得到：

$X_j=\underset{i=0}{\sum ^{n-1}}{x}_i{\omega }_n^{i\cdot j}(j=0,1,2,...,n-1)$

这正好就是第一个离散傅里叶变换计算式，而带入点值的过程就是多项式系数表示法转换为点值表示法的过程，于是我们就通过将这些 n 次单位复根代入多项式的未知数 x 从而在使其满足离散傅里叶变换式的同时得到了多项式相应的点值表示法，我们可以通过此法得到需要相乘的两个多项式的点值表示法，然后进行 $O(n)$ 的对应项相乘，即可得到这两个多项式乘积的点值表示法

如何将结果还原回系数表示法呢？

我们在前面已经证明过上面使用的正变换式以及其逆变换式的正确性了，于是我们可以直接使用逆变换式(即第 2 个式子)，将得到的 n 个相乘后的点值带入第 2 个式子左边的 $X_k$ 中，从而计算得到相乘后的多项式的系数表示法，这个逆变换过程与原来的正变换过程的差别只在于将式中单位复根的指数乘以了 -1，并将最终结果乘以了 $\frac{1}{n}$

还有最后一个问题，n 次单位复根 ${\omega }_n^j$ 在计算机中如何表示呢？

我们可以将其转换为三角函数的形式：

${\omega }_n^j=e^{j\cdot \frac{2\pi i}{n}}=e^{i\cdot \frac{2\pi j}{n}}=\cos (\frac{2\pi j}{n})+i\sin (\frac{2\pi j}{n})(j=0,2,...,n-1)$

于是我们就实现了使用离散傅里叶变换对多项式乘法进行计算，但是朴素计算时间复杂度显然还是 $O(n^2)$，由于需要三次变换，甚至增加了复杂度

4.5. 使用快速傅里叶变换加速对离散傅里叶变换计算式的计算

4.5.1. 快速傅里叶变换

现在我们就可以来分析如何在原计算式基础上通过递归分治思想来实现快速傅里叶变换了。

我们前面已经提到过快速傅里叶变换是利用离散傅里叶变换计算式中的复数即 n 次单位复根的一些性质来递归分治地进行计算的，n 次单位复根的这些性质我们在前面已经证明过了，现在来回顾一下：

对于任意正整数 n 和整数 k：

1. ${\omega }_n^n=1$
2. ${\omega }_n^k={\omega }_{2n}^{2k}$
3. ${\omega }_n^{k+\frac{n}{2}}=-{\omega }_n^k\to {\omega }_{2n}^{k+n}=-{\omega }_{2n}^k$

再回到离散傅里叶变换计算式：

1. $X_j=\underset{i=0}{\sum ^{n-1}}{x}_i{\omega }_n^{i\cdot j}(j=0,1,2,3,...,n-1)$

2. $x_i=\frac{1}{n}\underset{j=0}{\sum ^{n-1}}{X}_j{\omega }_n^{-i\cdot j}(i=0,1,2,3,...,n-1)$

前面我们提到过，逆变换式相比变换式只是将 n 次单位复根的指数乘以了-1，并将最终结果乘以了 $\frac{1}{n}$，其余的计算过程均一致，因此这里我们只单独考虑第 1 个式子，对第 2 个式子同样适用：

我们将式子展开：

$$\begin{array}{rl}{X}_j& =x_0{\omega }_n^{0\cdot j}+x_1{\omega }_n^{1\cdot j}+x_2{\omega }_n^{2\cdot j}+x_3{\omega }_n^{3\cdot j}+...+x_{n-1}{\omega }_n^{(n-1)\cdot j}(j=0,1,2,3,...,n-1)\\ & \Downarrow \\ X_j& =x_0({\omega }_n^j{)}^0+x_1({\omega }_n^j{)}^1+x_2({\omega }_n^j{)}^2+x_3({\omega }_n^j{)}^3+...+x_{n-1}({\omega }_n^j{)}^{n-1}\end{array}$$

不妨设 $n$ 为偶数，我们将 ${\omega }_n^j$ 按其指数的奇偶性分成两组：

$$\begin{array}{rl}{X}_j& =x_0({\omega }_n^j{)}^0+x_1({\omega }_n^j{)}^1+x_2({\omega }_n^j{)}^2+x_3({\omega }_n^j{)}^3+...+x_{n-1}({\omega }_n^j{)}^{n-1}\\ & =[x_0({\omega }_n^j{)}^0+x_2({\omega }_n^j{)}^2+x_4({\omega }_n^j{)}^4+...+x_{n-2}({\omega }_n^j{)}^{n-2}]+[x_1({\omega }_n^j{)}^1+x_3({\omega }_n^j{)}^3+x_5({\omega }_n^j{)}^5+...+x_{n-1}({\omega }_n^j{)}^{n-1}]\\ & =[x_0({\omega }_n^j{)}^0+x_2({\omega }_n^j{)}^2+x_4({\omega }_n^j{)}^4+...+x_{n-2}({\omega }_n^j{)}^{n-2}]+{\omega }_n^j[x_1({\omega }_n^j{)}^0+x_3({\omega }_n^j{)}^2+x_5({\omega }_n^j{)}^4+...+x_{n-1}({\omega }_n^j{)}^{n-2}]\\ & =[x_0({\omega }_n^{0j})+x_2({\omega }_n^{2j})+x_4({\omega }_n^{4j})+...+x_{n-2}({\omega }_n^{(n-2)j})]+{\omega }_n^j[x_1({\omega }_n^{0j})+x_3({\omega }_n^{2j})+x_5({\omega }_n^{4j})+...+x_{n-1}({\omega }_n^{(n-2)j})]\\ & =[x_0({\omega }_{\frac{n}{2}}^{0j})+x_2({\omega }_{\frac{n}{2}}^{1j})+x_4({\omega }_{\frac{n}{2}}^{2j})+...+x_{n-2}({\omega }_{\frac{n}{2}}^{\frac{(n-2)}{2}j})]+{\omega }_n^j[x_1({\omega }_{\frac{n}{2}}^{0j})+x_3({\omega }_{\frac{n}{2}}^{1j})+x_5({\omega }_{\frac{n}{2}}^{2j})+...+x_{n-1}({\omega }_{\frac{n}{2}}^{\frac{(n-2)}{2}j})]\\ & =[x_0({\omega }_{\frac{n}{2}}^j{)}^0+x_2({\omega }_{\frac{n}{2}}^j{)}^1+x_4({\omega }_{\frac{n}{2}}^j{)}^2+...+x_{n-2}({\omega }_{\frac{n}{2}}^j{)}^{\frac{n-2}{2}}]+{\omega }_n^j[x_1({\omega }_{\frac{n}{2}}^j{)}^0+x_3({\omega }_{\frac{n}{2}}^j{)}^1+x_5({\omega }_{\frac{n}{2}}^j{)}^2+...+x_{n-1}({\omega }_{\frac{n}{2}}^j{)}^{\frac{n-2}{2}}]\end{array}$$

令

$$\begin{gathered} f(x)=x_0{x}^0+x_1{x}^1+x_2{x}^2+...+x_{n-1}{x}^{n-1} \\ g(x)=x_0{x}^0+x_2{x}^1+x_4{x}^2+...+x_{n-2}{x}^{\frac{n-2}{2}} \\ h(x)=x_1{x}^0+x_3{x}^1+x_5{x}^2+...+x_{n-1}{x}^{\frac{n-2}{2}} \end{gathered}$$

则

$$\begin{array}{rl}{X}_j& =f({\omega }_n^j)=g({\omega }_{\frac{n}{2}}^j)+{\omega }_n^j\cdot h({\omega }_{\frac{n}{2}}^j)(9)\\ X_{j+\frac{n}{2}}& =f({\omega }_n^{j+\frac{n}{2}})=g({\omega }_{\frac{n}{2}}^{j+\frac{n}{2}})+{\omega }_n^{j+\frac{n}{2}}\cdot h({\omega }_{\frac{n}{2}}^{j+\frac{n}{2}})\\ & \Downarrow \\ X_{j+\frac{n}{2}}& =g({\omega }_{\frac{n}{2}}^j)-{\omega }_n^j\cdot h({\omega }_{\frac{n}{2}}^j)(10)\end{array}$$

即：

$$X_j=f({\omega }_n^j)=\{\begin{array}{ll}g({\omega }_{\frac{n}{2}}^j)+{\omega }_n^j\cdot h({\omega }_{\frac{n}{2}}^j),& j<\frac{n}{2}\\ g({\omega }_{\frac{n}{2}}^j)-{\omega }_n^j\cdot h({\omega }_{\frac{n}{2}}^j),& j\ge \frac{n}{2}\end{array}$$

于是我们就可以不断地将当前序列按奇偶项分成两个序列进行递归求解，直到当前每个序列中都只含有一个元素，此时 ${\omega }_1^0={\omega }_1^1=1$，由于每次递归每个序列中的元素个数减半，因此最多递归 ${\log }_2^n$ 次，回溯时每轮计算 n 个元素的值，每个元素的计算都只需要一次乘法和一次加减法，因此总体复杂度为 $O(n\log n)$

对于逆变换，将上述过程中所有单位复根的指数乘以-1，再将最后的计算结果乘以 $\frac{1}{n}$ 即可

具体递归过程如下图所示：

4.5.2. 代码实现

现在我们来一步一步实现这个算法。

首先，能一直递归下去的条件是 $n$ 是 $2$ 的整次幂，但是这个条件很容易满足，若多项式项数不满足这个条件，我们就为其补充系数为 0 的项即可，系数为 0 的项不影响计算。

我们可以通过计算 ${\log }_2(n+1)+1$ 来得到需要将项数补充到 2 的多少次幂：

int N = 1 << __lg(n + 1) + 1;

之所以是 n+1 而不是 n，是因为可能出现 0 次多项式，我们要避免计算 ${\log }_{2}0$

我们这里的计算大部分是在复数域上进行的，我们可以使用 <complex> 中提供的复数模板类：

typedef complex<double> comp; // 重命名一下以免麻烦

定义一些常量：

const comp I(0, 1);//辐角为90°的单位复根
const double PI = acos(-1); // 环境允许也可以使用M_PI

注意到离散傅里叶变换(DFT)与其逆变换之间之间只有很小的区别，所以我们可以只用一个函数实现两者的功能：

// 把F从频域变换到时域，或从时域变换到频域
void fft(comp F[], int N, int sgn = 1); // sgn = 1表示DFT，-1表示其逆变换

递归出口是 N=1，这时频域和时域相同（都只有一个数，且${\omega }_1^0={\omega }_1^1={\omega }_1^{-1}=1$）：

if (N == 1)
    return;

先把偶数项放到 F 的前半部分，奇数项放到后半部分：

memcpy(tmp, F, sizeof(comp) * N); // 先复制一个临时数组
for (int i = 0; i < N; i++)
    *(i % 2 ? F + i / 2 + N / 2 : F + i / 2) = tmp[i]; // 这里写得有点花哨，也可以展开成if-else

递归地对前半部分和后半部分进行变换，回溯时按轮计算每轮每个序列中每个元素的值：

fft(F, N / 2, sgn), fft(F + N / 2, N / 2, sgn);

然后按刚刚推导的公式进行计算即可：

comp *G = F, *H = F + N / 2;  // 用指针比较直观
//单位复根可以递推计算，逆变换要多乘一个-1
comp cur = 1, step = exp(2 * PI / N * sgn * I);
//使用 STL的complex可以调用exp函数求出ω_n，与使用欧拉公式来将ω_n转化为三角函数形式的复数等价

for (int k = 0; k < N / 2; k++)
{
    tmp[k] = G[k] + cur * H[k];
    tmp[k + N / 2] = G[k] - cur * H[k];
    cur *= step;//累乘单位复根
}
//当前轮元素计算完毕，将全部元素复制回原数组
memcpy(F, tmp, sizeof(comp) * N);
// 这里逆变换没有除以N，因为复数除法很慢，可以等要用到时先用real()求实部再除以N

于是现在就可以通过洛谷的模板题P3803了：

#include <iostream>
#include <complex>
#include <cstdio>
#include <cstring>
using namespace std;
inline int read()
{
    int ans = 0;
    char c = getchar();
    while (!isdigit(c))
        c = getchar();
    while (isdigit(c))
    {
        ans = ans * 10 + c - '0';
        c = getchar();
    }
    return ans;
}
typedef complex<double> comp;
const int MAXN = 1000005;
const comp I(0, 1);
const double PI = acos(-1);
comp A[MAXN * 3], B[MAXN * 3], tmp[MAXN * 3], ans[MAXN * 3]; // 数组开大一些
void fft(comp F[], int N, int sgn = 1)
{
    if (N == 1)
        return;
    memcpy(tmp, F, sizeof(comp) * N);
    for (int i = 0; i < N; i++)
        *(i % 2 ? F + i / 2 + N / 2 : F + i / 2) = tmp[i];
    fft(F, N / 2, sgn), fft(F + N / 2, N / 2, sgn);
    comp *G = F, *H = F + N / 2;
    comp cur = 1, step = exp(2 * PI / N * sgn * I);
    for (int k = 0; k < N / 2; k++)
    {
        tmp[k] = G[k] + cur * H[k];
        tmp[k + N / 2] = G[k] - cur * H[k];
        cur *= step;
    }
    memcpy(F, tmp, sizeof(comp) * N);
}
int main()
{
    int n = read(), m = read(), N = 1 << __lg(n + m + 1) + 1; // 补成2的整次幂
    for (int i = 0; i <= n; ++i)
        A[i] = read();
    for (int i = 0; i <= m; ++i)
        B[i] = read();
    fft(A, N), fft(B, N);
    for (int i = 0; i < N; ++i)
        ans[i] = A[i] * B[i];
    fft(ans, N, -1);
    for (int i = 0; i <= n + m; ++i)
        printf("%d ", int(ans[i].real() / N + 0.1)); // +0.1规避浮点误差
    return 0;
}

4.5.3. 快速傅里叶变换的优化

4.5.3.1. 蝴蝶变换

然而，这个实现方式需要把数组复制过来复制过去，不免还是有点慢，我们可以对其进行常数优化：

为了避免反复的拷贝，我们可以提前确定序列中每个元素最后的位置。把下标用二进制表示，对于 8 度多项式，模拟递归过程如下：

$$[a_{000}{a}_{001}{a}_{010}{a}_{011}{a}_{100}{a}_{101}{a}_{110}{a}_{111}]$$

对其按项的奇偶性进行一次递归划分：

$$[a_{000}{a}_{010}{a}_{100}{a}_{110}][a_{001}{a}_{011}{a}_{101}{a}_{111}]$$

再对递归后得到的两个序列再进行一次递归：

$$[a_{000}{a}_{100}][a_{010}{a}_{110}][a_{001}{a}_{101}][a_{011}{a}_{111}]$$

最后得到：

$$[a_{000}][a_{100}][a_{010}][a_{110}][a_{001}][a_{101}][a_{011}][a_{111}]$$

观察一下我们可以发现：最后得到的每个元素的下标，都是其最初时下标的二进制表示的翻转。

我们称这个变换为位逆序置换（bit-reversal permutation，国内也称蝴蝶变换），我们通过蝴蝶变换，就可以提前让序列就位，然后通过迭代而不是递归地进行求解，防止拷贝。

蝴蝶变换可以通过递推 $O(n)$ 地实现：

我们定义 $rev[i]$ 为数 $i$ 进行翻转之后得到的数，然后可以进行如下递推：

1. $rev[0]=0$

2. 对于当前数 $i$，将二进制下的 $i$ 表示为前缀 $head$ +末尾字符 $tail$ 的形式，则数 $i$ 翻转的结果即为 $head$ 翻转再右移一位之后再在其开头加上 $tail$ 之后的结果
   例如对于 $1101$，$head=0110,tail=1$，$head$ 的翻转为 $0110$，将其右移一位之后得到 $011$，将 $tail$ 加到其头部之后得到 $1011$，即为 $1101$ 的翻转

3. 对于当前数 $i$，其 $head$ 表示的数即为 $i>>1(\lfloor \frac{i}{2}\rfloor )$ ，而 $i>>1$ 翻转之后的结果 $rev[i>>1]$ 在之前已经被计算过了，因此我们只需要将 $tail$ 加到 $rev[i>>1]$ 的开头即可得到数 $i$ 翻转之后的结果

4. 若 $tail=0$，相当于不做任何处理，若 $tail=1$ ，则需要将 $tail$ 左移到最高位之后加在 $rev[i>>1]$ 上

5. 因此对 $rev[i]$ 的计算过程可用代码表示为：rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << bit) ，bit 为总位数-1( $2^x=$ 多项式补充后的总项数 N，$bit=x-1$)

还需要注意的是二进制数的翻转要么是两个数一对，要么是翻转后的数是其自身，若翻转后是其自身则不用改变位置，若是两个数一对相互翻转可以互相得到，则要注意只交换一次，防止将系数交换到它的最终位置后又将其交换回去

int bit = __lg(N) - 1;
for (int i = 0; i < N; ++i)
{
    rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << bit);
    if (i < rev[i])  // 这里要防止防止重复交换
        swap(a[i], a[rev[i]]);
}

由于每个数的位置都已经确定了，现在我们就完全可以用迭代的方式进行快速傅里叶变换的计算，代码如下：

// 非递归版FFT，确保N是2的整次幂
int rev[MAXN * 3];
const comp I(0, 1);//辐角为90°的单位复根
const double PI = acos(-1);
void fft(comp F[], int N, int sgn = 1)
{
    int bit = __lg(N)-1;
    for (int i = 0; i < N; ++i) // 蝴蝶变换
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << bit);
        if (i < rev[i])
            swap(F[i], F[rev[i]]);
    }
    for (int l = 1; l < N; l <<= 1) // 枚举合并前的区间长度
    {
        comp step = exp(sgn * PI / l * I);//当前需要累乘区间长度(2l)次单位复根(2pi/2l=pi/l)
        for (int i = 0; i < N; i += l * 2) // 遍历起始点
        {
            comp cur(1, 0);//从ω_(2l)^0开始
            for (int k = i; k < i + l; ++k)//遍历半个区间，两个两个计算区间元素
            {
                comp g = F[k], h = F[k + l] * cur;
                F[k] = g + h, F[k + l] = g - h;
                cur *= step;//累乘单位复根
            }
        }
    }
}
// 逆变换记得在外部把实部除以N

这个方法比朴素的递归 FFT 快了近一倍，也是最常用的版本。

4.5.3.2. "三步变两步"优化

对于多项式乘法，还有一个“三步变两步”优化：设 $P$ 和 $Q$ 是 实多项式，$F=P+Qi$ ，则 $F^2=P^2-Q^2+2PQi$ ，注意到我们要求的 $PQ$ 正是 $F$ 虚部的一半。这样只需要两次 FFT 就可以求出结果。

for (int i = 0; i <= n; ++i)
    A[i] = read();
for (int i = 0; i <= m; ++i)
    B[i] = read();
for (int i = 0; i <= max(n, m); ++i)
    F[i] = comp(A[i], B[i]);
fft(F, N);
for (int i = 0; i < N; ++i)
    F[i] = F[i] * F[i];
fft(F, N, -1);
for (int i = 0; i <= n + m; ++i)
    printf("%d ", int(F[i].imag() / (N * 2) + 0.1));

5. 代码模板

#include <iostream>
#include <complex>
using namespace std;
typedef complex<double> comp;
const int maxn = 1e5 + 5;
// 非递归版FFT，确保N是2的整次幂
int rev[maxn * 3];
const comp I(0, 1); //辐角为90°的单位复根
const double PI = acos(-1);
void fft(comp F[], int N, int sgn = 1)
{
    int bit = __lg(N) - 1;
    for (int i = 0; i < N; ++i) // 蝴蝶变换
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << bit);
        if (i < rev[i])
            swap(F[i], F[rev[i]]);
    }
    for (int l = 1; l < N; l <<= 1) // 枚举合并前的区间长度
    {
        comp step = exp(sgn * PI / l * I); //当前需要累乘区间长度(2l)次单位复根(2pi/2l=pi/l)
        for (int i = 0; i < N; i += l * 2) // 遍历起始点
        {
            comp cur(1, 0);                 //从ω_(2l)^0开始
            for (int k = i; k < i + l; ++k) //遍历半个区间，两个两个计算区间元素
            {
                comp g = F[k], h = F[k + l] * cur;
                F[k] = g + h, F[k + l] = g - h;
                cur *= step; //累乘单位复根
            }
        }
    }
}
// 逆变换记得在外部把实部除以N

6. 相关题目

6.1. 使用快速傅里叶变换加速计算大整数乘法

• 模板题：洛谷 P1919

对于每一个 $n$ 位的十进制数，我们可以将其看做一个 $n-1$ 次多项式 $F$，满足：

$$F(x)=a_0+a_1\times 10+a_2\times 10^2+...+a_{n-1}\times 10^{n-1}$$

于是我们就可以使用快速傅里叶变换加速大整数乘法了，要注意多项式的项数设置以及最终结果的进位处理

代码如下：

#include <iostream>
#include <complex>
#include <string>
using namespace std;
typedef complex<double> comp;
const int maxn = 1e6 + 5;
// 非递归版FFT，确保N是2的整次幂
int rev[maxn * 4];
const comp I(0, 1); //辐角为90°的单位复根
const double PI = acos(-1);
void fft(comp F[], int N, int sgn = 1)
{
    int bit = __lg(N) - 1;
    for (int i = 0; i < N; ++i) // 蝴蝶变换
    {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << bit);
        if (i < rev[i])
            swap(F[i], F[rev[i]]);
    }
    for (int l = 1; l < N; l <<= 1) // 枚举合并前的区间长度
    {
        comp step = exp(sgn * PI / l * I); //当前需要累乘区间长度(2l)次单位复根(2pi/2l=pi/l)
        for (int i = 0; i < N; i += l * 2) // 遍历起始点
        {
            comp cur(1, 0);                 //从ω_(2l)^0开始
            for (int k = i; k < i + l; ++k) //遍历半个区间，两个两个计算区间元素
            {
                comp g = F[k], h = F[k + l] * cur;
                F[k] = g + h, F[k + l] = g - h;
                cur *= step; //累乘单位复根
            }
        }
    }
}
comp A[maxn * 4], B[maxn * 4];
int ans[maxn * 4];
// 逆变换记得在外部把实部除以N
int main()
{
    string a, b;
    cin >> a >> b;
    int len1 = a.length(), len2 = b.length();
    int k1 = 0, k2 = 0;
    for (int i = len1 - 1; i >= 0; i--)
        A[k1++] = a[i] - '0';
    for (int i = len2 - 1; i >= 0; i--)
        B[k2++] = b[i] - '0';
    int N = 1;
    while (N < len1 + len2)
        N <<= 1;
    fft(A, N, 1);
    fft(B, N, 1);
    for (int i = 0; i < N; i++)
        A[i] = A[i] * B[i];
    fft(A, N, -1);
    bool flag = 0;
    int N2 = N;
    for (int i = 0; i < N2; i++)
    {
        ans[i] += A[i].real() / N + 0.1;
        if (ans[i] >= 10)
        {
            ans[i + 1] += ans[i] / 10;
            ans[i] %= 10;
        }
    }
    for (int i = N2 - 1; i >= 0; i--)
    {
        if (!flag && ans[i] != 0)
            flag = 1;
        if (flag)
            cout << ans[i];
    }
}