学习笔记：快速傅里叶变换
tsqtsqtsq0309
·
2023-08-30 19:50:52
·
算法·理论
快速傅里叶变换
引入
快速傅里叶变换 （fast Fourier transform） ，即利用计算机计算离散傅里叶变换（DFT）的高效、快速计算方法的统称，简称 FFT。快速傅里叶变换是1965年由J.W.库利和T.W.图基提出的。采用这种算法能使计算机计算离散傅里叶变换所需要的乘法次数大为减少，特别是被变换的抽样点数N越多，FFT 算法计算量的节省就越显著。
FFT（Fast Fourier Transformation） 是离散傅氏变换（DFT）的快速算法。即为快速傅氏变换。它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。
朴素高精度乘法的时间复杂度为 O(n^2) ，而 FFT 的时间复杂度为 O(n\log n) ，在高精度乘法中优化效果明显。
复数
引入
从方程的角度看，负实数能不能开平方，就是方程 x^2+a=0 (a>0) 有没有解，进而可以归结为方程 x^2+1=0 有没有解。
回顾已有的数集扩充过程，可以看到，每次扩充都与实际需求密切相关。例如，为了解决正方形对角线的度量，以及 x^2-2=0 这样的方程在有理数集中无解的问题，人们把有理数集扩充到了实数集。数集扩充后，在实数集中规定的加法运算、乘法运算，与原来在有理数集中规定的加法运算、乘法运算协调一致，并且加法和乘法都满足交换律和结合律，乘法对加法满足分配律。
依照这种思想，为了解决 x^2+1=0 这样的方程在实数系中无解的问题，我们设想引入一个新数 \mathrm{i} ，使得 x=\mathrm{i} 是方程 x^2+1=0 的解，即使得 \mathrm{i}^2=-1 。
思考：把新引进的数 \mathrm{i} 添加到实数集中，我们希望数 \mathrm{i} 和实数之间仍然能像实数那样进行加法和乘法运算，并希望加法和乘法都满足交换律、结合律，以及乘法对加法满足分配律。那么，实数系经过扩充后，得到的新数系由哪些数组成呢？
依照以上设想，把实数 b 与 \mathrm{i} 相乘，结果记作 b\mathrm{i} ；把实数 a 与 b\mathrm{i} 相加，结果记作 a+b\mathrm{i} 。注意到所有实数以及 \mathrm{i} 都可以写成 a+b\mathrm{i}(a,b\in \mathbf{R}) 的形式，从而这些数都在扩充后的新数集中。
定义
我们定义形如 a+b\mathrm{i} ，其中 a,b\in \mathbf{R} 的数叫做 复数 ，其中 \mathrm{i} 被称为 虚数单位 ，全体复数的集合叫做 复数集 ，记作 \mathbf{C} 。
复数通常用 z 表示，即 z=a+b\mathrm{i} 。这种形式被称为 复数的代数形式 。其中 a 称为复数 z 的 实部 ，记作 \operatorname{Re}(z) ，b 称为复数 z 的 虚部 ，记作 \operatorname{Im}(z) 。如无特殊说明，都有 a,b\in \mathbf{R} 。
对于一个复数 z ，当且仅当 b=0 时，它是实数，当 b\not = 0 时，它是虚数，当 a=0 且 b\not = 0 时，它是纯虚数。
性质与运算
几何意义
我们知道了 a+b\mathrm{i} 这样类似的形式的数被称为复数，并且给出了定义和分类，我们还可以挖掘一下更深层的性质。
我们把所有实数都放在了数轴上，并且发现数轴上的点与实数一一对应。我们考虑对复数也这样处理。
首先我们定义 复数相等 ：两个复数 z_1=a+b\mathrm{i},z_2=c+d\mathrm{i} 是相等的，当且仅当 a=c 且 b=d 。
这么定义是十分自然的，在此不做过多解释。
也就是说，我们可以用唯一的有序实数对 (a,b) 表示一个复数 z=a+b\mathrm{i} 。这样，联想到平面直角坐标系，我们可以发现 复数集与平面直角坐标系中的点集一一对应 。好了，我们找到了复数的一种几何意义。
那么这个平面直角坐标系就不再一般，因为平面直角坐标系中的点具有了特殊意义——表示一个复数，所以我们把这样的平面直角坐标系称为 复平面 ，x 轴称为 实轴 ，y 轴称为 虚轴 。我们进一步地说：复数集与复平面内所有的点所构成的集合是一一对应的 。
我们考虑到学过的平面向量的知识，发现向量的坐标表示也是一个有序实数对 (a,b) ，显然，复数 z=a+b\mathrm{i} 对应复平面内的点 Z(a,b) ，那么它还对应平面向量 \overrightarrow{OZ}=(a,b) ，于是我们又找到了复数的另一种几何意义：复数集与复平面内的向量所构成的集合是一一对应的（实数 0 与零向量对应） 。
于是，我们由向量的知识迁移到复数上来，定义 复数的模 就是复数所对应的向量的模。复数 z=a+b\mathrm{i} 的模 |z|=\sqrt{a^2+b^2} 。
于是为了方便，我们常把复数 z=a+b\mathrm{i} 称为点 Z 或向量 \overrightarrow {OZ} ，并规定相等的向量表示同一个复数。
并且由向量的知识我们发现，虚数不可以比较大小（但是实数是可以的）。
加法与减法
对复数 z_1=a+b\mathrm{i},z_2=c+d\mathrm{i} ，定义加法规则如下：
z_1+z_2=(a+c)+(b+d)\mathrm{i}
很明显，两个复数的和仍为复数。
考虑到向量的加法运算，我们发现复数的加法运算符合向量的加法运算法则，这同样证明了复数的几何意义的正确性。
同样可以验证，复数的加法满足 交换律 和 结合律 。即：
\begin{aligned}
z_1+z_2&=z_2+z_1\\
(z_1+z_2)+z_3&=z_1+(z_2+z_3)
\end{aligned}
减法作为加法的逆运算，我们可以通过加法法则与复数相等的定义来推导出减法法则：
z_1-z_2=(a-c)+(b-d)\mathrm{i}
这同样符合向量的减法运算。
乘法、除法与共轭
对复数 z_1=a+b\mathrm{i},z_2=c+d\mathrm{i} ，定义乘法规则如下：
\begin{aligned}
z_1z_2&=(a+b\mathrm{i})(c+d\mathrm{i})\\
&=ac+bc\mathrm{i}+ad\mathrm{i}+bd\mathrm{i}^2\\
&=(ac-bd)+(bc+ad)\mathrm{i}
\end{aligned}
可以看出，两个复数相乘类似于两个多项式相乘，只需要把 \mathrm{i}^2 换成 -1 ，并将实部与虚部分别合并即可。
复数的乘法与向量的向量积形式类似。
易得复数乘法满足 交换律 ，结合律 和 对加法的分配律 ，即：
z_1z_2=z_2z_1
(z_1z_2)z_3=z_1(z_2z_3)
z_1(z_2+z_3)=z_1z_2+z_1z_3
由于满足运算律，我们可以发现实数域中的 乘法公式在复数域中同样适用 。
除法运算是乘法运算的逆运算，我们可以推导一下：
\begin{aligned}
\frac{a+b\mathrm{i}}{c+d\mathrm{i}}&=\frac{(a+b\mathrm{i})(c-d\mathrm{i})}{(c+d\mathrm{i})(c-d\mathrm{i})}\\
&=\frac{ac+bd}{c^2+d^2}+\frac{bc-ad}{c^2+d^2}\mathrm{i} &(c+d\mathrm{i}\not =0)
\end{aligned}
由于向量没有除法，这里不讨论与向量的关系。
为了分母实数化，我们乘了一个 c-d\mathrm{i} ，这个式子很有意义。
对复数 z=a+b\mathrm{i} ，称 a-b\mathrm{i} 为 z 的 共轭复数 ，通常记为 \bar z 。我们可以发现，若两个复数互为共轭复数，那么它们 关于实轴对称 。
对复数 z,w ，复数共轭有如下性质
z\cdot\bar{z}=|z|^2
\overline{\overline{z}}=z
\operatorname{Re}(z)=\dfrac{z+\bar{z}}{2}$，$\operatorname{Im}(z)=\dfrac{z-\bar{z}}{2}
\overline{z\pm w}=\bar{z}\pm\bar{w}
\overline{zw}=\bar{z}\bar{w}
\overline{z/w}=\bar{z}/\bar{w}
辐角和辐角主值
如果设定实数单位 1 作为水平正方向，虚数单位 \mathrm{i} 作为竖直正方向，得到的就是直角坐标视角下的复平面。
表示复数 z 的位置，也可以借助于极坐标 (r, \theta) 确定。前文已经提到了 r 为复数 z 的模。
从实轴正向到 非零 复数 z=x+\mathrm{i}y 对应向量的夹角 \theta 满足关系：
\tan \theta=\frac{y}{x}
称为复数 z 的 辐角 ，记为：
\theta= \operatorname{Arg} z
任一个 非零 复数 z 有无穷多个辐角，故 \operatorname{Arg} z 事实上是一个集合。借助开头小写的 \operatorname{arg} z 表示 其中一个特定值 ，满足条件：
-\pi<\operatorname{arg} z \le \pi
称 \operatorname{arg} z 为 辐角主值 或 主辐角 。辐角就是辐角主值基础上加若干整数个（可以为零或负整数）2k\pi ，即 \operatorname{Arg} z = \{\operatorname{arg} z + 2k\pi \mid k\in \mathbf Z\} 。
需要注意的是两个辐角主值相加后不一定还是辐角主值，而两个辐角相加一定还是合法的辐角。
称模小于 1 的复数，在复平面上构成的图形为 单位圆 。称模等于 1 的复数为 单位复数 ，全体单位复数在复平面上构成的图形为 单位圆周 。在不引起混淆的情况下，有时单位圆周也简称单位圆。
在极坐标的视角下，复数的乘除法变得很简单。复数乘法，模相乘，辐角相加。复数除法，模相除，辐角相减。
多项式的系数表示法和点值表示法
FFT 其实是一个用 O(n\log n) 的时间将一个用系数表示法 表示的多项式转换成它的点值表示法 的形式的算法。
多项式的点指表示法 和系数表示法 可以相互转换。
系数表示法
一个 n-1 次的 n 项多项式 f(x) 可以表示为：
f(x)=\sum^{n-1}_{i=0}{a_ix^i}
也可以直接用每一项的系数表示为：
f(x)=\left\{a_1, a_2,\cdots,a_{n-1}\right\}
点值表示法
从某种意义上来说，一个 n 次多项式可以看作是一个 n 次函数。如果看成是 n 次函数的话，那就意味着我们可以考虑代入 n 个不同的 x ，根据这 n 个不同的 x 计算出对应的 f(x) 的值。每一个值都可以看作是平面直角坐标系上的点 (x, f(x)) ，这几个点就可以确定一个唯一的多项式 f(x) 。即有且仅有一个多项式满足：
\forall k \in \R,f(x_k)=y_k
每一个点都可以确定一个方程，将这几个方程联立起来求解就可以得出最终的多项式。那么多项式 f(x) 还可以用以下方式来表示：
f(x)=\left\{(x_0, f(x_0)), (x_1, f(x_1)), (x_2, f(x_2)), \cdots, (x_{n-1}, f(x_{n-1}))\right\}
这就是多项式的点值表示法。
高精度乘法下两种表示法的区别
令两个多项式分别为 f(x)=\displaystyle\sum^{n-1}_{i=0}{a_ix^i} 和 g(x)=\displaystyle\sum^{n-1}_{i=0}{b_ix^i} 。
对于多项式的系数表示法，我们在相乘的时候需要每一项都逐一相乘，时间复杂度为 O(n^2) 。即：
f(x)\cdot g(x)=\prod^n_{i=1}\prod^n_{j=1}a^ib^j
对于多项式的点值表示法，我们在相乘的时候需要将每一个对应的点值相乘，时间复杂度为 O(n) 。即：
f(x)=\left\{(x_0, f(x_0)), (x_1, f(x_1)), (x_2, f(x_2)), \cdots, (x_{n-1}, f(x_{n-1}))\right\} \\
g(x)=\left\{(x_0, g(x_0)), (x_1, g(x_1)), (x_2, g(x_2)), \cdots, (x_{n-1}, g(x_{n-1}))\right\} \\
f(x)\cdot g(x)=\left\{(x_0 \times x_0, f(x_0)\times g(x_0), (x_1 \times x_1, f(x_1)\times g(x_1)), (x_2 \times x_2, f(x_2)\times g(x_2)), \cdots, (x_{n-1} \times x_{n-1}, f(x_{n-1})\times g(x_{n-1}))\right\} \\ \\
朴素 系数转点值的算法叫 DFT（离散傅里叶变换） ，点值转系数叫 IDFT（离散傅里叶逆变换） 。
虽然朴素 DFT 和 IDFT 的时间复杂度为 O(n^2) ，但我们也许可以做一些简单的优化。
DFT
注意：从这里开始所有的 n 都默认为 2 的整数次幂。
对于任意的多项式系数表示法转点值表示法，可以随便代入几个值，但是暴力计算需要 O(n^2) 的时间复杂度。
实际上，可以代入一组特定的 x 以减少时间复杂度。假设有这样一组 x ，每个 x 的若干次方 都等于 1 ，我们就无需做全部次方运算了。
考虑以坐标系原点为圆心，$1$ 为半径的单位圆。那么位于单位圆上的点的若干次方等于 $1$。将这个圆进行 $n$ 等分，这里取 $n=8$，则总共有 $8$ 个点。对这 $8$ 个点逐一编号，从$(1,0)$ 开始逆时针编号。
记编号为 $n$ 的点所代表的复数的值为 $\omega^k_n$，由**模长相乘，极角相加**可知 $(\omega^1_n)^k=w^k_n$，其中 $\omega^k_n$ 称之为 $n$ 次单位根，且每一个 $\omega$ 都可以很好地求出。
$$
\omega^k_n=\cos\frac{k}{n}2\pi+i\sin\frac{k}{n}2\pi
$$
$\omega^0_n,\omega^1_n,\cdots,\omega^{n-1}_n$ 即为我们要代入的 $x_0, x_1,\cdots,x_{n-1}$。
## 单位根的一些性质
$$
\omega^k_n=\omega^{2k}_{2n} \\
\omega^{k+\frac{n}{2}}_n=-\omega^k_n \\
\omega^0_n=\omega^n_n
$$
文字说明：
1. 它们表示的点（或向量）表示的复数**是相同的**。
2. 它们表示的点**关于原点对称**，所表示的复数**实部相反**，所表示的向量**等大反向**。
3. 都等于 $1$，或 $1+0i$。
证明如下：
$$
\omega^k_n=\cos\frac{k}{n}2\pi+i\sin\frac{k}{n}2\pi=\cos\frac{2k}{2n}2\pi+i\sin\frac{2k}{2n}2\pi=\omega^{2k}_{2n} \\
\omega^{\frac{n}{2}}_n=\cos\frac{\frac{n}{2}}{n}2\pi+i\sin\frac{\frac{n}{2}}{n}2\pi=\cos\pi+i\sin\pi=-1 \\
$$
## FFT
虽然 DFT 处理出了一组 $\omega$ 作为代入多项式的 $x$，但还需要代入单位根暴力计算……
但 DFT 可以运用分治来做。
首先设一个多项式：
$$
f(x)=\sum^{n-1}_{i=0}a_ix^i=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1}
$$
按 $f(x)$ 下标的**奇偶性**进行拆分：
$$
f(x)=\sum^{n-1}_{i=0}a_ix^i=a_0+a_1x+a_2x^2+\cdots+a_{n-1}x^{n-1} \\
=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+a_5x^5+\cdots+a_{n-1}x^{n-1}) \\
=(a_0+a_2x^2+a_4x^4+\cdots+a_{n-2}x^{n-2})+x(a_1+a_3x^2+a_5x^4+\cdots+a_{n-1}x^{n-2}) \\
$$
可以看出来两边好像有一点相似，于是再设两个多项式，则有：
$$
f_1(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1} \\
f_2(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1} \\
f(x)=f_1(x^2)+xf_2(x^2) \\
$$
此时再设 $k<\frac{n}{2}$，把 $\omega^k_n$ 作为 $x$ 代入多项式进行计算：
$$
f(\omega^k_n)=f_1((\omega^k_n)^2)+\omega^k_nf_2((\omega^k_n)^2) \\
=f_1(\omega^{2k}_n)+\omega^k_nf_2(\omega^{2k}_n) \\
=f_1(\omega^k_{\frac{n}{2}})+\omega^k_nf_2(\omega^k_{\frac{n}{2}}) \\
$$
对于 $f(\omega^k_{\frac{n}{2}})$，代入 $\omega^{k+\frac{n}{2}}_n$：
$$
f(\omega^{k+\frac{n}{2}}_n)=f_1(\omega^{k+\frac{n}{2}}_n)+\omega^{k+\frac{n}{2}}_nf_2((\omega^{k+\frac{n}{2}}_n)) \\
=f_1(\omega^{2k}_n\omega^n_n)-\omega^k_nf_2(\omega^{2k}_n\omega^n_n) \\
=f_1(\omega^{2k}_n)-\omega^k_nf_2(\omega^{2k}_n) \\
=f_1(\omega^k_{\frac{n}{2}})-\omega^k_nf_2(\omega^k_{\frac{n}{2}}) \\
$$
这个操作叫**蝴蝶变换**。
由上述变换过程可知，$f(\omega^k_n)$ 和 $f(\omega^{k+\frac{n}{2}}_n)$ 两个多项式经过变换之后只有最后一项的系数的符号不同。
也就是说，要求 $f(\omega^k_n)$ 和 $f_(\omega^{k+\frac{n}{2}}_n)$ 的值，只要求 $f_1(\omega^k_{\frac{n}{2}})$ 和 $f_2(\omega^k_{\frac{n}{2}})$ 的值。
每一次回溯时只扫当前前面一半的序列，即可得出后面一半序列的答案。注意到 $n=1$ 时只有一个常数项，此时直接返回答案。
总的时间复杂度为 $O(n\log n)$。
## IFFT
我们现在已经将两个多项式进行相乘（亦可称之为求卷积），做了两遍 FFT 也求解出了积的多项式的点值表示法形式的值。然而，我们平时常常用的是系数表示法的多项式，点值表示法的值对于求解答案来说意义不大。那么如何把点值表示的多项式快速转回系数表示法呢？
有这样一个结论：**一个多项式在分治的过程中乘上单位根的共轭复数，分治完的每一项除以 $n$ 即为原多项式的每一项系数。**
也就是说 FFT 可以和 IFFT 一起处理。
## 朴素 FFT
很明显，FFT 只能处理 $n$ 为 $2$ 的整数次幂的多项式，所以在 FFT 前一定要把 $n$ 调到 $2$ 的次幂。
递归**常数太大**，要考虑优化…
## 迭代 FFT
**每个位置分治后的最终位置为其二进制翻转后得到的位置**
这样的话我们可以先把原序列变换好，把每个数放在最终的位置上，再一步一步向上合并。
```cpp
#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#define int long long
#define MAXN 5000005
using namespace std;
const double pi = acos(-1);
int n, m, limit = 1, len;
int c[MAXN];
string s;
struct complex{
double real, imag;
complex(double x = 0, double y = 0){
real = x;imag = y;
}
}a[MAXN], b[MAXN];
complex operator +(complex a,complex b){return complex(a.real + b.real, a.imag + b.imag);}
complex operator -(complex a,complex b){return complex(a.real - b.real, a.imag - b.imag);}
complex operator *(complex a,complex b){return complex(a.real * b.real - a.imag * b.imag, a.real * b.imag + a.imag * b.real);}
void FFT(complex *a, int op){
for(int i = 0 ; i < limit ; i ++){
if(i < c[i])swap(a[i], a[c[i]]);
}
for(int mid = 1 ; mid < limit ; mid <<= 1){
complex t(cos(pi / mid), op * sin(pi / mid));
for(int r = mid << 1 , j = 0 ; j < limit ; j += r){
complex w(1, 0);
for(int l = 0 ; l < mid ; l++ , w = w * t){
complex x = a[j + l], y = w * a[j + mid + l];
a[j + l] = x + y ; a[j + mid + l] = x - y;
}
}
}
}
signed main(){
ios::sync_with_stdio(0);
cin >> s;n = s.size() - 1;
for(int i = 0 ; i <= n ; i ++)a[n - i].real = s[i] - '0';
cin >> s;m = s.size() - 1;
for(int i = 0 ; i <= m ; i ++)b[m - i].real = s[i] - '0';
while(limit <= n + m){
limit <<= 1;
len++;
}
for(int i = 0 ; i < limit ; i ++)c[i] = (c[i >> 1] >> 1) | ((i & 1) << (len - 1));
FFT(a, 1);FFT(b, 1);
for(int i = 0 ; i <= limit ; i ++)a[i] = a[i] * b[i];
FFT(a, -1);
memset(c, 0, sizeof(c));
for(int i = 0 ; i <= n + m + 1 ; i ++)c[i]=a[i].real / limit + 0.5;
for(int i = 0 ; i <= n + m ; i ++){
c[i + 1] += c[i] / 10;
c[i] %= 10;
}
limit = n + m + 1;
while(c[limit] != 0){
c[limit + 1] += c[limit] / 10;
c[limit] %= 10;
limit++;
}
for(int i = limit - 1 ; i >= 0 ; i --)cout << c[i];
cout << endl;return 0;
}
```