zkx06111's Blog

FFT学习笔记

| Comments

一直听说FFT很厉害,有许多数学题的结尾都是“然后就可以FFT了,相信大家都会吧。”但我一直处于不明觉厉的状态,只是抄过Frank_c1博客上的一段代码。昨天,我在Miskcoo的博客上看到一篇题为从多项式乘法到快速傅里叶变换的文章,终于对它有了一点理解。

多项式的表示方法

对于多项式,可以用两种方法表示。
的系数向量表示,叫做系数表示。
或者用f(x)在n个位置的取值,,叫做点值表示。
我们通常看到的多项式,都是系数表示的,这样的话要做多项式乘法需要的复杂度。
而如果是点值表示,就可以直接对于同一个,求出,就是结果的点值表示。这样复杂度是的。
所以,只要能够较快地实现点值表示和系数表示的转化,就可以提高系数表示的多项式做乘法的速度。

单位根

对于,把复数称为次单位根。
次单位根共有个,它们恰好位于复平面的单位元上,并将其分成个部分,根据复数乘法,模相乘,幅角相加,那么第个单位根
,那么n个单位根就是
在转化为点值表示时把单位根代入原多项式。

单位根的性质

DFT

把单位根代入多项式求值的过程。暴力是的。

把它分成次数为奇和为偶的,再结合单位根的性质,对于可以得到


这样就把问题规模缩小了一半。而每次合并子问题是的,所以总复杂度是的。

DFT的递归实现

我比较弱,所以不会非递归实现...在UOJ上通过多项式乘法的时间是1786ms,比frank_c1博客上的非递归版本慢了一秒多,太辣鸡了。

首先要预处理出单位根。

init
void init(int n){
    double pi=acos(-1);
    REP(i,0,n-1) W[i]=(cpx){cos(2.0*pi*i/n),sin(2.0*pi*i/n)};
    REP(i,0,n-1) CW[i]=(cpx){W[i].x,-W[i].y};
}

然后根据前面描述的做。

dft
void DFT(cpx \*A,int n,int start,int step){
    if (n==1) return;
    DFT(A,n>>1,start,step<<1); // even case
    DFT(A,n>>1,start+step,step<<1); // odd case
    REP(k,0,n/2-1){
        int pos=2*k*step;
        T[k]=A[pos+start]+W[k*step]*A[pos+start+step];
        T[k+n/2]=A[pos+start]-W[k*step]*A[pos+start+step];
    }
    REP(k,0,n-1) A[k*step+start]=T[k];
}

IDFT

这里的I是指Inversion,由多项式的点值表示得到其系数表示。
这本质上是一个求解方程组的过程,可以把方程组用矩阵表示出来。

把这个式子记成
把左边的矩阵记为,考虑它的逆元。

那么,把与相乘后得到的矩阵记为。最后的系数矩阵可以表示出来。

下面重点来了【话说这也是miskcoo博客上没写,导致我一开始看不懂的部分。】
这其实相当于,把求出来的点值的当作系数,把原来单位根的共轭复数当成单位根代入,再做一遍DFT。
感觉实在太神辣。

未完待续

理解了FFT之后,当然要去学习各种数学姿势,比如多项式的各种运算,还有生成函数,毕竟应用才是最重要的。

Comments

comments powered by Disqus