快速傅里叶变换(FFT)模板

Posted 2020-12-28 mastervan

//有多少项N取多至少8~10倍 
#include <algorithm>
#include <cmath>
#include <complex>
#include <cstdlib>
using namespace std;
typedef complex<double> cn;
const int N=1e7+1;
const double Pi=3.1415926535;//Pi通常取10位即可（随意） 
cn a[N],b[N];
int r[N];
int mxb,bit;
int n,m;
//二进制翻转 
void Get_Rev(int bit) {
    for (int i=0;i<(1<<bit);i++)
    r[i]=(r[i>>1]>>1)|((i&1)<<(bit-1));
}
//傅里叶变换 
void DFT(cn *a,int mx,int inot) {//inot为1就是DFT，为-1就是IDFT 
    for (int i=0;i<mx;i++)
    if (i<r[i]) swap(a[i],a[r[i]]);
    for (int mlen=1;mlen<mx;mlen<<=1) {
        cn wn=exp(cn(0,inot*Pi/mlen));
        for (int len=mlen<<1,l=0;l<mx;l+=len) {
            cn wnk=cn(1,0);
            for (int k=l;k<l+mlen;k++,wnk*=wn) {
                cn x=a[k],y=wnk*a[k+mlen];
                a[k]=x+y;a[k+mlen]=x-y;//蝴蝶操作 
            }
        }
    }
    if (inot==-1)//IDFT的除n操作 
    for (int i=0;i<mx;i++) a[i]/=mx;
}

void FFT() {
    mxb=1;
    while (mxb<=n+m) mxb<<=1,bit++;
    Get_Rev(bit);
    DFT(a,mxb,1);DFT(b,mxb,1);
    for (int i=0;i<=mxb;i++) a[i]*=b[i];
    DFT(a,mxb,-1);
}

View Code