FFT（快速傅里叶变换）摘要

时间 2019-12-09

标签 fft 快速傅里叶变换摘要繁體版

原文原文链接

怎么说呢。。。此次的代码码风有点。。。html

从这篇博客开始，我终于会用LATEX了！撒花c++

注：如下涉及多项式的n均表示多项式的次数git

FFT用途

很简单，多项式相乘。ide

FFT原理

若是暴力乘，复杂度是$O(n^{2})$的，显然不行函数

因此考虑点值法。点值表示下的多项式相乘是$O(n)$的。这就是DFT（离散傅里叶变换）。ui

但很显然，随机带入n个值的复杂度是$O(n^{2})$的。spa

所以咱们考虑如下推导。翻译

设$$F(n)=a_{0}+a_{1}x+a_{2}x^{2}+...+a_{n}x^{n}$$那么咱们必定能够将其拆成两个函数，即$$G(n)=a_{0}+a_{2}x+a_{4}x^{2}+a_{6}x^{3}+...$$$$H(n)=a_{1}+a_{3}x+a_{5}x^{2}+a_{7}x^{3}+...$$因此$$F(n)=G(x^{2})+xH(x^{2})$$从这个式子咱们能够看出，若是咱们取这么一个点值集合，使得其中一半的值为全部值的平方，好比${-1,1}$，那么咱们就可使计算规模减半（对于$G(n)$和$H(n)$，咱们只需带入其中一半的值便可）。看来这是个很完美的策略。code

很不幸，若是咱们只在实数中这么思考，那复杂度仍为$O(n^{2})$，能够说没什么用。但这给了咱们一个启发，若是每次对函数进行拆分递归并使点值集合规模减半，咱们就能够在$O(nlogn)$的复杂度内求出一个多项式的点值表示。毫无疑问，问题的关键是这个点值集合。那是否存在这样一个点值集合呢？实数域固然没有，但复数域有，也就是单位根！（用复数作自变量真的让我很长时间没法接受）至于单位根么，出门右转问百度 htm

固然，既然用了FFT，天然要把多项式项数补成2的幂，具体作法为高次项补零。

递归式FFT的很好打，但常数大，且容易爆栈。所以咱们考虑迭代。

易发现，每次拆分多项式的时候，咱们老是按下标的奇偶性来拆分。这就致使了一个结果：原多项式通过递归拆分后的新多项式下标是原多项式下标的二进制反转。

 1 void unit() {  2     lg[1]=0;  3     fui(i,2,mxle,1)lg[i]=lg[i/2]+1;  4     bz[0]=1;  5     fui(i,1,mxlg,1)bz[i]=bz[i-1]<<1;  6     n=max(n,m);  7     n=bz[lg[n]+2]-1;  8     root[0].x=1;  9     root[1].x=cos((2*PI)/(n+1)); 10     root[1].y=sin((2*PI)/(n+1)); 11     fui(i,2,n,1)root[i]=root[i-1]*root[1];                                //求单位根
12     frot[0]=root[0]; 13     fui(i,1,n,1)frot[i]=root[n-i+1];                                    //求逆单位根（后面用
14     for(int i=0; i<=n; i++)fz[i]=(fz[i>>1]>>1)|((i&1)<<(lg[n+1]-1)) ;    //二进制反转 
15 }

预处理

而后就是FFT了

 1 void pre_FFT() {  2     fui(i,0,n,1)if(i<fz[i]) {  3         swp=fft[i];  4         fft[i]=fft[fz[i]];  5         fft[fz[i]]=swp;  6  }  7 }  8 void FFT() {  9     pre_FFT();                                                                        //对原多项式进行二进制反转
10     int zl=(n+1)>>1;                                                                //单位根增量
11     for(int i=1; i<=n; i=(i<<1)|1,zl>>=1) {                                            //每次处理的多项式长度
12         for(int j=0; j<n; j+=i+1) {                                                    //每次处理的多项式左端点
13             int mid=(j+j+i)>>1; 14             for(int k1=j,k2=mid+1,nw=0; k1<=mid; nw+=zl,k1++,k2++) {                //对该多项式的点值进行迭代处理
15                 x1=fft[k2]*root[nw]+fft[k1],x2=fft[k2]*root[nw+((n+1)>>1)]+fft[k1]; 16                 fft[k1]=x1; 17                 fft[k2]=x2; 18  } 19  } 20  } 21 }

FFT

记住，千万不要用stl里的复数，会被坑死的，强烈建议手打

IFFT（快速傅里叶逆变换）

FFT求出了点值，而后相乘，求出告终果多项式的点值表示。但这不是结果，咱们还要将其“翻译”回系数表示。因此就有了IFFT，能够叫其插值。

插值么，具体作法就是对结果求一遍相似FFT的东西，只不过把单位根改为{${{\omega}^{0},{\omega}^{-1},{\omega}^{-2},{\omega}^{-3},...,{\omega}^{-n}}$}，即上面预处理的时候求过的那个所谓逆单位根带入便可。证实么，我不会，反正做为一个OIer记结论就行了。具体证实参见自为风月马前卒大佬的博客。

 1 void IFFT() {  2  pre_FFT();  3     int zl=(n+1)>>1;  4     for(int i=1; i<=n; i=(i<<1)|1,zl>>=1) {  5         for(int j=0; j<n; j+=i+1) {  6             int mid=(j+j+i)>>1;  7             for(int k1=j,k2=mid+1,nw=0; k1<=mid; nw+=zl,k1++,k2++) {  8                 x1=fft[k2]*frot[nw]+fft[k1],x2=fft[k2]*frot[nw+((n+1)>>1)]+fft[k1];    //改为逆单位根带入
 9                 fft[k1]=x1; 10                 fft[k2]=x2; 11  } 12  } 13  } 14 }

IFFT

模板题

洛谷P3803 【模板】多项式乘法（FFT）

 1 #include<bits/stdc++.h>
 2 using namespace std;  3 #define INF 0x7fffffff
 4 #define ME 0x7f
 5 #define FO(s) freopen(s".in","r",stdin);freopen(s".out","w",stdout)
 6 #define fui(i,a,b,c) for(int i=(a);i<=(b);i+=(c))
 7 #define fdi(i,a,b,c) for(int i=(a);i>=(b);i-=(c))
 8 #define fel(i,a) for(register int i=h[a];i;i=b[i].ne)
 9 #define ll long long
10 #define ld double
11 #define MEM(a,b) memset(a,b,sizeof(a))
12 #define maxn 1000010
13 #define mxlg 21
14 #define mxle 4194310
15 const ld PI=3.1415926535897932384626; 16 int n,m,mn; 17 ld a0[mxle],b0[mxle]; 18 ld ans[mxle]; 19 int lg[mxle],bz[mxlg+3]; 20 int fz[mxle]; 21 struct Complex{ 22  ld x,y; 23  Complex(){} 24     Complex(ld __x,ld __y){x=__x;y=__y;} 25     Complex operator *(Complex xx)const{return (Complex){x*xx.x-y*xx.y,x*xx.y+y*xx.x};} 26     Complex operator *(ld xx)const{return (Complex){x*xx,y*xx};} 27     Complex operator +(Complex xx)const{return (Complex){x+xx.x,y+xx.y};} 28 }root[mxle],frot[mxle],a[mxle],b[mxle],fft[mxle],c[mxle],swp,x1,x2; 29 template<class T>
30 inline T read(T &n){ 31     n=0;int t=1;double x=10;char ch; 32     for(ch=getchar();!isdigit(ch)&&ch!='-';ch=getchar());(ch=='-')?t=-1:n=ch-'0'; 33     for(ch=getchar();isdigit(ch);ch=getchar()) n=n*10+ch-'0'; 34     if(ch=='.') for(ch=getchar();isdigit(ch);ch=getchar()) n+=(ch-'0')/x,x*=10; 35     return (n*=t); 36 } 37 template<class T>
38 T write(T n){ 39     if(n<0) putchar('-'),n=-n; 40     if(n>=10) write(n/10);putchar(n%10+'0');return n; 41 } 42 template<class T>
43 T writeln(T n){write(n);putchar('\n');return n;} 44 int inc(int &x){int in=(n+1)>>1;for(;x&in;in>>=1)x^=in;x|=in;return x;} 45 void unit(){ 46     lg[1]=0;fui(i,2,mxle,1)lg[i]=lg[i/2]+1; 47     bz[0]=1;fui(i,1,mxlg,1)bz[i]=bz[i-1]<<1; 48     n=max(n,m);n=bz[lg[n]+2]-1; 49     root[0].x=1; 50     root[1].x=cos((2*PI)/(n+1));root[1].y=sin((2*PI)/(n+1)); 51     fui(i,2,n,1)root[i]=root[i-1]*root[1]; 52     frot[0]=root[0]; 53     fui(i,1,n,1)frot[i]=root[n-i+1]; 54     fui(i,0,n,1)a[i]=root[0]*a0[i]; 55     fui(i,0,n,1)b[i]=root[0]*b0[i]; 56     for(int i=0;i<=n;i++)fz[i]=(fz[i>>1]>>1)|((i&1)<<(lg[n+1]-1)) ; 57 } 58 void pre_FFT(){ 59     fui(i,0,n,1)if(i<fz[i]){ 60         swp=fft[i];fft[i]=fft[fz[i]];fft[fz[i]]=swp; 61  } 62 } 63 void FFT(){ 64     pre_FFT();int zl=(n+1)>>1; 65     for(int i=1;i<=n;i=(i<<1)|1,zl>>=1){ 66         for(int j=0;j<n;j+=i+1){ 67             int mid=(j+j+i)>>1; 68             for(int k1=j,k2=mid+1,nw=0;k1<=mid;nw+=zl,k1++,k2++){ 69                 x1=fft[k2]*root[nw]+fft[k1],x2=fft[k2]*root[nw+((n+1)>>1)]+fft[k1]; 70                 fft[k1]=x1;fft[k2]=x2; 71  } 72  } 73  } 74 // pre_FFT();
75 } 76 void IFFT(){ 77     pre_FFT();int zl=(n+1)>>1; 78     for(int i=1;i<=n;i=(i<<1)|1,zl>>=1){ 79         for(int j=0;j<n;j+=i+1){ 80             int mid=(j+j+i)>>1; 81             for(int k1=j,k2=mid+1,nw=0;k1<=mid;nw+=zl,k1++,k2++){ 82                 x1=fft[k2]*frot[nw]+fft[k1],x2=fft[k2]*frot[nw+((n+1)>>1)]+fft[k1]; 83                 fft[k1]=x1;fft[k2]=x2; 84  } 85  } 86  } 87 } 88 int main(){ 89     read(n);read(m);mn=m+n; 90     fui(i,0,n,1)read(a0[i]);fui(i,0,m,1)read(b0[i]);unit(); 91     fui(i,0,n,1)fft[i]=a[i];FFT();fui(i,0,n,1)a[i]=fft[i]; 92     fui(i,0,n,1)fft[i]=b[i];FFT();fui(i,0,n,1)b[i]=fft[i]; 93     fui(i,0,n,1)c[i]=a[i]*b[i]; 94     fui(i,0,n,1)fft[i]=c[i];IFFT();fui(i,0,n,1)c[i]=fft[i]; 95     fui(i,0,mn,1)ans[i]=(c[i].x+0.5)/((n+1)); 96     fui(i,0,mn,1)printf("%.0f ",ans[i]); 97     return 0; 98 }

AC代码

常数太大，致使最后两个点2s多。。。