FFT的学习

来源：互联网发布：程序员必备生活装备编辑：程序博客网时间：2024/04/27 19:51

本文来自转载：

FFT算法学习笔记

Posted on July 14, 2015 by gatevin

写在前面

关于学习FFT算法的资料个人最推荐的还是算法导论上的第30章(第三版), 多项式与快速傅里叶变换, 基础知识都讲得很全面

这篇文章即是博主在阅读FFT算法的相关学习资料之后的一些总结, 作为日后复习查阅的材料来源, 同时也可以作为读者学习FFT算法的一份资料

FFT算法基本概念

FFT(FastFourierTransformation)

即快速傅里叶变换, 是离散傅里叶变换的加速算法, 可以在O(nlogn)的时间里完成DFT, 利用相似性也可以在同样复杂度的时间里完成逆DFT

DFT(DiscreteFourierTransform)

即离散傅里叶变换, 这里主要就是多项式的系数向量转换成点值表示的过程

在ACM-ICPC竞赛中, FFT算法常被用来为多项式乘法加速, 即在O(nlogn)

复杂度内完成多项式乘法, 当然实际应用不仅仅限于这些, 时常出现需要构造多项式相乘来进行计数的问题, 也需要用FFT算法来解决, 相关的几个问题在本文中也会提及

FFT算法需要的基础数学知识

多项式相关：

多项式相关的定义：一个以x

为变量的多项式定义在代数域F上, 将函数A(x)

表示为形式和

Σ j = 0 n - 1 a j x j

, 称aj为系数, 所有系数属于代数域F, 如果一个多项式的最高非零系数是ak, 那么成这个多项式的次数是k, 记做degree(A)=k, 任何一个严格大于一个多项式次数的整数都是该多项式的次数界

关于多项式的加法和乘法相信看到这篇博客的读者都会最基本的中学的算法, 在计算的时候, 如果采用传统的中学的计算方法, 多项式加法的时间复杂度是O(n)

, 乘法的时间复杂度是O(n2)

(n是两个多项式A和B的次数).

多项式的表示

在平常的学习中, 最常见的是多项式的系数表达方式,

次数界为n的多项式A(x)=a0+a1x+a2x2+...+an−1xn−1

的系数表达为一个由系数组成的向量a=(a1,a2,...an−1)

但是多项式还有一个比较常用的表示方法, 即多项式的点值表达方式

一个次数界为n的多项式A(x)=a0+a1x+a2x2+...+an−1xn−1

的点值表达就是一个由n个点对组成的集合

{(x 0, y 0), (x 1, y 1), . . ., (x n - 1, y n - 1)}

使得对任意的整数k=0,1,..n−1,xk各不相同, 且yk=A(xk)

关于点值表达的正确性证明：

对于任意n个点值对组成的集合{(x0,y0),(x1,y1),...,(xn−1,yn−1)}

, 如果存在一个次数界为n的多项式A(x)

过这n个点, 那么

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ 11 ⋮ 1 x 0 x 1 ⋮ x n - 1 x 20 x 21 ⋮ x 2 n - 1 \dots \dots ⋱ \dots x n - 1 0 x n - 1 1 ⋮ x n - 1 n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ a 0 a 1 ⋮ a n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ y 0 y 1 ⋮ y n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥

最左边这个n*n的矩阵称为范德蒙德矩阵, 可以用数学归纳法证明它的行列式值为

Π 0 \leq j < k \leq n - 1 (x k - x j)

当xk两两不相同时明显这个行列式的值不为0, 该矩阵可逆, 于是存在唯一解, 所以多项式的点值表达是合理的

相应的通过n个点的坐标直接确定多项式各个系数的值的方法是存在的, 感兴趣的读者可以查询拉格朗日公式的相关资料, 利用拉格朗日插值公式可以在O(n2)

的时间复杂度内得到多项式的系数表达, 这个也是算法导论中的一个习题

点值表达方式下多项式的乘法, 不难发现, 如果多项式A(x),B(x)

的点值表示分别是

{(x 0, y a 0), (x 1, y a 1), . . ., (x n - 1, y a n - 1)}

和

{(x 0, y b 0), (x 1, y b 1), . . ., (x n - 1, y b n - 1)}

那么如果多项式C(x)=A(x)B(x), 那么C(x)的点值表达是

{(x 0, y a 0 y b 0), (x 1, y a 1 y b 1), . . ., (x n - 1, y a n - 1 y b n - 1)}

卷积

对于两个多项式的系数向量a=(a0,a1,..,an−1)

和b=(b0,b1,..,bn−1), 两个多项式相乘得到的多项式的系数向量c=(c0,c1,..,c2n−2)满足cj=Σjk=0akbj−k, 称系数向量c是输入向量a和b的卷积, 记作c=a⊗b

在简单的多项式乘法的计算方法中, 每一个多项式的系数都通过系数表示方式下卷积的方式来进行计算, 时间复杂度是O(n2)

, 但是FFT是先将多项式的从系数表示法转换成点值表示法(可以在O(nlogn)的时间复杂度下完成, 也就是加速的DFT变换, 然后在点值表示法下进行乘积计算, 在O(n)的时间复杂度内得到结果的点值表示法, 然后进行逆DFT变换, 在O(nlogn)

的时间复杂度下完成逆DFT变换得到系数表示法

而要理解DFT, 则需要一定复数上的数学知识：

复数相关的基础知识：

单位复数根：n次单位复数根指的是满足ωn=1

的所有复数ω, n次单位复数根刚好有n个, 他们是e2kπni, 其中i是复数单位, k=0,1,2...n−1

, 在复平面上这n个根均匀的分布在半径为1的圆上, 关于复数指数的定义如下：

e u i = c o s (u) + s i n (u) i

其中ωn=e2iπn倍称为主n次单位根(这个定义好像接下来没用到= =)

关于复数根的几个定理和引理：

消去引理：对任何整数 n≥0,k≥0,d≥0

有ωdkdn=ωkn

证明： ωdkdn=(e2iπdn)dk=(e2iπn)k=ωkn

一个推论：对任意偶数 n > 0 有 ωn2n=ω2=−1

证明：设n = 2*k那么 ωn2n=ωk2k=ω2=eπ=cos(π)+sin(π)i=−1

折半引理：如果n > 0是偶数, 那么n个n次单位复数根的平方的集合就是n/2个n/2次单位复数根的集合

证明：实际上这个引理就是证明了(ωk+n2n)2=ω2k+nn=ω2knωnn=(ωkn)2

折半引理对于采用分治对多项式系数向点值表达的转换有很大作用, 保证了递归的子问题是原问题规模的一半

求和引理：对任意整数n≥1

和不能被n整除的非负整数k, 有

Σ j = 0 n - 1 (ω k n) j = 0

这个问题通过等比数列求和公式就可以得到：Σn−1j=0(ωkn)j=(ωkn)k−1ωkn−1=1k−1ωkn−1=0

DFT与FFT, 以及逆DFT

DFT

在DFT变换中, 希望计算多项式A(x)在复数根ω0n,ω1n,...,ωn−1n

处的值, 也就是求

y k = A (ω k n) = Σ j = 0 n - 1 a j ω k j n

称向量y=(y0,y1,...,yn−1)是系数向量a=(a0,a1,...,an−1)的离散傅里叶变换, 记为y=DFTn(a)

FFT

直接计算DFT的复杂度是O(n2)

, 而利用复数根的特殊性质的话, 可以在O(nlogn)

的时间内完成, 这个方法就是FFT方法, 在FFT方法中采用分治策略来进行操作, 主要利用了消去引理之后的那个推论

在FFT的策略中, 多项式的次数是2的整数次幂, 不足的话再前面补0, 每一步将当前的多项式A(x), 次数是2的倍数, 分成两个部分：

A [0] (x) = a 0 + a 2 x + a 4 x 2 + . . . + a n - 2 x n 2 - 1

A [1] (x) = a 1 + a 3 x 1 + a 5 x 2 + . . . + a n - 1 x n 2 - 1

于是就有了

A (x) = A [0] (x 2) + x A [1] (x 2)

那么我们如果能求出次数界是n2的多项式A[0](x)和A[1](x)在n个n次单位复数根的平方处的取值就可以了, 即在

(ω 0 n) 2, (ω 1 n) 2, (ω 2 n) 2, . . ., (ω n - 1 n) 2

处的值, 那么根据折半引理, 这n个数其实只有n2个不同的值, 也就是说, 对于每次分出的两个次数界n2的多项式, 只需要求出其n2个不同的值即可, 那么问题就递归到了原来规模的一半, 也就是说如果知道了两个子问题的结果, 当前问题可以在两个子问题次数之和的复杂度内解决, 那么这样递归问题的复杂度将会是O(nlogn)的, 用a=(a0,a1,...,an−1)表示系数向量,y=(y0,y1,...,yn−1)表示离散变换之后的向量, 这里给出将算导上的代码翻译出来的C++代码实现(以解决算法第三版大论第三十章的一个习题, 求(0, 1, 2, 3)的DFT为例)：

/* * Author: Gatevin * Created Time:  2015/7/12 20:08:47 * File Name: FFT.cpp */#include<iostream>#include<sstream>#include<fstream>#include<vector>#include<list>#include<deque>#include<queue>#include<stack>#include<map>#include<set>#include<bitset>#include<algorithm>#include<cstdio>#include<cstdlib>#include<cstring>#include<cctype>#include<cmath>#include<ctime>#include<iomanip>using namespace std;const double eps(1e-8);typedef long long lint;const double PI = acos(-1.0);/* * 这是一个递归实现FFT的测试, 测试习题中求DFT(0, 1, 2, 3) */struct Complex{    double real, image;    Complex(double _real, double _image)    {        real = _real;        image = _image;    }    Complex(){}};Complex operator + (const Complex &c1, const Complex &c2){    return Complex(c1.real + c2.real, c1.image + c2.image);}Complex operator - (const Complex &c1, const Complex &c2){    return Complex(c1.real - c2.real, c1.image - c2.image);}Complex operator * (const Complex &c1, const Complex &c2){    return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real);}Complex* RecursiveFFT(Complex a[], int n)//n表示向量a的维数{    if(n == 1)        return a;    Complex wn = Complex(cos(2*PI/n), sin(2*PI/n));    Complex w = Complex(1, 0);    Complex* a0 = new Complex[n >> 1];    Complex* a1 = new Complex[n >> 1];    for(int i = 0; i < n; i++)        if(i & 1) a1[(i - 1) >> 1] = a[i];        else a0[i >> 1] = a[i];    Complex *y0, *y1;    y0 = RecursiveFFT(a0, n >> 1);    y1 = RecursiveFFT(a1, n >> 1);    Complex* y = new Complex[n];    for(int k = 0; k < (n >> 1); k++)    {        y[k] = y0[k] + w*y1[k];        y[k + (n >> 1)] = y0[k] - w*y1[k];        w = w*wn;    }    delete a0;    delete a1;    delete y0;    delete y1;    return y;}int main(){    Complex* a = new Complex[10];    a[0] = Complex(0, 0);    a[1] = Complex(1, 0);    a[2] = Complex(2, 0);    a[3] = Complex(3, 0);    Complex* ans = new Complex[10];    ans = RecursiveFFT(a, 4);    for(int i = 0; i < 4; i++)        cout<<ans[i].real<<"+"<<ans[i].image<<"i"<<endl;    delete a;    delete ans;    return 0;}

100

101

102

* Author: Gatevin

* Created Time: 2015/7/12 20:08:47

* File Name: FFT.cpp

#include<iostream>

#include<sstream>

#include<fstream>

#include<vector>

#include<list>

#include<deque>

#include<queue>

#include<stack>

#include<map>

#include<set>

#include<bitset>

#include<algorithm>

#include<cstdio>

#include<cstdlib>

#include<cstring>

#include<cctype>

#include<cmath>

#include<ctime>

#include<iomanip>

usingnamespace std;

const doubleeps(1e-8);

typedeflong longlint;

constdouble PI= acos(-1.0);

* 这是一个递归实现FFT的测试, 测试习题中求DFT(0, 1, 2, 3)

struct Complex

{

doublereal,image;

Complex(double_real,double _image)

{

real= _real;

image= _image;

}

Complex(){}

};

Complexoperator +(constComplex &c1,const Complex&c2)

{

returnComplex(c1.real+ c2.real,c1.image+ c2.image);

}

Complex operator- (constComplex &c1,const Complex&c2)

{

returnComplex(c1.real- c2.real,c1.image- c2.image);

}

Complexoperator * (constComplex &c1,const Complex&c2)

{

returnComplex(c1.real*c2.real- c1.image*c2.image,c1.real*c2.image+ c1.image*c2.real);

}

Complex* RecursiveFFT(Complexa[],int n)//n表示向量a的维数

{

if(n== 1)

returna;

Complex wn= Complex(cos(2*PI/n),sin(2*PI/n));

Complexw =Complex(1,0);

Complex*a0 =new Complex[n>> 1];

Complex*a1 =new Complex[n>> 1];

for(inti =0;i <n;i++)

if(i& 1)a1[(i- 1)>> 1]= a[i];

elsea0[i>> 1]= a[i];

Complex *y0,*y1;

y0= RecursiveFFT(a0,n >>1);

y1= RecursiveFFT(a1,n >>1);

Complex*y =new Complex[n];

for(intk =0;k <(n>> 1);k++)

{

y[k]= y0[k]+ w*y1[k];

y[k+ (n>> 1)]= y0[k]- w*y1[k];

w= w*wn;

}

deletea0;

delete a1;

deletey0;

delete y1;

returny;

}

int main()

{

Complex*a =new Complex[10];

a[0]= Complex(0,0);

a[1]= Complex(1,0);

a[2]= Complex(2,0);

a[3]= Complex(3,0);

Complex*ans =new Complex[10];

ans= RecursiveFFT(a,4);

for(inti =0;i <4;i++)

cout<<ans[i].real<<"+"<<ans[i].image<<"i"<<endl;

deletea;

delete ans;

return0;

}

可以得到系数向量(0, 1, 2, 3)经过DFT之后是(6, -2-2i, -2, -2+2i)

在上面这个视线中需要注意的就是RecursiveFFT中的y[k] = y0[k] + w*y1[k];和y[k + (n >> 1)] = y0[k] - w*y1[k];的变化, 正是利用了对于偶数的n，有ωn2n=−1

成立, 得到数组y的所有值

另外有一个概念性的东西： yk=y[0]k+ωkny[1]k

和yk+n2=y[0]k−ωkny[1]k都成立, 在这其中将ωkn

称为旋转因子

在进行了DFT操作之后, 将多项式的系数表达成功转为了点值表达接下来是如何在O(nlogn)

的时间内完成逆DFT的问题, 通过逆DFT将点值表达还原为系数表达

逆DFT

根据DFT得到的向量y和系数向量a之间的关系, 可以用矩阵乘积的形式来表达他们之间的关系, 即y=Vna

, 也就是

⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ y 0 y 1 y 2 y 3 ⋮ y n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ = ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ 1111 ⋮ 1 1 ω n ω 2 n ω 3 n ⋮ ω n - 1 n 1 ω 2 n ω 4 n ω 6 n ⋮ ω 2 (n - 1) n 1 ω 3 n ω 6 n ω 9 n ⋮ ω 3 (n - 1) n \dots \dots \dots \dots ⋱ \dots 1 ω n - 1 n ω 2 (n - 1) n ω 3 (n - 1) n ⋮ ω (n - 1) (n - 1) n ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎡ ⎣ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ ⎢ a 0 a 1 a 2 a 3 ⋮ a n - 1 ⎤ ⎦ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥ ⎥

那么要将DFT变化得到的向量y还原成向量a的话, 只需要用Vn

的逆矩阵V−1n

乘上向量y即可, 这里需要用到一个定理

在矩阵Vn

中, 不难发现对于任意的j,k=0,1,...,n−1,Vn(k,j)=ωkjn, 那么可以找到这样一个矩阵V−1n, 对于任意的j,k=0,1,...,n−1,V−1n的(j,k)出的元素为ω−kjnn,V−1n是矩阵Vn

的逆矩阵

证明如下：要证明这两个矩阵互逆, 证明其积为单位矩阵即可, 考虑两个矩阵的乘积在(j,j′)

出的元素, 可以发现这个元素是

Σn−1k=0(ω−kjnn)(ωkj′n)=

Σn−1k=0ωk(j′−j)nn

当j′=j

时, 这个和是1, 否则根据求和引理, 这个和是0, 故这两个矩阵互逆

那么根据这个逆矩阵可以发现要计算DFT−1n(y)

的话, 有关系式 aj=1nΣn−1k=0ykω−kjn比较这个式子和之前DFT里面y和a的关系式子, 可以发现只需要用ω−1n替换掉ωn即可, 最后结果需要除以n, 所以计算逆DFT的方法和计算DFT和相似, 都可以在O(nlogn)

的时间复杂度内解决

卷积定理

对任意两个长度为n的向量a和b, 其中n是2的幂, 有

a \otimes b = D F T - 1 2 n (D F T 2 n (a) \cdot D F T 2 n (b))

其中向量a和b用0填充, 使其长度达到2n, 并用⋅表示两个2n个元素组成向量的点乘(也就是每一维上的数相乘)

这个式子实际上就是多项式的系数表达在乘法时进行的卷积运算得到的结果, 等同于通过将其系数进行DFT变换变成点值表达之后相乘再换回来的过程

关于FFT算法的迭代实现

在递归实现DFT过程的FFT算法中, 我们每次将系数向量a分成两个部分利用折半引理来降低计算的规模, 可以发现在每次分组当中他们满足这样一个完全二叉树的分组(n是2的幂)：

FFT递归系数分组

通过上图的流程可以看出, 最后一层的子节点下标的顺序实际上就是其下标转换成二进制串的倒序的字符串按照字典序排列的顺序

比如a0~a7得到的序列a0, a4, a2, a6, a1, a5, a3, a7下标的二进制是000, 100, 010, 110, 001, 101, 011, 111对应的串的倒序是000, 001, 010, 011, 100, 101, 110, 111这个倒序的二进制刚好是0, 1, 2, 3, 4, 5, 6, 7的二进制表示, 于是我们可以在O(nlogn)

的复杂度内得到做下面一层的下标顺序, 然后可以根据子节点的结果向上迭代得到父亲结点的值, 这样计算的话直接避免了递归, 如果直接递归的话在一些OJ上可能会造成爆栈的错误, 所以还是采用迭代的方式进行比较好.

关于迭代形式的FFT算法, C++代码实现如下(同样以求(0, 1, 2, 3)的DFT变换为例)：

这段代码的话同时也进行了逆DFT, DFT和逆DFT的过程相似, 加上一个标记判断当前执行的是哪一种就行了

/* * Author: Gatevin * Created Time:  2015/7/14 14:15:27 * File Name: IterativeFFT.cpp */#include<iostream>#include<sstream>#include<fstream>#include<vector>#include<list>#include<deque>#include<queue>#include<stack>#include<map>#include<set>#include<bitset>#include<algorithm>#include<cstdio>#include<cstdlib>#include<cstring>#include<cctype>#include<cmath>#include<ctime>#include<iomanip>using namespace std;const double eps(1e-8);typedef long long lint;const double PI = acos(-1.0);struct Complex{    double real, image;    Complex(double _real, double _image)    {        real = _real;        image = _image;    }    Complex(){}};Complex operator + (const Complex &c1, const Complex &c2){    return Complex(c1.real + c2.real, c1.image + c2.image);}Complex operator - (const Complex &c1, const Complex &c2){    return Complex(c1.real - c2.real, c1.image - c2.image);}Complex operator * (const Complex &c1, const Complex &c2){    return Complex(c1.real*c2.real - c1.image*c2.image, c1.real*c2.image + c1.image*c2.real);}int rev(int id, int len){    int ret = 0;    for(int i = 0; (1 << i) < len; i++)    {        ret <<= 1;        if(id & (1 << i)) ret |= 1;    }    return ret;}//当DFT= 1时是DFT, DFT = -1则是逆DFTComplex* IterativeFFT(Complex* a, int len, int DFT)//对长度为len(2的幂)的数组进行DFT变换{    Complex* A = new Complex[len];//用A数组存储数组a分组之后新的顺序    for(int i = 0; i < len; i++)        A[rev(i, len)] = a[i];    for(int s = 1; (1 << s) <= len; s++)    {        int m = (1 << s);        Complex wm = Complex(cos(DFT*2*PI/m), sin(DFT*2*PI/m));        for(int k = 0; k < len; k += m)//这一层结点的包含数组元素个数都是(1 << s)        {            Complex w = Complex(1, 0);            for(int j = 0; j < (m >> 1); j++)//折半引理, 根据两个子节点计算父亲节点            {                Complex t = w*A[k + j + (m >> 1)];                Complex u = A[k + j];                A[k + j] = u + t;                A[k + j + (m >> 1)] = u - t;                w = w*wm;            }        }    }    if(DFT == -1) for(int i = 0; i < len; i++) A[i].real /= len, A[i].image /= len;    return A;}int main(){    Complex* a = new Complex[4];    a[0] = Complex(0, 0);    a[1] = Complex(1, 0);    a[2] = Complex(2, 0);    a[3] = Complex(3, 0);    a = IterativeFFT(a, 4, 1);    cout<<"----------After DFT----------"<<endl;    for(int i = 0; i < 4; i++)        printf("%.9f + %.9fi\n", a[i].real, a[i].image);    cout<<"----------After DFT-1----------"<<endl;    a = IterativeFFT(a, 4, -1);    for(int i = 0; i < 4; i++)        printf("%.9f + %.9fi\n",a[i].real, a[i].image);    delete a;    return 0;}

100

101

102

103

104

105

106

107

108

109

110

111

112

* Author: Gatevin

* Created Time: 2015/7/14 14:15:27

* File Name: IterativeFFT.cpp

#include<iostream>

#include<sstream>

#include<fstream>

#include<vector>

#include<list>

#include<deque>

#include<queue>

#include<stack>

#include<map>

#include<set>

#include<bitset>

#include<algorithm>

#include<cstdio>

#include<cstdlib>

#include<cstring>

#include<cctype>

#include<cmath>

#include<ctime>

#include<iomanip>

usingnamespace std;

const doubleeps(1e-8);

typedeflong longlint;

constdouble PI= acos(-1.0);

structComplex

{

doublereal,image;

Complex(double_real,double _image)

{

real= _real;

image= _image;

}

Complex(){}

};

Complex operator+ (constComplex &c1,const Complex&c2)

{

returnComplex(c1.real+ c2.real,c1.image+ c2.image);

}

Complexoperator -(constComplex &c1,const Complex&c2)

{

returnComplex(c1.real- c2.real,c1.image- c2.image);

}

Complex operator *(constComplex &c1,const Complex&c2)

{

returnComplex(c1.real*c2.real- c1.image*c2.image,c1.real*c2.image+ c1.image*c2.real);

}

intrev(intid,int len)

{

intret =0;

for(inti =0;(1<< i)< len;i++)

{

ret<<=1;

if(id& (1<< i))ret |=1;

}

returnret;

}

//当DFT= 1时是DFT, DFT = -1则是逆DFT

Complex*IterativeFFT(Complex*a,int len,int DFT)//对长度为len(2的幂)的数组进行DFT变换

{

Complex*A =new Complex[len];//用A数组存储数组a分组之后新的顺序

for(inti =0;i <len;i++)

A[rev(i,len)]= a[i];

for(ints =1;(1<< s)<= len;s++)

{

intm =(1<< s);

Complexwm =Complex(cos(DFT*2*PI/m),sin(DFT*2*PI/m));

for(intk =0;k <len;k +=m)//这一层结点的包含数组元素个数都是(1 << s)

{

Complexw =Complex(1,0);

for(intj =0;j <(m>> 1);j++)//折半引理, 根据两个子节点计算父亲节点

{

Complext =w*A[k+ j+ (m>> 1)];

Complexu =A[k+ j];

A[k+ j]= u+ t;

A[k+ j+ (m>> 1)]= u- t;

w= w*wm;

}

if(DFT== -1)for(inti =0;i <len;i++)A[i].real/= len,A[i].image/= len;

returnA;

}

intmain()

{

Complex*a =new Complex[4];

a[0]= Complex(0,0);

a[1]= Complex(1,0);

a[2]= Complex(2,0);

a[3]= Complex(3,0);

a= IterativeFFT(a,4,1);

cout<<"----------After DFT----------"<<endl;

for(inti =0;i <4;i++)

printf("%.9f + %.9fi\n",a[i].real,a[i].image);

cout<<"----------After DFT-1----------"<<endl;

a= IterativeFFT(a,4,-1);

for(inti =0;i <4;i++)

printf("%.9f + %.9fi\n",a[i].real,a[i].image);

deletea;

return0;

}

通过上面这个代码的示例, FFT算法的实现基本没有什么问题了, 另外算法导论中的习题有一些很不错, 便于熟悉这一算法的很多细节, 这里就不一一提及了

经过这些学习之后, 进行一些实战演练是很有必要的, 接下来是相关习题的练习部分

相关练习：

HDU 1402 A*B Problem Plus 大整数乘法, 我的解答：戳我

HDU 4609 3-idiots FFT计数, 解答：戳我

UVALive 4671 K-neighbor Substrings FFT算字符串Hamming距离解答：戳我

UVA 12298 Super Poker II FFT计数, long double, 解答：戳我

URAL 1996 Cipher Massage 3 FFT + KMP 解答：戳我

CodeChef COUNTARI Arithmetic Progressions FFT + 分块, 解答：戳我

ZOJ 3856 Goldbach FFT计数, 解答：戳我

UVALive 6886 Golf Bot FFT模板题, 解答：戳我

HDU 4093 Xavier is Learning to Count 容斥原理 + FFT,(2011年上海现场赛C题) 解答：戳我

0 0