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

作者丨星夜@知乎（已授权）
来源丨https://zhuanlan.zhihu.com/p/347091298
编辑丨极市平台

前言

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

FFT用于用来加速多项式乘法，可以以 $O(nlogn)$ 的时间复杂度解决朴素高精度 $O(n^2)$ 复杂度能够解决的多项式乘法。

复数

星夜：虚数与复数与欧拉公式zhuanlan.zhihu.com

系数与点值表示法

对于任意的 $n$ 项 $n-1$ 次整系数多项式 $f(x)$ ，显而易见我们可以将其表示为 $a_0x^{n-1}+a1x^{n-2}+...+a{n-1}x^0$ 的形式。

我们同样可以用每一项的系数表示 $f(x)$ ，即：

$$
f(x)=\left{ a_0,a_1,a2,...,a{n-1} \right}\quad\
$$

这就是系数表示法。

当我们将多项式 $f(x)$ 看作一个函数时，如果我们将 $n$ 个不同的数 $x_0,x1,...,x{n-1}$ 带入，就能得到 $n$ 个不同的点。

这 $n$ 个不同的点同样可以看作是 $n$ 个不同的一元 $n-1$ 次方程联立，这 $n$ 个点或者方程也能反过来唯一确定 $f(x)$ 。

所以$f(x)$也可以表示为

$$
f(x)=\left{ (x_0,f(x_0)),(x_1,f(x1)),...,(x{n-1},f(x_{n-1})) \right}\quad\
$$

这就是点值表示法。

我们假设两个多项式 $f(x),g(x)$ 分别为：

$$
f(x)=\left{ (x_0,f(x_0)),(x_1,f(x1)),...,(x{n-1},f(x_{n-1})) \right}\quad\
$$

$$
g(x)=\left{ (x_0,g(x_0)),(x_1,g(x1)),...,(x{n-1},g(x_{n-1})) \right}\quad \
$$

其乘积为：

$$
h(x)=\left{ (x_0,f(x_0)\cdot g(x_0)),(x_1,f(x_1)\cdot g(x1)),...,(x{n-1},f(x{n-1})\cdot g(x{n-1})) \right}\quad \
$$

显然点值表示法完成多项式乘法的时间复杂度仅有 $O(n)$ ，所以剩下的问题就是如何将系数表示法转为点值表示法了。

单位根

以单位圆点为起点，单位圆的 $n$ 等分点为终点，在单位圆上可以得到 $n$ 个复数，设幅角为正且最小的复数为 $\omega_n$ ，称为 $n $ 次单位根。

由于复数乘法模数相乘，辐角相加的性质，编号为 $k$ 的点 $\omega_n^k=(\omega_n)^k$ 。

由欧拉公式可得

$$
\omega_{n}^{k}=\cos\ k \cdot\frac{2\pi}{n}+i\sin k\cdot\frac{2\pi}{n}\quad\
$$

显然 $\omega_n^0=\omega_n^n$ 。

单位根同样具有以下性质：

$$
\omega_{rn}^{rk}=\cos rk\frac{r\pi}{rn}+i\sin rk\frac{r\pi}{2n}=\omega_n^k \quad r \in \mathbb{Z}^+\
$$

$$
\begin{align}\omega _{n}^{k+\frac{n}{2}}&=\omega_n^k \cdot (\cos\frac{n}{2}\frac{2\pi}{n}+i\sin\frac{n}{2}\frac{2\pi}{n} )\&=\omegan^k \cdot (\cos \pi +i\sin \pi)\&=-\omega {n}^{k}\end{align}\quad\
$$

上述性质从角的角度看就很容易理解。

快速傅里叶变换

我们假设多项式 $f(x)=a_0+a1*x+ \dots+a{n-2}x^{n-2}+a_{n-1}x^{n-1}$ ，那么我们可以按 $a$ 下标的奇偶性将 $f(x)$ 分为两半。

也就是 $\begin{align}f(x)&=(a_0+a_2{x^2}+a4{x^4}+\dots+a{n-2}x^{n-2})\&+(a_1x+a_3{x^3}+a5{x^5}+ \dots+a{n-1}x^{n-1})\end{align}\quad\ $

我们令 $f_1(x)=a_0+a_2{x^1}+a4{x^2}+\dots+a{n-2}x^{\frac{n}{2}-1}\quad\$ $f_2(x)=x(a_1+a_3{x}+a5{x^2}+ \dots+a{n-1}x^{\frac{n}{2}-1})\quad\$

显然 $f(x)=f_1(x^2)+xf_2(x^2)$ 。

带入 $\omega_n^k(k<\frac{n}{2}) $ 可得 $\begin{align}f(\omega_n^k)&=f_1(\omega_n^{2k})+\omega_n^kf_2(\omega_n^{2k})\&=f1(\omega{\frac{n}{2}}^{k})+\omega_n^kf2(\omega{\frac{n}{2}}^{k}) \end{align}\quad\$

带入 $\omega_n^{k+\frac{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}\omega_n^n)-\omega_n^kf_2(\omega_n^{2k}\omega_n^n)\&=f_1(\omega_n^{2k})-\omega_n^kf_2(\omega_n^{2k}) \end{align}\quad\$

我们发现上面的两个式子只有一个常数项不同，所以当我们求出第一个式子的值后可以以 $O(1)$ 的时间复杂度求出第二个式子的值。

所以我们就成功将原来的问题缩小到了之前的一半，而缩小后的问题还能继续缩小。

所以显然FFT的时间复杂度即为 $O(nlogn)$ 。

快速傅里叶逆变换

我们把两个多项式相乘（或者说求卷积），我们跑完FFT并相乘后我们就得到了积的多项式的点值表示，现在我们需要将点值表示转为系数表示，这个转换的过程我们称为离散傅里叶逆变换（IDFT）。

在某种意义上，我们可以将DFT的操作视作以下操作：

$$
\begin{bmatrix}1 & 1 & 1 & \cdots & 1 \1 & x_1 & x_1^2 & \cdots & x_1 ^ n\1 & x_2 & x_2^2 & \cdots & x2 ^ n \ \vdots & \vdots &\vdots & \ddots & \vdots \ 1 & x{n-1} & x{n-1}^2 & \cdots & x{n-1} ^ n\end{bmatrix} \cdot \begin{bmatrix}a0 \ a 1 \ a 2 \ \vdots \ a n\ \end{bmatrix} = \begin{bmatrix}y0 \ y 1 \ y 2 \ \vdots \ y n\ \end{bmatrix}\quad\
$$

其中 $A=\begin{bmatrix}a0 \ a 1 \ a 2 \ \vdots \ a n\ \end{bmatrix}\quad\$ 意义为多项式的系数表示法，$B=\begin{bmatrix}y0 \ y 1 \ y 2 \ \vdots \ y n\ \end{bmatrix}\quad\$ 意义为多项式的点值表示法。

我们假设有 $C=\begin{bmatrix}c0 \ c 1 \ c 2 \ \vdots \ c n\ \end{bmatrix}\quad\$ 满足 $ck=\sum{i=0}^{n-1}y_i(\omega_n^{-k})^i$ 。

也就是说， $C$ 为 $B$ 代表一个多项式的系数表示法时这个多项式的点值表示法。

我们可以推出 $\begin{align}ck&=\sum{i=0}^{n-1}y_i(\omegan^{-k})^i\&=\sum{i=0}^{n-1}(\sum_{j=0}^{n-1}a_j(\omega_n^j)^i(\omegan^{-k})^i)\&=\sum{i=0}^{n-1}\sum_{j=0}^{n-1}a_j(\omegan^{j-k})^i \&=\sum{j=0}^{n-1}aj(\sum{i=0}^{n-1}(\omega_n^{j-k})^i) \end{align}\quad\$

我们设 $g(x)=\sum_{i=0}^{n-1}x^i$ ，那么有 $\omega_n^kS(\omega_n^k)-S(\omega_n^k)=(\omega_n^k)^{n}-1\quad\$

变形得 $\begin{align}g(\omega_n^k)&=\frac{(\omega_n^k)^{n}-1}{\omega_n^k-1}\&=\frac{(\omega_n^n)^{k}-1}{\omega_n^k-1}\&=\frac{1-1}{\omega_n^k-1}\end{align}\quad\$

也就是 $g(\omega_n^k)=0(k\ne 0)\quad\$

而当 $k=0$ 时，显然 $g(\omega_n^0)=n\quad\$

所以由于 $ck=\sum{j=0}^{n-1}aj(\sum{i=0}^{n-1}(\omega_n^{j-k})^i)\quad\$ 有$c_k=na_k\quad\$

也就是 $a_k=\frac{c_k}{n}\quad\$

所以一个多项式在分治的过程中乘上单位根的共轭复数，分治完的每一项除以 $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>>1;
    complex<double> A1[mid+1],A2[mid+1];
    for(int i=0;i<=n;i+=2){//拆分多项式
        A1[i>>1]=a[i],
        A2[i>>1]=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>>1)]=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的二的正整数次方
    FFT(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(int i=0;i<=n+m;++i)//输出
        printf("%.0f ",a[i].real()/len+eps);//记得要除n，eps为解决掉精度问题
    return 0;
}