快速傅里叶变换入门

什么是 FFT

先不说 FFT 在傅里叶分析等领域的运用，我们直接讨论此算法在信息学竞赛中的应用。

FFT 是 Fast Fourier Transformation 的缩写，其意为” 快速的傅里叶变换”。那么也就是说，原本有一个算法叫做 “傅里叶变换”，后来人们发明了一个更快的算法。其实，基本的离散傅里叶变换 (DFT) 的时间复杂度为 $O(n^2)$，而快速离散傅里叶变换 (FFT) 的时间复杂度为 $O(nlogn)$。

下面，我们将围绕 FFT 的基本过程和其内部机理展开介绍。本文中的 FFT 也许不太正宗，因为这是专门针对信息学竞赛的 FFT。

离散傅里叶变换 (DFT) 的作用

有时我们需要求解一类卷积问题：如离散型函数 F(x) 和 G(x) 的卷积：
$$
(F*G)(x)=\sum\limits_{i=0}^{n-1}F(i)G(i)
$$
直观的理解上述过程，我们会发现，这个过程很像小学数学中的竖式乘法。例如:
$$
\begin{matrix}
& 1 & 2 & 3\\
\times & & 2 & 3\\
\hline
& 3 & 6 & 9\\
2 & 4 & 6 & \\
\hline
2 & 3+4 & 6+6 & 9
\end{matrix}
$$
也就是说这类卷积的本质就是竖式乘法，也可以认为是多项式乘法。

那么说，我们如何利用手段将 DFT，也就是竖式乘法，多项式乘法，优化成 O(nlogn) 的复杂度呢？这或许需要利用一些高深的手段。

FFT 中的数学

• 将 F 和 G 转化为点值表达式—(DFT)
• 将 F 和 G 的点值表达式相乘得到 (F*G) 的点值表达式
• 通过 (F*G) 的点值表达式求出 (F*G) 的系数表达式—(IDFT)

那么点值表达式是什么呢？先看一个引理：

引理 1：恰当的 n 个点可以确定一个 n-1 次函数，每个 n-1 次函数都可以确定至少 n 个恰当点

证明 1：略

前面我们讲到，函数的卷积，就是竖式乘法，就是多项式乘法。因为一个值域有 n 个值的离散型函数一定可以用一个 n-1 次函数描述 (显然)，所以我们可以用一个 n-1 次的函数 (多项式) 去代替这个函数。而又因为这个函数和函数上恰当的 n 个点一一对应，所以我们又可以用 n 个点去代替这个函数。

点值表达式也就是用 $P_i(x_i,F(x_i))$代替 F，用 $Q_i(x_i,(Gx_i))$代替 G，这样一定有 $R_i(x_i,F(x_i)G(x_i))$在函数 $(F*G)$上。这也就是步骤 2。

理论准备 (公式预警)

如果我们枚举恰当的 n 个点去代替 F 和 G 的话，需要将每个 x 值带入 F 和 G 求值。运算量依旧是 $O(n^2)$。因此我们需要在 n 个点的取值上做做手脚，让它在代入时可以 $O(logn)$算出 f 和 G 的取值。

根据科学家们的智慧，我们直接写出结论：一些特殊的复数可以做到这一点。

引理 2： 复数是 $(x\cdot1+y\cdot i)$这样的数

证明 2： 略。当我们学习实数时，我们利用的是数轴。而这种有两个单位元的数，我们就可以用复

​ 平面来理解。复数 $(x,yi)$就可以看成复平面上 (x,y) 这个点。

引理 3： $e^{i\theta}=cos(\theta)+i\cdot sin(\theta)$(也就是著名的欧拉公式)

证明 3： (欧拉通过泰勒公式发现欧拉公式的过程)

​ 以下的公式根据 sin(x) 和 cos(x) 的泰勒展开得出，具体的证明和正确性我怎么可能清楚…
$$
\begin{align}
e^x &=\frac{1}{0!}x^0+\frac{1}{1!}x+\frac{1}{2!}x^2+\frac{1}{3!}x^3+\dots\\
cos(x) &=\frac{1}{0!}x^0-\frac{1}{2!}x^2+\frac{1}{4!}x^4-\dots\\
sin(x) &=\frac{1}{1!}x-\frac{1}{3!}x^3+\frac{1}{5!}x^5-\dots\\
\end{align}
$$
​ 将 $x=i\theta$代入 e 的式子得：
$$
\begin{align}
e^{i\theta} &=cos(\theta)+i\cdot sin(\theta)
\end{align}
$$
​ 也可以把$e^{i\theta}$理解成复平面上点 (1,0) 绕原点旋转$\theta$度后的结果。这样可以更方便地理解其他的引理。

我们需要的复数其实是 $x^n=1$的所有根。把 $e^{i\theta}$理解成复平面上点 (1,0) 绕原点旋转 $\theta$弧度后的结果后，我们可以知道这 n 个根 W 是：
$$
\large W_n^k=e^{\frac{2\pi i k}{n}} \text{(k=0,1…n-1)}
$$
下面我们分析以下 W 的特性。这些特性被用于步骤 1 和步骤 3。

引理 4： 消去引理
$$
\large W_{dn}^{dk}=W_n^k
$$
证明 4：
$$
\large W_{nk}^{dk}=e^{\frac{2\pi i dk}{dn}}=W_n^k
$$
引理 5： 周期引理
$$
\large W_n^k=W_n^{k+n}
$$
证明 5： 将欧拉公式理解为点的旋转

引理 6： 折半引理（重要）

​ 当 n 为偶数时有：
$$
\large (W_n^k)^2=W_{\frac{n}{2}}^k
$$
证明 6：
$$
\large (W_n^k)^2=(e^{\frac{2\pi i k}{n}})^2=e^{\frac{2\pi i k}{n/2}}=W_{n/2}^{k}
$$
引理 7： 求和引理
$$
\sum\limits_{j=0}^{n-1}(W_n^k)^j=0 (n>=1)
$$
证明 7： 运用等比数列求和公式
$$
\sum\limits_{k=0}^{n-1}(W_n^j)^k=W_n^j\frac{1-(W_n^j)^n}{1-W_n^j}=0
$$

步骤一：将 F 和 G 转化为点值表达式—DFT

在 DFT 变换中，我们希望计算函数 F(x) 在 $W_n^0,W_n^1\dots W_n^{n-1}$处的值，也就是
$$
A(k)=F(W_n^k)=\sum\limits_{j=0}^{n-1}a_jW_n^{kj}
$$
下面我们利用 “6: 折半引理” 二分计算 A(k)，使其达到 O(nlogn) 的复杂度。在此之前，我们先确保一个事情：a 序列的长度是偶数。然后我们就可以二分了。
$$
\begin{align}
令 A^{[0]}(x)&=a_0+a_2x+a_4x^2+a_6x^3+\dots+a_{n-2}x^{\frac{n}{2}-1}\\
令 A^{[1]}(x)&=a_1+a_3x+a_5x^2+a_7x^3+\dots+a_{n-1}x^{\frac{n}{2}-1}\\
于是有 A(x)&=A^{[0]}(x^2)+x\cdot A^{[1]}(x^2)
\end{align}
$$
也就是说：
$$
A(W_n^k)=A^{[0]}(W_{\frac{n}{2}}^k)+W_n^kA^{[1]}(W_{\frac{n}{2}}^{k})
$$

复杂度证明： 这样，当我们递归求解 A(x) 时，递归到第 d 层，数列就被分成了 $2^d$个子数列，每个子数列要代入的 W 都是 $\large W_{\frac{n}{2^d}}^k$，因此每层的状态数都是 $2^d\cdot \frac{n}{2^d}=n$，总层数又是 $logn$，所以总状态数就是 $nlogn$。这样我们就有了 $O(nlogn)$的递归算法。

步骤二：求 (F*G) 的点值表达式

直接将 F 和 G 的点值表达式的对应项相乘即可。
$$
C[i]=A[i]\cdot B[i]
$$

步骤三：根据 (F*G) 的点值表达式求 (F*G) 的系数表达式—IDFT

根据 DFT 得到的多项式点值表达式可以看做是一个矩阵和多项式的系数相乘

$C$为点值，$c$为系数。
$$
\begin{bmatrix}
C[0]\\
C[1]\\
C[2]\\
\vdots\\
C[n]
\end{bmatrix}=
\begin{bmatrix}
W_n^0 & W_n^0 & W_n^0 & \cdots & W_n^{0\cdot (n-1)}\\
W_n^0 & W_n^1 & W_n^2 & \cdots & W_n^{1\cdot (n-1)}\\
W_n^0 & W_n^2 & W_n^4 & \cdots & W_n^{2\cdot (n-1)}\\
\vdots & \vdots & \vdots & \ddots & \vdots \\
W_n^{n\cdot0} & W_n^{n\cdot 1} & W_n^{n\cdot 2} & \cdots & W_n^{(n-1)\cdot(n-1)}
\end{bmatrix}\times
\begin{bmatrix}
c_0\\
c_1\\
c_2\\
\vdots\\
c_{n-1}
\end{bmatrix}
$$
由于我们需要求最右边的矩阵，因此我们希望把上面的矩阵乘法变成这样的形式：
$$
C\times V^{-1}=c
$$
也就是说，我们要求出中间那个大矩阵 $V$的逆矩阵 $V^{-1}$，就可以求得系数矩阵 c。如果求得的逆矩阵也有着非常好的特性，我们就可以照样 $O(nlogn)$求出系数矩阵 c

引理 8： $V^{-1}$中第 j 行第 l 列的值为 $\large \frac{W_n^{-jl}}{n}$

证明 8： 观察 V 矩阵不难发现，$V[j][l]=W_n^{jl}$。根据逆矩阵的定义，$V\times V^{-1}$为单位矩阵。

​ 也就是说：
$$
\sum\limits_{k=0}^{n-1}V[j][k]V^{-1}[k][l]=
\begin{cases}
0 & (j\not= l)\\
1 & (j=l)\\
\end{cases}
$$
​ 而根据求和引理，$\large V^{-1}[j][l]=\frac{W_n^{-jl}}{n}$恰好满足上述限制。

​ 因为：
$$
\sum\limits_{k=0}^{n-1}(\frac{W_n^{-jk}}{n})(W_n^{kl})=\frac{1}{n}\sum\limits_{k=0}^{n-1}{W_n^{k(l-j)}}=\frac{1}{n}\sum\limits_{k=0}^{n-1}({W_n^{l-j}})^k
$$
​ 参考”7：求和引理” 即可证明这是对的。故这样的 $V^{-1}$的确是 V 的逆矩阵。

现在我们来找这个逆矩阵的” 很好的性质”

我们已经知道了这个 $V$矩阵的逆矩阵 $V^{-1}$的每一个元素的值，现在我们把 $V^{-1}$和 C 相乘的结果写出来：
$$
c_i=\frac{1}{n}\sum_{j=0}^{n-1}C_jW_n^{-ij}
$$
这个式子其实很像之前我们的 DFT 的公式：
$$
C_j=\sum_{i=0}^{n-1}c_iW_n^{ij}
$$
唯二的区别就在于多了个 $\frac{1}{n}$，以及 $W$的幂变成了负的。因此既然我们能够 $O(nlogn)$计算下面那个式子，我们也能 $O(nlogn)$计算上面那个式子。上面那个式子的计算就是 IDFT。

总结以下：

回到最初的目的上来，我们要求的是 (A*B) 的系数表达式。

如果我们默认 A 和 B 都是长度为 n 多项式，且 n 为 2 的幂。那么有：
$$
A*B=IDFT_{2n}(DFT_{2n}(A)\cdot DFT_{2n}(B))
$$
其中那个点表示向量的点乘，也就是说把 A 和 B 经过 DFT 后的点值对应相乘。

这也就是我们常说的 FFT 的流程。

FFT 算法实现

对数组进行奇偶划分

假设我们现在在计算多项式 $a[0]x^0+a[1]x^1+\cdots+a[2^n-2]x^{2^n-2}$的值。现在我们用 a 数组的每一项 a[i] 保存 x 取 $W_n^i$时多项式的取值。

一个简单的对这个多项式递归 FFT 的方式就是：计算 $a[0]x^0+a[2]x^2+\cdots+a[2^n-2]x^{2^n-2}$和 $a[1]x^1+a[3]x^3+\cdots+a[2^n-1]x^{2^n-1}$在 $x=W_{n/2}^i$处的取值，然后利用这两个多项式的共 n 个取值，通过上文的方法求出原多项式的 n 个取值。

具体的做法是先将偶数项放到数组左边，奇数项放到数组右边，递归计算出 a 数组的取值，其中 a[l~mid] 表示偶数项的 n/2 个取值，a[mid+1~r] 表示奇数项的 n/2 个取值。进一步，我们可以知道 (负号是由于 $W_n^{\frac{n}{2}}=-1​$而来的)：
$$
A(W_n^k)=A^{[0]}(W_{\frac{n}{2}}^k)+W_n^kA^{[1]}(W_{\frac{n}{2}}^{k}) (k<\frac{n}{2})\\
A(W_n^{k+\frac{n}{2}})=A^{[0]}(W_{\frac{n}{2}}^k)-W_n^kA^{[1]}(W_{\frac{n}{2}}^{k}) (k<\frac{n}{2})
$$
所以，我们用分为奇偶的 a[k] 和 a[mid+k] 就可以算出新的 a[k] 和 a[mid+k]。这个操作也叫一个蝶形运算。

注意到 FFT 时递归将 A 数组按奇偶划分为两个数列，其划分方式大致如下图所示：

比如 a0~a7 的序列经过如上的划分得到 0,4,2,6,1,5,3,7 的序列。也就是其下标的二进制数位倒转后的升序排列。所以我们可以事先将系数数组 a 排成划分后的样子，用 FFT 自底向上计算，结果就是取值数组 a 的最终样子。

下面给出 FFT 求两个十进制数乘法的代码，细节全部在代码中体现。

#include <algorithm>
#include <iostream>
#include <cstring>
#include <cstdio>
#include <cmath>
#define pi acos(-1)
#define MX 262145
#define reg register

using namespace std;

typedef struct com_t        //定义实数类
{
    double r,i;
    inline com_t operator + (const com_t t)const{return (com_t){r+t.r,i+t.i};}
    inline com_t operator - (const com_t t)const{return (com_t){r-t.r,i-t.i};}
    inline com_t operator * (const com_t t)const{return (com_t){r*t.r-i*t.i,r *t.i+i*t.r};}
}com;
com A[MX],B[MX],C[MX];
int ord[MX];
int n,m,l;                  //分别为十进制数补为 2 的幂后的长度、长度之和，将长度补为 2 的几次幂
inline void dft(com *a,int f,int len)   //f=-1 表示求 IDFT
{
    for(reg int i=0;i<len;i++)if(i<ord[i])swap(a[i],a[ord[i]]);
    for(reg int i=1;i<len;i<<=1)
    {
        com wn=(com){cos(pi/i),f*sin(pi/i)};    //将 2pi/2i 约分后的结果
        for(reg int j=0;j<len;j+=i<<1)
        {
            com w=(com){1,0};
            for(reg int k=0;k<i;k++,w=w*wn)
            {
                com x=a[j+k],y=w*a[j+k+i];      //因为原址 FFT 会修改数组的值，必须先将值取出
                a[j+k]=x+y;
                a[j+k+i]=x-y;
            }
        }
    }
    if(f==-1)for(reg int i=0;i<len;i++)a[i].r/=len;
}

inline void fft(com *a,com *b,com *c)   //fft 的流程
{
    dft(a,1,n);
    dft(b,1,n);
    for(reg int i=0;i<n;i++)c[i]=a[i]*b[i];
    dft(c,-1,n);
}

void output(com *c)
{
    int ans[MX],st=0;
    for(int i=0;i<n;i++)ans[i]=int(c[i].r+0.5);
    for(int i=0;i<n;i++)ans[i]>=10?ans[i+1]+=ans[i]/10,ans[i]%=10:0;
    for(int i=n-1;i>=0;i--)if(ans[i]){st=i;break;}
    for(int i=st;i>=0;i--)printf("%d",ans[i]);
}

void init()
{
    int l1,l2;
    char str1[MX],str2[MX];
    scanf("%*d");
    scanf("%s%s",str1,str2);
    l1=strlen(str1),l2=strlen(str2);
    reverse(str1,str1+l1);
    reverse(str2,str2+l2);
    for(int i=0;i<l1;i++)A[i].r=str1[i]-'0';
    for(int j=0;j<l2;j++)B[j].r=str2[j]-'0';
    m=l1+l2;                    //求长度之和
    for(n=1;n<m;n<<=1)l++;      //补为 2 的幂
    for(int i=0;i<n;i++)ord[i]=(ord[i>>1]>>1)|((i&1)<<l-1); //求每个二进制数反转二进制位后的值
}

int main()
{
    init();
    fft(A,B,C);
    output(C);
    return 0;
}