算法简介：无穷动！无穷

快速傅里叶变换，是一种快速计算多项式乘法的算法。它通过快速的转换多项式的表现形式，从而快速计算多项式乘法。推导过程较为繁杂，代码也晦涩难懂，所以这篇文章旨在帮助各位快速理解快速傅里叶变换的一些知识。

友情链接：

【算法介绍】拉格朗日插值法 | 祝馀宫

算法基础：狂想者，呜咽

首先说，为了方便下面的推导，我来简单说一下一会要用到的几个数学工具。

三角函数相关

首先是三角函数，三角函数是为了我们后续的复数运算推导做准备的。三角函数，想必大家都不陌生，常用的就是锐角三角函数 sin , cos , tan 。假设我们有下面的直角三角形 $Rt_\triangle ABC$ （ $\angle \theta$ 是点 $A$ 所在的角，直角的顶点是 $B$，剩下一个节点是 $C$ ）

假设我们分别要求 $\sin\theta,\cos\theta$ 和 $\tan\theta$（ $0\lt \angle\theta \lt 90^\circ$），那么我们就设 $\angle \theta$ 所对的这条边为 $b$ 一般称作对边，夹在 $\angle\theta$ 和直角之间的边叫做邻边 $a$ ，最长的一条边叫做斜边 $c$ 。那么这三个三角函数的值分别为：

然后的话，顺便就给几个常见的三角函数值罢：

角度 $\theta$ $\sin\theta$ $\cos\theta$ $\tan\theta$
$0^\circ$ $0$ $1$ $0$
$30^\circ$ $\cfrac 1 2$ $\cfrac{\sqrt 3} 2$ $\cfrac{\sqrt{3}} 3$
$45^\circ$ $\cfrac{\sqrt 2} 2$ $\cfrac{\sqrt 2} 2$ $1$
$60^\circ$ $\cfrac{\sqrt 3} 2$ $\cfrac 1 2$ $\sqrt 3$
$90^\circ$ $1$ $0$ 不定义

然后来说几个三角函数常用性质：

  1. $\tan\theta=\cfrac{b}{a}=\cfrac{b}{c}\div\cfrac{a}{c}=\cfrac{\sin\theta}{\cos\theta}$
  2. 三角函数的奇偶性：$\cos(-x)=\cos(x),\ \sin(-x)=-\sin(x),\ \tan(-x)=\cfrac{\sin(-x)}{\cos(-x)}=\cfrac{-\sin(x)}{\cos(x)}=-\tan(x)$
  3. 欧拉公式：$e^{i\theta}=\cos\theta+\operatorname i \times\sin\theta$

锐角三角函数加减法公式推导

首先来推加法公式，也就是：$\cos(\theta_1+\theta_2)$ 和 $\sin(\theta_1+\theta_2)$ 当然还有一个 $\tan(\theta_1+\theta_2)$，但是说 $\tan$ 可以从另外两个的比值得出，所以先不着急

首先掏出大宝贝欧拉公式，$e^{i(\theta_1+\theta_2)}=\cos(\theta_1+\theta_2)+\operatorname i\times\sin(\theta_1+\theta_2)$ 这里面有我们要求的部分，可以从这里入手。直接开推，不等你们了

观察第一行式子的形式和最后一行式子的形式，得出结论：

然后我们回来看这个 $\tan(\theta_1+\theta_2)$ 直接大力带入开始算：

接下来推减法公式，其实有了加法公式，减法公式就特别简单，毕竟我们都学过一个式子 $a-b=a+(-b)$，也就是说减法可以看出加法对待，我们就浅浅推一推，三角函数的奇偶性我已经写在上面了。

至此搞定了三角函数的加法公式和减法公式，这个会在下面讲复数运算的时候用到

复数及其运算

然后就是复数了。首先来说一下复数的概念，我们定义一个复数 $z=a+b\operatorname i$，其中 $\operatorname i=\sqrt{-1},\operatorname i^2=-1$ 也就是我们常说的虚数单位。在一个复数 $z=a+b\operatorname i$ 中，我们称 $a$ 为实部、$b$ 为虚部，在复平面上我们可以把一个复数看成这个坐标系上的点 $(a,b)$ 。

复数表示法

像上图这样，我们同样可以把一个复数看作极坐标系上的的一个从原点指向点 $(a,b)$ 的向量，其中这个点到达原点的距离被称为模长，向量与 $x$ 轴正半轴所夹的角被称为辐角 $\theta$ 。我们同时也可以这样表示一个复数：$z=r(\cos\theta +\operatorname i\times \sin\theta)$ 。

可以尝试推一下为什么这样表示，点 $(a,b)$ 距离横轴距离为 $b$，距离纵轴距离为 $a$，这时候我们设模长为 $r$，那么上面的式子就可以表示为 $z=r\left(\cfrac a r + \operatorname i \times \cfrac b r\right)=a+b\operatorname i$，发现和一开始表示的方式是一样的，所以这样表示是正确的。

然后就是一些相关于复数的运算法则，假设我们有两个复数 $z_1,z_2$ ，那么对于他们之间的运算我们分别有：

  • $z_1+z_2 = (a_1+b_1\operatorname i)+(a_2+b_2\operatorname i)=(a_1 +a_2) + (b_1+b_2)\operatorname i$

    这个很好证明，也没什么好说的，就是正经的交换律和结合律

  • $z_1-z_2=(a_1+b_1\operatorname i) - (a_2+b_2\operatorname i)=(a_1-a_2)+(b_1-b_2)\operatorname i$

    这个和加法一样，没什么好说的

接下来才是重头戏，乘除法。乘除法的话，我希望从两个方面来讲，几何意义和代数意义。先来讲一个总的概念，复数的乘法，在几何上来说，是模长相乘辐角相加的过程。那么除法作为逆运算，就应该是模长相除辐角相减，这是为什么，我们可以先来推一推。

复数乘法的推导

首先我们一样先从代数入手，设有两个复数 $z_1=a_1+b_1\operatorname i,z_2=a_2+b_2\operatorname i$ 。
那么这两个式子相乘的结果应该是

别忘了 $\operatorname i$ 的定义 $\sqrt{-1}$ ，也就是说 $\operatorname i^2=-1$ 。

那么代数推到完了，我们应该就剩下从几何上搞定他了，同样的，我们沿用两个复数的定义，不过这次改为另外一种表示方法 $z_1=r_1(\cos\theta_1+\sin\theta_1\operatorname i), z_2=r_2(\cos\theta_2+\sin\theta_2\operatorname i)$，在根据我们刚刚说过的，模长相乘辐角相加，那么最后的结果应该是这样：

然后这时候我们就需要用到一个关于三角函数的公式：

这时候我们再用这个加法公式把上面的集合意义的式子展开来，就可以进一步推导：

似乎式子越来越玄学了，不过别忘了我们还有三角函数的定义，如果我们把三角函数的定义带入进去的话，式子就会更加有趣：

然后把这个式子带入到刚刚的展开式，就可以得到答案了：

这个就是我们说的集合形式理解复数乘法，显而易见的，这个和代数推导出来的式子完全一样，也就是说所谓模长相乘，辐角相加也是正确的了。

然后说我们就要看看所谓复数除法怎么算了，如果时复数除复数的话，显而易见是非常不好处理的，设我们有两个复数 $z_1=a_1+b_1\operatorname i,\ z_2=a_2+b_2\operatorname i$ 那么复数除法可以表示为：

我们马上就敏锐的注意到，这个的计算似乎不是我们能化简的，因为学二次根式的时候就说过，分母里头没根号根好里头没分母才算化简完成，但是如果只是简单的上下同时平方，无法消除掉这个 $\sqrt{-1}$ ，这时候，我们灵光乍现，想到了小学时老师教的平方差公式：$(a+b)(a-b)=a^2-b^2$ 这样可以做到把 $\sqrt{-1}$ 消去，说干就干：

我们上下同时乘上的 $a_2-b_2\operatorname i$ 也是个复数，我们将其称为复数 $z_2$ 的共轭复数，记作 $\overline{z_2}$ 。那么从几何意义上来讲，复数除法也是可以推导的，也就是模长相除，辐角相减，由于已经写过了乘法，所以说这个我就写的简单一点点了：

复数除法几何推导：

同样的，我们设拥有复数 $z_1=a_1+b_1\operatorname i=r_1(\cos\theta_1+\operatorname i\times\sin\theta_1),\ z_2=a_2+b_2\operatorname i=r_2(\cos\theta_2+\operatorname i\times\sin\theta_2)$，那么可以推出下式：

以上应该就是我们这里要用到的所有复数的知识点了，接下来就是一些常见的基础知识。

多项式表示法

我们最常见的，应该是多项式的系数表示法，如果我们用 $A(x)$ 表示一个未知数为 $x$ 的多项式，并将其命名为 $A$ 。那么如果是系数表示法的话，我们可以用下面的式子来表示：

其中这个多项式的的系数列表 $(a_0,a_1,\dots,a_n)$ 就是这个多项式的系数表示法。

然后有一个特殊的表示方式，如果我们把一个多项式看作一个函数的话，也就是说对于 $A(x)=y$ 来说，平面直角坐标系上的的点 $(x,y)$ 一定是在这个多项式的函数图像上的，那么众所周知，想要确定一个常数函数（也就是 $y=a$ 其中 $a$ 为定值）我们只需要一个点；要确定一个一次函数，我们需要两个点（在函数图像上的）；二次函数的话至少需要三个点……

也就是说，想要确定一个 $k$ 次函数，就需要知道 $k+1$ 个点。这时候我们说过，如果把一个多项式看成函数的话，假设这个多项式的最高项是 $n$ 那么我们就需要 $n+1$ 个点，其实这个确定多项式的过程，就是求解这个多项式系数的线性方程的过程：

友情链接：

【算法介绍】拉格朗日插值法 | 祝馀宫

这时候用来表示计算这个多项式的若干个点 $(x_0,y_0),\dots,(x_n,y_n)$ 就是这个多项式的点值表示法。以上就是这个快速傅里叶变换的全部基础知识了，接下来就是更加玄学的快速傅里叶推导过程和优化，默哀罢。

算法推导：我赞美，即兴

那么接下来就是快速傅里叶变换的内容了，我们先从基础开始。

概念与方向

众所周知，快速傅里叶变换是用于快速计算两个多项式乘积的算法，设我们有多项式 $A(x),B(x)$ 这两个多项式的次数中较高的是 $n$（因为如果没到次数可以用系数 $0$ 补齐）且他们用系数表示法分别是数列 $(a_0,\dots,a_n)$ 和数列 $(b_0,\dots,b_n)$，那么多项式乘法可以表示为：

显而易见的我们如果这样暴力的去处理多项式乘积是没有意义的，复杂度过高达到了 $O(n^2)$ 的级别。显而易见是不好的，因为这样的话我们所有人都会，于是就有了一个高明的做法。

我们把目光看到点值表示法，如果我们设第三个多项式 $C(x)=A(x)B(x)$ 来表示这两个多项式的乘积，那么这个多项式 $C(x)$ 的每一个点都应该是：

这时候我们就可以不再单独的去求解 $C(x)$ 的系数，而是把 $A(x)B(x)$ 看作一个整体求解，这样解出来的多项式就是我们想要的 $C(x)$ 这就是我们说的点值表示法来体现多项式相乘。

不难发现这个过程是 $O(n)$ 的，也就是说，我们只需要一个快速转换系数表示法和点值表示法的算法，这就是快速傅里叶变换。如果你看过了拉格朗日插值的话，你可能会说用高斯消元解方程，但是这样的复杂度也高，而且似乎不支持点值表示向系数表示的转换，只好作罢。

快速傅里叶变换的原理

定义与性质

鱼逝，聪明的傅里叶站出来了，他发明了傅里叶离散，所谓傅里叶离散，就是使用特殊值带入多项式，使得这个值具有带入的特殊值的特殊性质，从而得到答案。

那么我说的这个特殊值，就是一个叫做单位根的东西，这个单位根是什么呢？我们用 $\omega_n$ 表示一个 $n$ 次单位根，这是一个负数，且他在复平面上表示的向量，是从原点出发，指向半径为一的圆上任意一点的，$n$ 次单位根的含义是把这个圆评分为 $n$ 分。比如说下图就是一个八次单位根：

八次单位根

这个 $\omega_8^0$ 就是横轴正半轴的 $1$ ，然后往后每加上一个 $1$ 就是辐角增加 $\frac{360^\circ}{n}$，其实用欧拉公式表示出来大概就是 $\omega_n^i=e^{\operatorname i\theta},\ \theta=\frac{360^\circ i}{n}=\frac{2\pi i}{n}$（注意这里的虚数单位 $\operatorname i$ 和变量 $i$ 的区别）

当然如果要是以几何的角度来讲，这里的 $\omega_n^k$ 其实就是模长相乘辐角相加，这里的模长都是一所以相乘没有大碍，辐角相加的话，如果我们设 $\theta=\operatorname{Arg}(\omega_n^1)=\frac{360^\circ}{n}$，这样的其实就对应了一个小扇形，其他 $\omega$ 的辐角都是由多个小扇形组成的，而他们的关系是加法，所以说我们可以得到一个经典的式子 $\left(\omega_n^1\right)^k=\omega_n^k$，可以认为这是单位根的递推式。

然后就是推导快速傅里叶变换时需要的重要性质了。

单位根性质推导

一：

  • $\omega_n^0=\omega_n^n=1$

显而易见的，这个就是显而易见的，看上面八次单位根就可以知道。

二：

  • $\omega_n^0,\omega_n^1,\omega_n^2,\cdots,w_n^{n-1}$ 互不相同

虽然是显而易见的，但是我们还是严格证明一下；

反证法：假设存在一对 $a,b\ (a\neq b)$ 使得 $\omega_n^a=\omega_n^b$ 成立，则根据单位根的定义：

上面的式子没看懂的：$\omega_n^a=(\omega_n^1)^a=(e^{\frac{2\pi\operatorname i}{n}})^a$

然后一路推下去：

打了 $(1)$ 的那句话的意思大概是这样，我们把等式两边同时取自然对数，然后的话，这里的 $2\pi\operatorname i$ 是表示角的周期性，这个应该都知道，意思就是角转一圈赚回来了，然后 $k\in \mathbb Z$ 是表示转的圈是整圈。

然后我们继续来看条件：

  1. 由于 $0 \leq a, b < n$，因此 $a - b = kn$ 。
  2. 因为 $|a - b| < n$，而 $kn$ 是 $n$ 的倍数，因此 $k$ 只能为 $0$ 。
  3. 当 $k = 0$ 时：$a=b$ 和一开始的规定不符

所以说 $\omega_n^0\sim \omega_n^{n-1}$ 是互不相同的。

三：

  • $\omega_n^2=\omega_{n/2}^1,\ \omega_n^{n/2+k}=-\omega_n^k$

这也是我们要推快速傅里叶变换时需要的重要公式，这个推导还真不难，可以带入自己推：

这个的推导没什么难度，然后就是最后一个了。

亖：

  • $\sum\limits_{k=0}^{n-1}(\omega_n^{j-i})^k=\begin{cases}0, & (i\neq j)\\n, & (i=j)\end{cases}$

为了证明 $\sum\limits{k=0}^{n-1}(\omega_n^{j-i})^k=\begin{cases}0,&(i\neq j)\\n,&(i=j)\end{cases}$ ，我们需要考虑两种情况： $i=j$ 和 $i\neq j$ 。

情况1: $i=j$

当 $i=j$ 时，我们有：

因此，求和式变为：

所以，当 $i=j$ 时，求和结果为 $n$ 。

情况2: $i\neq j$

当 $i\neq j$ 时，我们有：

此时， $\omega_n^{j-i}$ 是 $n$ -次单位根，即 $\omega_n=e^{2\pi i/n}$ 。我们知道 $n$ -次单位根满足：

因此， $(\omega_n^{j-i})^n=1$ 。

求和式 $\sum{k=0}^{n-1}(\omega_n^{j-i})^k$ 是一个等比数列的和，其中首项 $a=1$ ，公比 $r=\omega_n^{j-i}$ 。等比数列的和公式为：

代入 $a=1$ 和 $r=\omega_n^{j-i}$ ，我们得到：

由于 $(\omega_n^{j-i})^n=1$ ，上式变为：

所以，当 $i\neq j$ 时，求和结果为 $0$ 。

综上所述，这个性质是成立的。

有了单位根的定义，接下来就可以进入正式的快速傅里叶变换推导了。

计算与推导

然后就是快速傅里叶变换的推导了，首先是通过已知的 $\omega_n^0\sim\omega_n^{n-1}$ 算出所有的 $A(\omega_n^0)\sim A(\omega_n^{n-1})$ 这个是我们所谓的 快速傅里叶变换，这就是上面我画的快速傅里叶变换的第一步 FFT 。

快速傅里叶变换的原理

但是很显然的是，这一步 FFT 的复杂度显而易见是小于 $O(n^2)$ 的不然我们肯定不会去用它。而我们最常见的算法基本上都是 $O(n\log n)$ 级别的，所以我们就可以略略尝试一下每次砍半。

但是我们就会发现，如果是每次对半分的话，似乎没什么滑头：

这个东西既没有减少计算规模，也没有特殊性质可以继续往下推，所以说不可取。这里说的计算规模是指可以把 $A(x)$ 分为两个多项式来解决，但是这里的话拆出来第一个多项式的次数是 $\frac{n}{2}$，第二个多项式是 $n$ 次，相较于原问题没有减少多项式的次数，也就是说解决这个问题需要的方程数量没有减少，所以说这并不算减小问题规模。

那么如果说我们要想到把多项式次数降低的话，我们可以想到一个很经典的思路，提取公因式。这个怎么提取呢，但就次数来说，原来的多项式应该是这样的：

这样就可以做到优雅的降次，但是别忘了系数，我们只提出来了 $x$ 并没有提出或者修改任何的系数，也就是说这两边拆出来的多项式只是次数相同，系数不一样的，分别是 $(a_0,a_2,\dots,a_{n-1})$ 和 $(a_1,a_3,\dots,a_n)$ 不能混为一谈，但是这样确实减少了次数也减少了计算规模，所以这个大方向是对的。

为什么说是大方向呢？因为显而易见上面这个推导是口胡的，非常不严谨而且存在略略的错误。首先我们要明确的问题是这个多项式的次数。如果我们的多项式是 $n$ 次的话，就意味着有 $n + 1$ 项，但是我们不确定这个 $n+1$ 是奇数还是偶数，这样就不利于分治。鱼逝，我们想到的一个好方法就是把项数补齐到 $n=2^k\ (k\in \operatorname N)$ 这样的话多项式就变成了 $n-1$ 次多项式，一共有 $n$ 项，可以一直砍半直到变为 $1$ 。但是多出来的项放在那里？原来的多项式被我们补齐到了 $2^k$ 项，多出来的项似乎会影响多项式的值，这个很简单，如果我们把这些项的系数设为 $0$ 就可以了（我们说的这个 $2^k-1$ 次多项式其实是构想出来补齐用的，实际上还是原来的次数）。

解决了这个项数的问题，我们就来看看按照刚刚的思路递推会是什么样的，设我们有 $n-1$ 次多项式 $A(\omega_n^k)$，其中 $n=2^k$ ：

这里如果你觉得遇到了匪夷所思的变换，就可以去上面单位根的四条性质里面找找，都是能找到依据的。

这个就是快速傅里叶变换的半个过程了，为什么说是半个，上面的图中给出了快速傅里叶变换的全过程，这里只是把系数表示法转化为了点值表示法，接下来就是把点值两两相乘，这个很简单，大家都会。但是后面的有一个叫做 IFFT 是什么东西？这个是离散傅里叶逆变换，就是把点值表示法转化为系数表示法的一个算法。我们得出结果以后，接下来就是要回到这个系数表示法。

我们上面就说过，其实这个从点值表示法转化为系数表示法的过程，就是解方程的过程，有些人已经开始想高斯消元了，这是不必要的。一开始用单位根不仅可以方便递归快速计算，在逆变换的时候也可以起到关键作用。下面的矩阵分别是每一个单位根（当然这里是被拿来当系数用了），$a$ 是我们要求的未知数矩阵，$y$ 是每一个点对应的 $y$ 坐标。

那么显然的根据定义，$\vec y=\operatorname{DFT}_n(\vec a),\ \vec a=\operatorname{IDFT}_n(\vec y)$ 这个就是我们上面说的点值转系数，系数转点值的离散傅里叶变换。如果我们把这边的 $\omega$ 系数矩阵除到右边，那么就可以得到要求的矩阵 $\vec a$ 了。那么众所周知，我们要除以一个矩阵，其实就是乘上这个矩阵的逆矩阵。

现在的问题变成了如何求这个系数矩阵的逆矩阵。这时候我们从逆矩阵的定义出发看看，我们注意到，假设我们有矩阵 $A$，那么它的逆矩阵就是 $A^{-1}$，他们相乘的结果就是一个单位矩阵，单位矩阵就是除了 $i=j$ 的点为 $1$，其他点都是 $0$ 的矩阵。这时候敏感的人应该会发现关于这个 $i = j$ 为 $1$ 、$i\neq j$ 为 $0$ 的格式有点眼熟，没错就是我们上面说过的单位根的性质亖：

可是你一看，这里 $i=j$ 的值不是 $n$ 吗？这就是傻了，我们可以让这个矩阵整体乘上 $\frac 1 n$ 啊，这样不就是单位矩阵 $I$ 了吗？所以说我们可以定义系数矩阵 $V$ 的“逆矩阵”为 $D$，那么就会有下面的式子：

那么有脑子想一想就可以知道，系数矩阵 $V$ 真正的逆矩阵应该是 $\frac 1 n D$，所以最后只需要在刚刚的矩阵乘法的等式两边同时乘上 $\frac 1 n D$ 上。那么我们也该考虑一下如何计算 $D$ 了。从性质四下手，讨论 $E$ 矩阵每一位置应是什么：

而我们知道 $\omega_n^{kj}$ 其实是系数矩阵 $V$ 的内容，也就是说实际上的 $D_{i,j}$ 应该就是 $\omega_n^{-ik}$ 写出来是这样的：

然而即便如此我们似乎仍然不好计算这个矩阵，鱼逝详细推算一下这个 $\omega_n^{-k}$ 和 $\omega_n^k$ 什么关系：

不记得这些符号的就往上翻罢，这里最后一步用的是三角函数奇偶性，结果发现就是原来这玩意的共轭复数。这样的话大体的思路已经差不多了，所以我们就来看看快速傅里叶变换递归版本的代码罢。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
// complex 复数类型，n 长度，inv = 1/-1 i*sin的符号，决定了是 FFT 还是 IFFT
const double pi = acos(-1);
void FFT(complex<double >*a, int n, int inv) {
    if(n == 1) return ;
    int mid = n / 2;
    complex<double >a0[mid + 1], a1[mid + 1];
    for(int i = 0; i <= n; i += 2) {
        a0[i / 2] = a[i], a1[i / 2] = a[i + 1];
    }
    FFT(a0, mid, inv), FFT(a1, mid, inv);
    complex<double >w0(1,0), wn(cos(2 * pi / n), inv * sin(2 * pi / n));
    for(int i = 0; i < mid; i++, w0 *= wn) {
        a[i] = (a0[i] + w0 * a1[i]) / n;
        a[i + mid] = (a0[i] - w0 * a1[i]) / n;
    }
}

// how to use? example a*b
FFT(a, n, 1), FFT(b, n, 1);
for(int i = 0; i <= n; i++)
    c[i] = a[i] * b[i];
FFT(c, n, -1);

算法优化：把宣叙呈献给

其实的话，这里就很简单了，我们常见的快速傅里叶变换不是递归形式而是递推形式，这其实也很好理解，递归形式费力费时费空间，我们理解了递归的思想以后改用递推就可以。

首先是常数优化，快速傅里叶变换上面的递归版本也写了，我们要三次 FFT，常数巨大，于是呢我们有一种优化是可以把三次 FFT 改为两次，这个优化比较逆天，大概是下面这个意思，我们观察一个虚数的平方项结构：

然后就发现虚部（也就是含 $\operatorname i$ 项去除掉 $\operatorname i$ 的部分）刚刚好是 $a\times b$，于是就拥有了一个大胆的想法，我们让这个虚数的实部存储 $A(x)$ 虚部存储 $B(x)$ 这样的话我们只需要做一遍快速傅里叶变换后平方，然后对其进行快速傅里叶逆变换最后取虚部除以二就是答案（

另外还有一个内存上的优化是把递归改成线性结构，我们注意到一个性质，假设我们有 $16$ 个数，那么他们最终的分组应该是下图这样的：

如果你仔细观察的话，就会发现，括号里面的二进制倒过来，是相邻的数字！就比如说 $0$ 和 $8$ 他们倒过来就是 0000 和 0001 相邻的；$4$ 和 $12$ 倒过来是 0010 和 0011 也是相邻的数字。网上回溯一层到步骤三，就会发现 $0,4,8,12$ 倒过来是连续的 $0\sim 3$ 等等等等。这个步骤叫做蝴蝶变换，把他们的二进制倒过来，就可以转化为连续的区间问题。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
inline void FFT(complex<double> *a, int n, int inv) { // 参数基本同上
    for (int i = 0; i < n; ++i) //先将a变成最后的样子
        if (i < rev[i]) swap(a[i], a[rev[i]]);
    for (int i = 1; i < n; i <<= 1) { //合并区间的长度的二分之一，或者是合并区间的中点
        complex<double> wn(cos(pi / i), inv * sin(pi / i)); //原根
        for (int j = 0; j < n; j += (i << 1)) { //枚举具体区间，j也就是区间右端点
            complex<double> w0(1, 0);
            for (int k = 0; k < i; ++k, w0 *= wn) { //合并
                complex<double> x = a[j + k], y = w0 * a[i + j + k];
                a[j + k] = x + y;
                a[i + j + k] = x - y;
            }
        }
    }
}

算法应用：只有今晚奏鸣

那么也应该正经的来做一题了，模板就别想了，当然类似于模板的模板还是有的。

题目大意

题目传送门

高精度乘法，不过数据范围很大，不可以 $O(n^2)$ 解决。

题目思路

摆明了就是傅里叶变换，多项式的 $x$ 是 $10$，系数就是每一项的那个数字。最后变幻出来还要进位。

题目代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
#include<bits/stdc++.h>
using namespace std;

typedef long long ll;
// typedef __int128 int128;
typedef pair<int , int > pii;
typedef unsigned long long ull;

namespace FastIO 
{
    // char buf[1 << 20], *p1, *p2;
    // #define getchar() (p1 == p2 && (p2 = (p1 = buf) + fread(buf, 1, 1 << 20, stdin), p1 == p2) ? EOF : *p1++)
    template<typename T> inline T read(T& x) {
        x = 0;
        int f = 1;
        char ch;
        while (!isdigit(ch = getchar())) if (ch == '-') f = -1;
        while (isdigit(ch)) x = (x << 1) + (x << 3) + (ch ^ 48), ch = getchar();
        x *= f;
        return x;
    }
    template<typename T, typename... Args> inline void read(T& x, Args &...x_) {
        read(x);
        read(x_...);
        return;
    }
    inline ll read() {
        ll x;
        read(x);
        return x;
    }
};
using namespace FastIO;

const int N = (1 << 21);
const double pi = acos(-1);

complex<double >a[N];
string s;

int n;

inline void Input() {
    cin >> s; n += s.length();
    for(int i = s.length() - 1, cnt = 0; i >= 0; i--, cnt++) {
        a[cnt].real(s[i] - '0');
    }
    cin >> s; n += s.length();
    for(int i = s.length() - 1, cnt = 0; i >= 0; i--, cnt++) {
        a[cnt].imag(s[i] - '0');
    }
}

int rev[N];

inline void FFT(complex<double >*a, int n, int inv) {
    for(int i = 0; i < n; i++)
        if(i < rev[i]) swap(a[i], a[rev[i]]);
    for(int i = 1; i < n; i <<= 1) {
        complex<double> wn(cos(pi / i), inv * sin(pi / i));
        for(int j = 0; j < n; j += (i << 1)) {
            complex<double> w0(1, 0);
            for(int k = 0; k < i; k++, w0 *= wn) {
                complex<double> x = a[j + k], y = w0 * a[i + j + k];
                a[j + k] = x + y;
                a[i + j + k] = x - y;
            }
        }
    }
    for(int i = 0; i < n; i++) {
        if(inv == -1) a[i] /= n;
    }
}

int c[N];

inline void Work() {
    int len = 1;
    while(len < n) len <<= 1;
    n = len;
    for(int i = 0; i <= n; i++) {
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (max((int)ceil(log2(n)), 1) - 1));
    }
    FFT(a, n, 1);
    for(int i = 0; i <= n; i++) 
        a[i] = a[i] * a[i];
    FFT(a, n, -1);
    for(int i = 0; i <= n; i++) {
        c[i] += a[i].imag() / 2 + 0.49;
        c[i + 1] += c[i] / 10;
        c[i] %= 10;
    }
    len = N - 1;
    while(c[len] == 0 && len > 0) len--;
    for(int i = len; i >= 0; i--) {
        printf("%d", c[i]);
    }
}

int main() {
    int T = 1;
    while(T--) {
        Input();
        Work();
    }
    return 0;
}

算法实战：回旋，悄悄地

查看相关算法标签喵！

参考资料（以下排名不分先后）