[数据结构与算法]快速傅里叶变换FFT简明教程

问题

给出长度为$n$的多项式$f(x)={\sum }_{i=0}^{n?1}{a}_i{x}^i$和长度为$m$多项式$g(x)={\sum }_{i=0}^{m?1}{b}_i{x}^i$。求解多项式$h(x)=f(x)\times g(x)$。
$(1\le n\le 10^5)$

系数表示法和点值表示法

多项式的两种表示方法，系数表示法和点值表示法。
系数表示法：我们存储$f(x)$的每一项$x^i$的系数$a_i$。
点值表示法：由大家都知道知，$N$个点可以确定一个最高次为$N?1$多项式，所以我们可以存储多项式上任意$N$个点，来确定这个多项式。

点值表示法有个性质是，对于多项式$f$的一个点值$(x,y_1)$，和多项式$g$的一个点值$(x,y_2)$，可以得到两个多项式相乘的多项式$h$在横坐标为$x$处的点为$(x,y_1\times y_2)$。

所以在对于长度为$n$的多项式$f$和长度为$m$的多项式$g$做乘法时，我们可以算出$f$的$n+m?1$个点值和$g$的$n+m?1$个点值，且要求两个多项式点值所取的$x$相同。
我们就可以通过遍历算出这$n+m?1$个横坐标的点值纵坐标相乘，来得到相乘后多项式的$n+m?1$个点值，也就是新多项式的点值表示法。
(注意：虽然我们存储$f$只需要$n$个点值，但是为了计算出长度为$n+m?1$的$h$，我们仍要先得到$f$的(n+m-1)个点值)

因此，假如我们可以通过$o(nlogn)$的办法，将系数表示法转化为点值表示法，然后将两个多项式点值$o(n)$相乘，最后通过$o(nlogn)$的办法将点值表示法转化回系数表示法，就可以解决这个问题了。

后面来介绍如何将系数表示发和点值表示法进行转换。

离散傅里叶变换（DFT）

对于没有周期的函数，我们可以认为其周期为$\infty$，用无穷个周期函数去拟合。
同理，对于一个长度为$N$的离散多项式$x(n)(n=0,1,2..N?1)$，可以通过$N$个周期函数来拟合。（$N$个周期函数只能拟合$N$个点，不能用来拟合连续多项式函数）
由此我们得到离散傅里叶级数展开式如下。
$$X(k)=\sum _{n=0}^{N?1}x(n)e^{?2\pi ik\frac{n}{N}}$$
其中$e^{?2\pi ik\frac{n}{N}}=cos(?2\pi k\frac{n}{N})+isin(?2\pi k\frac{n}{N})$是关于$n$的周期为$N$的函数。（由欧拉公式$e^{ix}=cosx+isinx$得）

离散傅里叶逆变换（IDFT）

对于变换后得到的$X(n)$，我们可以通过离散傅里叶逆变换公式得到原函数$x(n)$。与$DFT$公式差不多，在原来的基础上$i$变成了$?i$，前面乘上一个$\frac{1}{N}$。公式如下。
$$x(n)=\frac{1}{N}\sum _{k=0}^{N?1}X(k)e^{2\pi ik\frac{n}{N}}$$

快速傅里叶变换（FFT）

由前面讲过的，我们要解决的核心问题就是，算出多项式的$N$个点值。
大致做法是，我们求出函数在离散傅里叶变换后的N个点值，由点值乘出新多项式后，再将其逆变换回系数表示法。

我们先来观察式子$DFT$的式子。
我们为了简化式子，我们把$e^{\frac{2\pi ik}{N}}$用$W_N$代替。($W_N$是一个关于$k$的函数)
得到式子如下。
$$X(k)=\sum _{n=0}^{N?1}x(n)W_N^n(k)$$
然后观察$W_N$如下。
$$W_N(k)=e^{\frac{2\pi i}{N}}=cos(\frac{2\pi k}{N})+isin(\frac{2\pi k}{N})$$
发现$W_N$是一个周期为$N$的函数。

由此，我们发现变换后的函数$X(k)$是由$N$个周期为$N$的函数的幂次，乘上一个系数，相加而成。且这$N$个函数的实部和虚部都是三角函数（$cos$和$sin$）。

我们很轻易地发现了$W_N(k)=W_N(k+N)$，但是这对我们来说不重要。更重要的性质源于三角函数的另一个性质——对称性。
我们把$[0,2\pi ]$分成两部分$[0,\pi )$，和$(\pi ,2\pi ]$，可以发现我们算出左边 $t$ 处的值，就可以得到右边$\pi +t$处的值。
说白了就是$sin(t)=?sin(\pi +t)$和$cos(t)=?cos(\pi +t)$。

对应到上面就得到了一个结论，我们只要算出了$W_N(t)$，就可以得到$W_N(\frac{N}{2}+t)$。即$W_N(t)=?W_N(\frac{N}{2}+t)$。

上面只是简单的性质推导，正式的傅里叶变换在下面。

总而言之就是要分治，这个式子可以和分治有机结合。（我也不知道他怎么想出来的，公式推导在《数字信号处理》第四版的P122）
我们把枚举的奇数项和偶数项分离开。
$$\begin{gathered} X(k)=\sum _{n为偶数}x(n)W_N^n(k)+\sum _{n为奇数}x(n)W_N^n(k) \\ =\sum _{r=0}^{\frac{N}{2}?1}x(2r)W_N^{2r}(k)+\sum _{r=0}^{\frac{N}{2}?1}x(2r+1)W_N^{2r+1}(k) \\ =\sum _{r=0}^{\frac{N}{2}?1}{x}_1(r)W_N^{2r}(k)+W_N(k)\sum _{r=0}^{\frac{N}{2}?1}{x}_2(r)W_N^{2r}(k) \\ =\sum _{r=0}^{\frac{N}{2}?1}{x}_1(r)W_{\frac{N}{2}}^r(k)+W_N(k)\sum _{r=0}^{\frac{N}{2}?1}{x}_2(r)W_{\frac{N}{2}}^r(k) \\ =X_1(k)+W_N(k)X_2(k) \end{gathered}$$
其中$X_1$是$x_1$的傅里叶变换，$X_2$是$x_2$的傅里叶变换。
其中有一步$W_N^{2r}(k)=W_{\frac{N}{2}}^r$前面忘记证了，想必读者应该一下子就能证出来了。

到此位置，所有准备工作都做完了，容易发现只要分治做好两个小多项式的点值，大的自然就出来了。所以可以分治做。

最后式子如下。
$$\begin{gathered} X(k)=X_1(k)+W_N^k{X}_2(k),k=0,1,\dots ,\frac{N}{2}?1 \\ X(k+\frac{N}{2})=X_1(k)?W_N^k{X}_2(k),k=0,1,\dots ,\frac{N}{2}?1 \end{gathered}$$

但是有一个细节问题，计算$W_N(k)$的复杂度是$o(logn)$。而在分治的每一层我们需要求$N$次$W$函数，总共分治$logn$层，最终导致我们的复杂度变成了$o(nlognlogn)$。
那么解决办法是前面提的$W_N(t)=?W_N(\frac{N}{2}+t)$，我们每层只需要计算一半的$W$，也就是$\frac{N}{2}$个$W$，于是在计算$W$上的总开销就变成了$o(\frac{n}{2}lognlogn)$。有锅巴用。

然后又发现一个可以优化的地方。既然$W_N(k)=?W_N(\frac{N}{2}+k)$来减少$W$的计算次数（主要是减少$sin$和$cos$的计算次数），主要是来源于三角函数知道$[0,\pi ]$可以推演出$[\pi ,2\pi ]$。这样只需要在每层算$\frac{N}{2}$次$cos$和$sin$，所以我们再往下细化一层，使得每层只需要计算$\frac{N}{4}$次$cos$和$sin$。
那么同理，三角函数知道$[0,\frac{\pi }{2}]$就可以推算出$[\frac{\pi }{2},2\pi ]$。
所以，我们开始推$W_N(k+\frac{N}{4})$。
$$W_N(k+\frac{N}{4})=cos(\frac{2\pi k}{N}+\frac{\pi }{2})+isin(\frac{2\pi k}{N}+\frac{\pi }{2})=?sin(\frac{2\pi k}{N})+icos(\frac{2\pi k}{N})$$
然后，我们发现，我们对于某个$k$计算出的$sin$和$cos$函数，可以沿用至$k+\frac{N}{4}$，$k+\frac{N}{2}$，$k+\frac{3N}{4}$。
所以，在每一层，我们只需要计算$\frac{N}{4}$次$W$。
然后复杂度被优化到了$o(\frac{n}{4}lognlogn)$。仍然没有任何卵用。

其实我们只需要算$W_N(1)$，有$W_N(k+1)=W_N(k)\times W_N(1)$，这才是真的有用。证明可以考虑这个复数的几何意义，就是把位于圆心的单位圆划等分成了$N$个部分，取出第一个部分就相当于是一个旋转矩阵，幂次就是旋转的次数。
于是，我们只要在开始用三角函数算出$W_N(1)$，后面一直乘上它就行了。

然后就做完了。

关于转换回系数，也就是逆变换，发现式子和这个差不多，所以做法也基本一样。记得最后要除以$N$。这里就不重复介绍了。

为了方便起见，一般先把多项式长度拉到2的幂次，然后再进行变换，方便递归到底层为1。反过来说，要使用FFT板子时，注意长度应为2的幂次。
虽然是分治，但为了常数考虑，一般都用非递归版的写法。

递归版
递归版常数很大，强烈建议找个常数小的非递归版代码
#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod=1e9+7;
const double pi = acos(-1);
const int MAXN = 3000005;
struct cp
{
    double s,v;
    cp(double ss=0,double vv=0){
        s=ss;v=vv;
    }
};
cp operator * (cp a,cp b)
{
    return cp(a.s*b.s-a.v*b.v,a.s*b.v+a.v*b.s);
}
cp operator + (cp a,cp b)
{
    return cp(a.s+b.s,a.v+b.v);
}
cp operator - (cp a,cp b)
{
    return cp(a.s-b.s,a.v-b.v);
}
cp operator / (cp a,int b)
{
    return cp(a.s/b,a.v/b);
}

//预处理出2的幂次的cos和sin
double coss[MAXN],sinn[MAXN];

void fft(cp *a,int len,int opt)
{
    if(len==1)return ;
    static cp b[MAXN];
    int mid=len>>1,mid_2=mid>>1;
    for(int i=0;i<len;i++)b[(i>>1)|(i&1?mid:0)]=a[i];
    for(int i=0;i<len;i++)a[i]=b[i];
    fft(a,mid,opt);fft(a+mid,mid,opt);
    cp uni(coss[mid],opt*sinn[mid]),w(1,0);
    for(int k=0;k<len;k++)
    {
        if(k<mid)b[k]=a[k]+w*a[k+mid];
        else b[k]=a[k-mid]+w*a[k];
        w=w*uni;
    }
    for(int i=0;i<len;i++)a[i]=b[i];
}
void mul(cp *a,cp *b,int len)
{
    fft(a,len,1);
    fft(b,len,1);
    for(int i=0;i<len;i++)a[i]=a[i]*b[i];
    fft(a,len,-1);
    for(int i=0;i<len;i++)a[i]=a[i]/len;
}
string sa,sb,sans;
cp a[MAXN],b[MAXN];
int ans[MAXN];
int main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    cin>>sa>>sb;
    int n=sa.size(),m=sb.size();
    int len=1;
    while(len<n+m-1)coss[len]=cos(pi/len),sinn[len]=sin(pi/len),len<<=1;
    for(int i=0;i<len;i++)
    {
        a[i]=i<n?cp(sa[n-i-1]-'0',0):cp(0,0);
        b[i]=i<m?cp(sb[m-i-1]-'0',0):cp(0,0);
    }
    
    mul(a,b,len);
    
    for(int i=0;i<len;i++)
    {
        ans[i]+=round(a[i].s);
        ans[i+1]+=ans[i]/10;
        ans[i]%=10;
    }
    int flag=0;
    for(int i=len-1;i>=0;i--)
    {
        if(ans[i]==0&&flag==0)continue;
        else
        {
            flag=1;
            sans+=(char)(ans[i]+'0');
        }
    }
    cout<<sans<<endl;
    return 0;
}

非递归版

阅读世界，共赴山海

423全民读书节，邀你共读