快速傅里叶变换入门

什么是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数组按奇偶划分为两个数列,其划分方式大致如下图所示:

FFT系数划分方式

比如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;
}
分类: 文章

2 条评论

蒟蒻XZY · 2017年11月19日 12:12 下午

卧槽绝对的本站精品啊!!!!
Orz真是感谢Boshi大佬
看起来很舒服~
可惜不能置顶。。。
本蒟蒻有空还是要好好研读一下,现在还是好好搞常规算了。。。

    yyb · 2017年12月9日 11:46 上午

    orz

发表评论

电子邮件地址不会被公开。 必填项已用*标注

你是机器人吗? =。= *