FFT 快速傅里叶变换

前置知识

复数

多项式

形如

\[F(x)=\sum_{i=0}^n a_i x^i \]

的式子。其中 \(n\) 为非负整数。我们只需知道所有 \(a_i\)，就可以确定一个多项式。这就是多项式的系数表示法。

为了方便运算，我们引入多项式的点值表示法。

对于一个 \(n\) 次多项式 \(F(x)\)，我们可以用 \(n+1\) 个互不相同的点 \(\{(x_0,F(x_0)),(x_1,F(x_1)),\cdots,(x_{n+1},F(x_{n+1}))\}\) 来确定这个多项式。我们下文称这个集合为 \(S_{F(x)}\)。

系数表示法转为点值表示法的过程叫做 DFT，反之叫 IDFT。

卷积

\(H(x)=F(x)\times G(x)\) 表示 \(H(x)\) 为 \(F(x)\) 和 \(G(x)\) 的卷积。

刚才的多项式乘法求得的式子并不是标准的多项式形式。我们设 \(H(x)=F(x)\times G(x)=\sum_{k=0}^{n+m}c_kx^k\)。则每项的系数 \(c_k=\sum_{i+j=k}a_ib_j\)，变为更容易计算的形式

\[c_k=\sum_{i=\max(0,k-m)}^{\min(n,k)}a_ib_{k-i} \]

若使用点值表示法，设

\[S_{F(x)}=\{(x_i,y_i)\} \]

\[S_{G(x)}=\{(x_i,y'_i)\} \]

则有

\[S_{H(x)}=\{(x_i,y_iy'_i)\} \]

单位根

根据代数基本定理^[1]，\(x^n=1\) 有 \(n\) 个根，这 \(n\) 个根都称为单位根。记作 \(\{\omega_n^k\mid k=0,1,\cdots,n-1\}\)，其中，\(\omega_n^0=1\)。在复平面上，它们刚好将单位圆 \(n\) 等分。一般说的单位根 \(\omega_n\)，指从 \((1,0)\) 开始逆时针方向上的第一个根。

一些下面会用到的小式子：

对于偶数次单位根，有 \(\omega_n^i=-\omega_n^{i+\frac{n}{2}}\)（其实就是在复平面上关于原点中心对称）。
\(\omega_{2n}^{2k}=\omega_n^k\)
\((\omega_n^k)^2=\omega_n^{2k}\)

FFT

显然，朴素算法求解 \(F(x) * G(x)\) 的时间复杂度为 \(O(nm)\)，而 FFT 可以让我们在 \(O(n\log n)\) 的时间复杂度内计算两个 \(n\) 次多项式的乘法。基本思想是分治。

我们只需对 \(F(x)\) 和 \(G(x)\) 进行 DFT，计算 \(S_{H(x)}=S_{F(x)}*S_{G(x)}\)，最后再对 \(S_{H(x)}\) IDFT 即得 \(H(x)\)。

接下来说 DFT 的过程。

对于 \(F(x)\)，将其划分为奇次与偶次两部分。

\[F(x)=\sum_{i=0}^{\frac{n}{2}}a_{2i}x^{2i}+\sum_{i=0}^{\frac{n}{2}}a_{2i+1}x^{2i+1} \]

将右半部分提出一个 \(x\)

\[F(x)=\sum_{i=0}^{\frac{n}{2}}a_{2i}x^{2i}+x\sum_{i=0}^{\frac{n}{2}}a_{2i+1}x^{2i} \]

将前后两部分用新的多项式表示

\[F_1(x)=\sum_{i=0}^{\frac{n}{2}}a_{2i}x^{i} \]

\[F_2(x)=\sum_{i=0}^{\frac{n}{2}}a_{2i+1}x^{i} \]

\[F(x)=F_1(x^2)+xF_2(x^2) \]

这时我们代入 \(\omega_n^k\)，可得

\[\begin{aligned} F(\omega_n^k)&=F_1((\omega_n^k)^2)+\omega_n^kF_2((\omega_n^k)^2)\\ &=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{aligned}\]

同理，代入 \(\omega_n^{k+\frac{n}{2}}=-\omega_n^k\)，得

\[F(\omega_n^{k+\frac{n}{2}})=F_1(\omega_\frac{n}{2}^k)-\omega_n^kF_2(\omega_\frac{n}{2}^k) \]

所以，我们可以根据 \(F_1(\omega_\frac{n}{2}^k)\) 和 \(F_2(\omega_\frac{n}{2}^k)\) 求出 \(F(\omega_n^k)\) 和 \(F(\omega_n^{k+\frac{n}{2}})\)。这种做法只能处理长度为 \(2\) 的正整次幂的多项式，所以我们要把高次系数补为 \(0\)。

接下来是 IDFT。它的操作与 DFT 极像，就是将 \(\omega_n^k\) 变为 \(\omega_n^{-k}\)，并在最后乘 \(\dfrac{1}{n}\)。

现实计算中，递归处理效率较低，我们使用位逆序置换优化和蝶形运算优化，直接将值排列为特定的顺序，避免了递归和额外的临时数组。

位逆序置换与蝶形运算

我们要想避免递归，就要将需要一起计算的部分放在一起。

我们以 \(7\) 次多项式为例，有 \(8\) 个 \(a_i\)，具体划分方式如图

规律很难注意到，我直接说了，就是将每个下标的二进制反转，以反转后的数为新下标。例如，\(3\) 的二进制是 \((011)_2\)，反转后为 \((110)_2\)，即 \(6\)，从图上来看，\(a_3\) 确实到了原 \(a_6\) 的位置。

位逆序置换后，我们可以直接计算 \(F_1(x)\) 和 \(F_2(x)\) 而无需临时数组，因为计算要用到的数与计算完成后的数应当被存在相同的下标内，直接覆盖原数就行了。具体地，计算 \(F(\omega_n^k)\) 和 \(F(\omega_n^{k+\frac{n}{2}})\) 时 \(F_1(\omega_\frac{n}{2}^k)\) 的值在下标 \(k\)，\(F_2(\omega_\frac{n}{2}^k)\) 的值在下标 \(k+\dfrac{n}{2}\)。

#include<iostream>
#include<cmath>
#include<vector>
using namespace std;
const double PI=4*atan(1);//这个 π 的精确度非常高
struct Complex{
    double real,imag;
    friend Complex operator+(const Complex &x,const Complex &y){
        return {x.real+y.real,x.imag+y.imag};
    }
    friend Complex operator-(const Complex &x,const Complex &y){
        return {x.real-y.real,x.imag-y.imag};
    }
    friend Complex operator*(const Complex &x,const Complex &y){
        return {x.real*y.real-x.imag*y.imag,x.real*y.imag+y.real*x.imag};
    }
};
void DFT(vector<Complex> &f,int n,int rev){
    for(int i=0,j=0;i<n;i++){
        if(i<j) swap(f[i],f[j]);
        for(int k=n>>1;(j^=k)<k;k>>=1);
    }
    for(int k=2;k<=n;k<<=1){
        double arg=2*PI*rev/k;
        Complex wn={cos(arg),sin(arg)};
        for(int i=0;i<n;i+=k){
            Complex w={1,0};
            for(int j=0;j<k/2;j++){
                Complex f1=f[i+j];
                Complex f2=f[i+j+k/2];
                f[i+j]=f1+w*f2;
                f[i+j+k/2]=f1-w*f2;
                w=w*wn;
            }
        }
    }
    if(rev==-1)
        for(int i=0;i<n;i++)
            f[i]={f[i].real/n,f[i].imag/n};
}
int n,m,len;
vector<Complex> a,b,c;
int main(){
    ios::sync_with_stdio(0),cin.tie(0),cout.tie(0);
    cin>>n>>m;
    for(int i=0,x;i<=n;i++){
        cin>>x;
        a.push_back({double(x),0});
    }
    for(int i=0,x;i<=m;i++){
        cin>>x;
        b.push_back({double(x),0});
    }
    len=1;while(len<=n+m) len<<=1;
    a.resize(len),b.resize(len),c.resize(len);
    DFT(a,len,1);
    DFT(b,len,1);
    for(int i=0;i<len;i++)
        c[i]=a[i]*b[i];
    DFT(c,len,-1);
    for(int i=0;i<=n+m;i++)
        cout<<int(round(c[i].real))<<' ';
    return 0;
}

代码实现的一些注意事项：