快速傅里叶变换（FFT）超详解

前言

快速傅里叶变换 (Fast Fourier Transform)，即利用计算机计算离散傅里叶变换（DFT）的高效、快速计算方法的统称，简称FFT，于1965年由J.W.库利和T.W.图基提出。

对多项式 f(x)=\sum_{i=0}^{n}a_ix^i,g(x)=\sum_{i=0}^{n}b_ix^i ，定义其乘积 fg 为 (fg)(x)=\left(\sum_{i=0}^{n}a_ix^i\right)\left(\sum_{i=0}^{n}b_ix^i\right)\\ 显然我们可以以 \mathrm O(n^2) 的复杂度计算这个乘积的每一项的系数。

但FFT可以以 \mathrm O(n\log n) 的时间复杂度来计算这个乘积。

复数多项式的系数与点值表示

对于任意 n+1 次多项式 f(x)=\sum_{i=0}^{n}a_ix^i\in X_n ，只需知道每一项的性价比手机系数就可以唯一确定 f(x) ，故 f(x) 可以被写作 ( a_0,a_1,a_2,...,a_{n} ) 的形式，即其系数的有序排列，这就是多项式的系数表示。

我们还有以下结论：

n 次多项式上 n+1 个不同的点能唯一确定这个多项式

即将 n+1 个不同的数 x_0,x_1,...,x_{n} 带入 f(x) ，就能得到 n+1 个不同的数值 y_i=f(x_i) ，进而存在唯一的满足\forall 0\leq i\leq n,f(x_i)=y_i\\ 的 n 次多项式 f 。上述结论可以通过因式定理与代数基本定理，或者用范德蒙特行列式的性质证明。

从而 f(x) 也可以被写作

\left\{ (x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n},f(x_{n})) \right\}\\这就是多项式的点值表示，在 x_0,x_1,...,x_{n} 确定时，令 \tau(f) 表示 f 的点值表示。

我们假设两个多项式 f(x),g(x) 的点值表示法分别为

\tau(f)=\left\{ (x_0,f(x_0)),(x_1,f(x_1)),...,(x_{n},f(x_{n})) \right\}\\ \tau(g)=\left\woot{ (x_0,g(x_0)),(x_1,g(x_1)),...,(x_{m},g(x_{m})) \right\} 则

\tau(fg)=\left\{ (x_0,f(x_0)\cdot g(x_0)),(x_1,f(x_1)\cdot g(x_1)),\cdots,(x_{n+大学数学题m},f(x_{n+m}))\right\}\quad \\ 于是多项式的乘法在点值表示法下可以以 \mathrm O(n) 的复杂度计算，所以我们如果能够在较低的时间复杂度内将系数表示法转化为点值表示法，再将点值表示法转回系数表示法，就能以较低的时间复杂度计算多项式的乘法。

单位根

以单位圆点为起点，单位圆的 n 等分点为终点，在单位圆上可以得到 n 个复数，设幅角为正且最小的复数为 \omega_n ，称为 n 次单位根，即 \omega_{n}=\cos \frac{2\pi}{n}+{\rm i}\sin \frac{2\pi}{n}\\ 由欧拉公式

\omega_{n}^{k}=\cos \frac{2k\pi}{n}+{\rm i}\sin \frac{2k\pi}{n}\\ 特别地， \omega_n^0=\omega_n^n=1 。

可以发现单位根具有以下性质：

\begin{align}\omega_{rn}^{rk}&=\cos \frac{2rk\pi}{rn}+{\rm i}\sin \frac{2rk\依然爱丽丝pi}{rn}=\omega_n^k &r \in \mathbb{N}^+\\ \omega _{n}^{k+\frac{n}{2}}&=\omega_n^k \cdot \left(\cos\left(\frac{n}{2}\cdot\frac{2\pi}{n}\right)+{\rm i}\sin\left(\frac{n}{2}\cdot \frac{2\pi}{n} \right)\right)&k\in\mathbb N^+\\ &=\omega_n^k \cdot (\cos \pi +{\rm i}\sin \pi)\\&=-\omega _{n}^{k}\\ \color{red}{\overline{\omega_n^k}}&\color{red}{=\cos\frac{2k\pi}n-{\rm i}\sin\frac{2k\pi}n}\\ &\color{red}{=\cos\left(2\pi-\frac{2k\pi}n\right)+{\rm i}\sin\left(2\pi-\frac{2k\pi虚拟信用卡}n\right)}\\ &\color{red}{=\omega_n^{n-k}}\end{align}\\

快速傅里叶变换

对多项式 f(x)=\sum_{i=0}^{n-1}a_ix^i\in X_{n-1} ，发际线调整不失一般性的，设 n=2^s\quad s\in\mathbb N （由于在多项式的乘法中我们可以将一个多项式等价地看作是次数更高的高次项系数均为零的多项式，故可以将 n 看作第一个大于等于它的 2 的整数次幂），考虑按 a_i 下标的奇偶性将 f(x) 中的项分为两部分，即 \begin{align}f(x)&=(a_0+a_2{x^2}+a_4{x^4}+\dots+a_{n-2}x^{n-2})\\&+(a_1x+a_3{x^3}+a_5{x^5}+ \dots+a_{n-1}x^{n-1})\end{align}\\ 令 \begin{align} f_1(x)&=a_0+a_2{x^1}+a_4{x^2}+\dots+a_{n-2}x^{\frac{n}{2}-1}\\ f_2(x)&=a_1+a_3{x}+a_5{x^2}+ \dots+a_{n-1}x^{\frac{n}{2}-1} \end{align}\\

则 f(x)=f_1(x^2)+xf_2(x^2)\\ 带入 x=\omega_n^k\quad(k<\frac{n}{2}) 可得 \begin{align}f(\omega_n^k)&=f_1(\omega_n^{2k})+\omega_n^kf_2(\omega_n^{2k})\\ &=f_1(\omega_{\frac{n}{2}}^{k})+\omega_n^kf_2(\omega_{\frac{n}{2}}^{k}) \end{align}\tag1 带入 \omega_n^{k+\f巴克利奥尼尔rac{n}{2}}(k<\frac{n}{2}) 可得 \begin{align}f(\omega_n^{k+\frac{n}{2}})&=f_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}}f_2(\omega_n^{2k+n})\\ &=f_1(\omega_n^{2k}\cdot\omega_n^n)-\omega_n^kf_2(\omega_n^{2k}\cdot\omega_n^n)\\&=f_1(\omega_n^{2k})-\omega_n^kf_2(\omega_n^{2k}) \\ &=f_1(\omega_{\frac{n}{2}}^{k})-\omega_n^kf_2(\omega_{\frac{n}{2}}^{k}) \end{align}\tag2 观察 (1),(2) 两式的结构，我们只需要求出 f_1(\omega_{\frac{n}{2}}^{k}),f_2(\omega_{\frac{n}{2}}^{k}) 即可 \mathrm O(1) 求出 (1),(2) 两式的值，进而再经过类似的步骤，我们可以以 \mathrm O(1) 的时间复杂度将问题继续转化为求 f_1(\omega_{\frac{n}{4}}^{k}),f_2(\omega_{\frac{n}{4}}^{k}) ……最终问题被转化为求 f_1(\omega_{1}^{k})=f_2(\omega_{1}^{k})=1\\ 故以 \mathrm O(\log n) 的时间复杂度可以求出 f(\omega_n^k),f(\omega_n^{k+\frac n2}) ，进而可以以 \mathrm O(n\log n) 的时间复杂度求出所有的 f(\omega_n^k) ，即我们以 \mathrm O(n\log n) 的时间复杂度将 f 的系数表示法转化为了点值表示法。

快速傅里叶逆变换

跑完FFT后我们就得到了多项式乘积的点值表示，现在我们需要将点值表示转回系数表示，这个转换的过程被称为离散傅里叶逆变泰国人妖电影换（IDFT）。

如果我们用矩阵将DFT的过程封装，那么DFT就相当于求

\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & (\omega_n^1)^1 &水晶面膜amp; (\omega_n^1)^2 & \cdots & (\omega_n^1) ^ {n-1}\\1 & (\omega_n^2)^1 & (\omega_n^2)^2 & \cdots & (\omega_n^2) ^ {n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ 1 & (\omega_n^{n-1})^1 & (\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1}) ^ {n-1}\end{bmatrix} \begin{bmatrix}a_0 \\ a _ 1 \\ a _ 2 \\ \vdots \\ a _ {n-1}\\ \end{bmatrix} = \begin{bmatrix}y_0 \\ y _ 1 \\ y _ 2 \\ \vdots \\ y _ {n-1}\\ \end{bmatrix}\quad\\

其中 \begin{bmatrix}a_0 \\ a _ 1 \\ a _ 2 \\ \vdots \\ a _ {n-1}\\ \end{bmatrix}\\ 意义为多项式的系数表示法，\begin{bmatrix}y_0 \\ y _ 1 \\ y _ 2 \\ \vdots \\ y _ {n-1}\\ \end{bmatrix}\\ 意义为多项式的点值表示法， y_i 即为 f(\omega_n^i) 。

注意到 x^n-1=(x-1)(x^{n-1seku}+x^{n-2}+\cdots+1)\\ 于是在 x=1 时 x^{n-1}+x^{n-2}+\cdots+1=n\\ 否则在 x^n-1=0 时 x^{n-1}+x^{n-2}+\cdots+1=0\\ 其中 x 为 n 次单位根。

受 n 次单位根这一特性的启发，注意到 x^{n-1}+x^{n-2}+\cdots+1 在其形式上与矩阵乘法的关系，可以发现 \begin{ali斯巴鲁森林人油耗gn}\begin{bmatrix}1&(\omega_n^{n-i})^1&({\omega_n^{n-i}})^2&\cdots&(\omega_{n}^{n-i})^{n-1}\end{bmatrix}\begin{bmatrix}1\\\ (\omega_n^{j})^1\\(\o垄断mega_n^{j})^2\\\vdots\\(\omega_n^{j})^{n-1}\end{bmatrix}&=\sum_{k=0}^{n-1}(\omega_n^{n-i})^k(\omega_n^j)^k\\ &=\sum_{k=0}^{n-1}(\omega_n^{n-i+j})^k\\ &=\begin{cases}n&i=j\\0&i\ne j\end{cases}\end{align}\\ 这表明 \mathbf C=\frac1n\begin{bmatrix}1&1&1&\cdots&1\\ 1&(\omega_n^{n-1})^1&({\omega_n^{n-1}})^2&\cdots&(\omega_{n}^{n-1})^{n-1}\\ 1&(\omega_n^{n-2})^1&({\omega_n^{n-2}})^2&\cdots&(\omega_{n}^{n-2})^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&(\o无限团mega_n^1)^1&{(\omega_n^1)}^2&\cdots&(\omega_{n}^1)^{n-1}\end{bmatrix}\\ 为 \begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & (\omega_n^1)^1 & (\omega_n^1)^2 & \cdots & (\omega_n^1) ^ {n-1}\\1 & (\omega_n^2)^1 &livall (\omega_n^2)^2 & \cdots & (\omega_n^2) ^ {n-1} \\ \vdots萌芽论坛 & \vdots &\vdots & \ddots & \vdots \\ 1 & (\omega_n^{n-1})^1 & (\omega_n^{n-1})^2 & \cdots & (\omega_n^{n-1}) ^ {n-1}\end{bmatrix}\\ 的逆矩阵。

由上文标红部分 \begin{align} \mathbf C&=\frac1n\begin{bmatrix}1&1&1&\cdots&1\\ 1&(\omega_n^{n-1})^1&({\omega_n^{n-1}})^2&\cdots&(\omega_{n}^{n-1})^{n-1}\\ 1&(\omega_n^{n-2})^1&({\omega_n^{n-2}})^2&\cdots&(\omega_{n}^{n-2})^{n-1}\\ \vdots&\vdots&\vdots&\ddots&\vdots\\ 1&(\omega_n^1)^1&{(\omega_n^1)}^2&\cdots&(\omega_{n}^1)^{n-1}\end{bmatrix}\\ &=\frac1n\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \\1 & (\overline{\omega_n^1})^1 & (\overline{\omega_n^1})^2 & \cdots & (\overline{\omega_n^1}) ^ {n-1}\\1 & (\overline{\omega_n^2})^1 & (\overline{\omega_n^2})^2 & \cdots & (\overline{\omega_n^2}) ^ {n-1} \\ \vdots & \vdots &\vdots & \ddots & \vdots \\ 1 & \overline{(\omega_n^{n-1}})^1 & (\overline{\omega_n^{n-1}})^2 & \cdots & (\overline{\omega_n^{n-1}}) ^ {n-1}\end{bmatrix}\\\end{align}\\

所以将一个多项式在分治的过程中乘上的单位根变为其共轭复数，分治完的每一项除以 n 即可得到原多项式的每一项系数。

具体实现（递归版）//LuoguP3803#include <bits/stdc++.h>const int NR = 1 << 22;const double eps = 1e-6, pi = acos(-1.0);using namespace std;complex<double> a[NR], b[NR]; //complex为C++自带虚数int n, m;void FFT(complex<double> *a, int n, int inv) //inv为虚部符号，inv为1时FFT，inv为-1时IFFT{ if (n == 1) //如果需要转换的只有一项就直接返回 return; int mid = n / 2; complex<double> A1[mid + 1], A2[mid + 1]; for (int i = 0; i <= n; i += 2) //拆分多项式 { A1[i / 2] = a[i]; A2[i / 2] = a[i +肛检 1]; } FFT(A1, mid, inv); //递归分治 FFT(A2, 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] = A1[i] + w0 * A2[i]; a[i + n / 2] = A1[i] - w0 * A2[i]; }}int main(){ scanf("%d%d", &n, &m); for (int i = 0; i <= n; ++i) //输入第一个多项式 { double x; scanf("%lf", &x); a[i].real(x); //complex类型变量.real(x)意味着将实数部赋为x，real()返回实数部值 } for (int i = 0; i <= m; ++i) //输入第二个多项式 { double x; scanf("%lf", 考研英语复习资料&x); b[i].姐妹节real(x); } int len = 1 << max((int)ceil(log2(n + m)), 1); //由于FFT需要项数为2的整数次方倍，所以多项式次数len为第一个大于n+m的二的正整数次方 F电线电缆规格FT(a, len, 1); //系数表达转点值表达 FFT(b, len, 1); for (int i = 0; i <= len; ++i) a[i] = a[i] * b[i]; //O(n)乘法 FFT(a, len, -1); //点值表达转系数表达 for (i四大丑女nt i = 0; i <= n + m; ++i) //输出 printf("%.0f ", a[i].real() / len + eps); //记得要除n，eps为解决掉精度问题 return 0;}FFT的优化快速数论变换（NTT）

关于某没有精度损失还可以取模的比FFT快的小常数多项式乘法.jpg

例题