FFT详细推导

FFT(傅里叶快速变换)一.前置知识1.复数和单位根2.单位根的三个引理3.多项式二.FFT(快速傅里叶变换推导)三.IFFT四.FFT求解多项式乘积模板代码1.递归版2.非递归版(这个更快，省去了递归时间)五.视频资源六.FFT题目集

FFT(傅里叶快速变换)

FFT在实际工程中有着非常的广泛，尤其是在信号领域，在ACM算法竞赛领域主要可以用来快速计算多项式的乘积

一.前置知识

1.复数和单位根

有人会觉得FFT怎么会用到复数的知识，想必是有点古怪，后面在推导过程中会介绍到它的用处。

2.单位根的三个引理

3.多项式

前置知识学完了之后下面我们来推导FFT

二.FFT(快速傅里叶变换推导)

DFT的作用是把多项式的系数表示法，转化为点值表示法,复杂度O(n2)O(n^2)O(n2)，而FFT则是快速版的DFT，作用是一样的，FFT的复杂度是O(nlogn)O(nlogn)O(nlogn).

上面我们说过n次多项式需要n+1个点来唯一确定，那么我们找点的时候为什么带入的x是复数(也即wi)而不是实数呢？这个视频讲的很详细FFT推导过程。

FFT模板：

/*这里的opt=1*/void FFT(Complex *a,int lim,int opt){if(lim==1) return ;Complex a0[lim>>1],a1[lim>>1];/*我们把多项式的奇数项和偶数项拆开*/ for(int i=0;i<lim;i+=2){a0[i>>1]=a[i],a1[i>>1]=a[i+1];}FFT(a0,lim>>1,opt);//然后递归拆解奇数项 FFT(a1,lim>>1,opt);//递归拆解偶数项 Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));//单位根Complex w = Complex(1,0);//第一个根w0 for(int k=0;k<(lim>>1);k++){Complex t=w*a1[k];//这样会少一次复数运算； a[k]=a0[k]+t;//最后我们把求得的值代回 a[k+(lim>>1)]=a0[k]-t;w=w*wn;//想当于次幂+1 }return ; }

三.IFFT

我们要解决多项式a和多项式b的乘积运算，FFT把多项式转化为点值表示之后，我们对两个多项式进行乘积运算，得到一个新的多项式C的点值表示，我们还需要把C的点值表示转化为系数表示，这样才能得到每一项的系数，这时就需要用到FFT的逆过程，IFFT(快速傅里叶逆变换)，也就是IDFT(离散傅里叶逆变换)的快速版。

IFFT模板：

/*这里的opt=-1*/void FFT(Complex *a,int lim,int opt){if(lim==1) return ;Complex a0[lim>>1],a1[lim>>1];/*我们把多项式的奇数项和偶数项拆开*/ for(int i=0;i<lim;i+=2){a0[i>>1]=a[i],a1[i>>1]=a[i+1];}FFT(a0,lim>>1,opt);//然后递归拆解奇数项 FFT(a1,lim>>1,opt);//递归拆解偶数项 Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));//单位根Complex w = Complex(1,0);//第一个根w0 for(int k=0;k<(lim>>1);k++){Complex t=w*a1[k];//这样会少一次复数运算； a[k]=a0[k]+t;//最后我们把求得的值代回 a[k+(lim>>1)]=a0[k]-t;w=w*wn;//想当于次幂+1 }return ; }

通过FFT和IFFT的对比我们发现只有opt参数不一样其余地方全是相同的，FFT的opt=1，IFFT的opt=-1，而opt就只作用于下面这一段代码块。

Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));

这就是我们上面这张图的解释

最后a/n，因为a是一个Complex，所以a/n指的是a的实部除以n，也即a.r/n，因为我们的系数是存在Complex的实部里面的，虚部只是用于FFT和IFFT的递归运算过程。

四.FFT求解多项式乘积模板代码

1.递归版

下面给出整个代码：

#include<bits/stdc++.h>using namespace std;const double PI=acos(-1);//3.1415926... const int Maxn=1e5;struct Complex{double r,i;Complex() {r=0.0,i=0.0;}//不带参数的构造函数 Complex(double real,double imag):r(real),i(imag){}//带参数的构造函数 }F[Maxn],G[Maxn];Complex operator + (Complex a,Complex b){return Complex(a.r+b.r,a.i+b.i);}Complex operator - (Complex a,Complex b){return Complex(a.r-b.r,a.i-b.i);}Complex operator * (Complex a,Complex b){return Complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}/*快速傅里叶变换*/void FFT(Complex *a,int lim,int opt){if(lim==1) return ;Complex a0[lim>>1],a1[lim>>1];/*我们把多项式的奇数项和偶数项拆开*/ for(int i=0;i<lim;i+=2){a0[i>>1]=a[i],a1[i>>1]=a[i+1];}FFT(a0,lim>>1,opt);//然后递归拆解奇数项 FFT(a1,lim>>1,opt);//递归拆解偶数项 Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));//单位根Complex w = Complex(1,0);//第一个根w0 for(int k=0;k<(lim>>1);k++){Complex t=w*a1[k];//这样会少一次复数运算； a[k]=a0[k]+t;//最后我们把求得的值代回 a[k+(lim>>1)]=a0[k]-t;w=w*wn;//想当于次幂+1 }return ; } int main(){int n,m,lim=1;scanf("%d %d",&n,&m);for(int i=0;i<=n;i++) scanf("%lf",&F[i].r);for(int i=0;i<=m;i++) scanf("%lf",&G[i].r);int len = n+m;while(lim<=len) lim<<=1;// FFT(F,lim,1);//我们先通过FFT，把a多项式的系数表示转化为点对表示 FFT(G,lim,1);//我们先通过FFT，也把b多项式的系数表示转化为点对表示/*这时候我们对a，b进行卷积操作，就算F,G没有lim那么多项，但是我们给复数默认传参是0+0i，所以卷积对结果没有影响 */ for(int i=0;i<=lim;i++) F[i]=F[i]*G[i]; FFT(F,lim,-1);//然后我们通过IFFT，也就是快速傅里叶逆变换，把点对表示转化为系数表示 for(int i=0;i<=n+m;i++) printf("%d ",(int)(F[i].r/lim+0.5));//最后的每一个系数/n,然后四舍五入就是答案 }

我们这里测试一下

第一行 我们输入n，m分别是3，3表示多项式a，多项式b的最高次次数。

第二行 我们输入多项式a各项的系数(默认从低次到高次项)

1 2 3 4

第三行 我们输入多项式b各项的系数(默认从低次到高次项)

1 2 3 4

我们的结果是

1 4 10 20 25 24 16

上面想当于解决了多项式A×多项式B的问题，也即A⨁BA\bigoplus BA⨁B

其中A(x)=1+2x+3x2+4x3A(x) = 1+2x+3x^2+4x^3A(x)=1+2x+3x2+4x3,B(x)=1+2x+3x2+4x3B(x) = 1+2x+3x^2+4x^3B(x)=1+2x+3x2+4x3

计算得到C(x)=1+4x+10x2+20x3+25x4+24x5+16x6C(x) = 1+4x+10x^2+20x^3+25x^4+24x^5+16x^6C(x)=1+4x+10x2+20x3+25x4+24x5+16x6

2.非递归版(这个更快，省去了递归时间)

非递归版需要用到蝴蝶效应的知识，本蒟蒻还没有学习完，待更…

五.视频资源

如果看了笔记还是不明白的话，这里推荐一个非常详细的视频讲解，看不懂直接来打我！

建议两个都仔细看一遍，有必要做一下笔记尝试推导

1.FFT的推导细节说明

2.全网最详细FFT,DFT,IDFT，IFFT讲解。

六.FFT题目集

1.模板题，这道题虽然可以之前可以python，现在数据加强了不知道题目建议用FFT，注意高精度进位就可以了，还有就是，数组开到5e7不然会re

P1919 【模板】A*B Problem升级版（FFT快速傅里叶）

AC代码：

#include<bits/stdc++.h>using namespace std;const double PI=acos(-1);//3.1415926... const int Maxn=5e6+5;struct Complex{double r,i;Complex() {r=0.0,i=0.0;}//不带参数的构造函数 Complex(double real,double imag):r(real),i(imag){}//带参数的构造函数 }F[Maxn],G[Maxn];Complex operator + (Complex a,Complex b){return Complex(a.r+b.r,a.i+b.i);}Complex operator - (Complex a,Complex b){return Complex(a.r-b.r,a.i-b.i);}Complex operator * (Complex a,Complex b){return Complex(a.r*b.r-a.i*b.i,a.r*b.i+a.i*b.r);}/*快速傅里叶变换*/void FFT(Complex *a,int lim,int opt){if(lim==1) return ;Complex a0[lim>>1],a1[lim>>1];/*我们把多项式的奇数项和偶数项拆开*/ for(int i=0;i<lim;i+=2){a0[i>>1]=a[i],a1[i>>1]=a[i+1];}FFT(a0,lim>>1,opt);//然后递归拆解奇数项 FFT(a1,lim>>1,opt);//递归拆解偶数项 Complex wn = Complex(cos(2.0*PI/lim),opt*sin(2.0*PI/lim));//单位根Complex w = Complex(1,0);//第一个根w0 for(int k=0;k<(lim>>1);k++){Complex t=w*a1[k];//这样会少一次复数运算； a[k]=a0[k]+t;//最后我们把求得的值代回 a[k+(lim>>1)]=a0[k]-t;w=w*wn;//想当于次幂+1 }return ; } char a[Maxn],b[Maxn]; int ans[Maxn];int main(){int n=0,m=0,lim=1;scanf("%s",a);scanf("%s",b);n=strlen(a)-1;m=strlen(b)-1;for(int i=n;i>=0;i--) F[n-i].r = a[i]-'0';for(int i=m;i>=0;i--) G[m-i].r = b[i]-'0';int len = n+m;while(lim<=len) lim<<=1;FFT(F,lim,1);//我们先通过FFT，把a多项式的系数表示转化为点对表示 FFT(G,lim,1);//我们先通过FFT，也把b多项式的系数表示转化为点对表示for(int i=0;i<=lim;i++) F[i]=F[i]*G[i]; FFT(F,lim,-1);//然后我们通过IFFT，也就是快速傅里叶逆变换，把点对表示转化为系数表示 for(int i=0;i<=n+m;i++) ans[i]=(int)(F[i].r/lim+0.5);//最后的每一个系数/n,然后四舍五入就是答案 int num=0;for(int i=0;i<=n+m;i++){ans[i]+=num;num=ans[i]/10;ans[i]%=10;}while(num){ans[++len] = num%10;num/=10;}for(int i=len;i>=0;i--) cout<<ans[i];return 0;}