快速傅里叶变换FFT理解及其c++实现

快速傅里叶变换(FFT)是离散傅里叶变换DFT(Discrete Fourier Transform)的一种高效算法

时间复杂度

FFT

代码复杂度

Prerequisite

*线性代数

↑复变函数

reference

kuangbin模板

问题描述

给你两个多项式函数f(x)g(x)的系数表达式,f(x)*g(x)的系数表达式,

e.g.



分析



如果选的n+1个点都相同,则求 f(x)*g(x) 的点值表达式
只需要算(n+1)个值的乘积;





DFT

DFT过程对于两个多项式选取的n+1个点相同




FFT通过取某些特殊的 x
来加速DFT


命名方式:从实数1开始,逆时针依次命名






FFT使用点值表达式


并用分治方法处理

分治divide and conquer



方便起见,下文的n均是补全后的n,比如f5补全后,n-1次多项式,n=8







分治的具体过程:





  1. 如果把上面视作四个序列的话(忽略括号)

第一个序列的
索引
二进制依次为
000 001 010 011 100 101 110 111
最后一个序列
索引
二进制依次为
000 100 010 110 001 101 011 111

注意到上面的索引是翻转的,

如果采用最后一个序列,数组X的变换方法,见下面代码

template
void change(T (&x)[N]){

int i=1;//i是索引 001

int j=N/2;//ji二进制反过来的索引 100

while(i<N-1){

if(i<j) swap(x[i],x[j]);

int k=N/2;

while(j>= k){

j-=k;

k/=2;

}

j+=k;

i++;

}

}

*原理探究

IDFT

现在我们通过DFT得到了一个n项的方程组

设系数ai的矩阵为A

用矩阵表示 A是下面的增广矩阵的解


记上面增广矩阵的形式为 [Ω|Y]

Y=ΩX


*现成的结论逆矩阵形式如下



*试图证明

所以说IDFE的操作是得到上面的逆矩阵,可以说与DFT的操作过程完全相同,只是一个传入π,一个传入并且/n

总代码

//仅供学习,不是acm模板

//代码中的N就是上文中提到的
补全后的n

// m.h 文件参见
https://gist.github.com/Haozun/ad7520ee9f3468b19b3618bc12942adf

namespace fft{

typedef complex dcomplex;

template

void change(dcomplex (&X)[N]){//将序列x转变成另一个序列

int i=1;

int j=N/2;

while(i<N-1){

if(i<j) swap(X[i],X[j]);

int k=N/2;

while(j>= k){

j-=k;

k/=2;

}

j+=k;

i++;

}

}

template

void fft(dcomplex (&X)[N],int anti){

change(X);

for(int h=2;h <= N;h<<=1){//分治过程

dcomplex w_n=polar(1.0,-anti*2*PI/h);

for(int j=0;j<N;j+=h){

dcomplex w=1;

for(int k=j;k<j+h/2;k++){

dcomplex u=X[k];

dcomplex t=w*X[k+h/2];

X[k]=u+t;

X[k+h/2]=u-t;

w*=w_n;

}

}

}

if(anti)for(auto e:X){e/=N;}

}

}

 

Test1


#include “m.h”
using namespace FFT;

int main(){

int len=pow2(1+ceil(log2(max(5,3))));//补全

dcomplex X1[len]={{1,0},{0,0},{1,0}};

dcomplex X2[len]={{1,0},{0,0},{1,0},{0,0},{1,0}};

dcomplex X1_X2[len];

fft(X1,len,1);

fft(X2,len,1);

rep(i,len) X1_X2[i]=X1[i]*X2[i];

fft(X1_X2,len,-1);

Da(X1_X2)

}

Application:

输入两个整数,输出他们的乘积


#include “m.h”

using namespace FFT;

int main(){

string str1,str2

cin>>str1>>str2;

int len1= str1.size();

int len2= str2.size();

int len = (pow2(ceil(log2(max(len1,len2)))));

len*=2;//X1*X2 的最大长度

dcomplex X1[len],X2[len],X1_X2[len];

rep(i,len1) X1[i]=dcomplex(str1[len1-1-i]-‘0’,0);

rep(i,len2) X2[i]=dcomplex(str2[len2-1-i]-‘0’,0);

fft(X1,len,1);

fft(X2,len,1);

rep(i,len) X1_X2[i]=X1[i]*X2[i];

fft(X1_X2,len,-1);

Da(X1_X2)

dcomplex z(2,0);

long double ansapp=0;//这里为了方便直接加起来了

rep(i,len){

ansapp+=round(X1_X2[i].real())*pow(10,i);

}

cout<<setprecision(17)<<(ansapp);

}

 

Advertisements