FFT快速傅里叶变换

FFT快速傅里叶变换

先放个模板咕一会儿

#include <bits/stdc++.h>
using namespace std;

//&Start
#define lng long long
#define lit long double
#define re register
#define kk(i,n) " \n"[i==n]
const int inf=0x3f3f3f3f;
const lng Inf=0x3f3f3f3f3f3f3f3f;
const lit Pi=acos(-1.0);

//&Complex
struct cp{
    double x,y;
    cp(double X=0,double Y=0){x=X,y=Y;}
    friend cp operator+(const cp a,const cp b){return cp(a.x+b.x,a.y+b.y);}
    friend cp operator-(const cp a,const cp b){return cp(a.x-b.x,a.y-b.y);}
    friend cp operator*(const cp a,const cp b){return cp(a.x*b.x-a.y*b.y,a.x*b.y+a.y*b.x);}
};

//&Data
const int N=4e6;
int n,m;
cp f[N+10];

//&FFT
int r[N+10],lim=1;
void FFT(cp a[],int t){
    for(re int i=1;i<=lim;i++)
        if(i<r[i]) swap(a[i],a[r[i]]);
    re cp rt,w,x,y;
    for(re int mid=1;mid<lim;mid<<=1){
        rt=cp(cos(Pi/mid),sin(Pi/mid)*t);
        for(re int j=0;j<lim;j+=(mid<<1)){
            w=cp(1,0);
            for(re int k=0;k<mid;k++){
                x=a[j|k];
                y=w*a[j|k|mid];
                a[j|k]=x+y;
                a[j|k|mid]=x-y;
                w=w*rt;
            }
        }
    }
    if(t==-1)
    for(re int i=0;i<=lim;i++)
        a[i].x/=lim,a[i].y/=lim;
}

//&Main
int main(){
    scanf("%d%d",&n,&m);
    for(int i=0;i<=n;i++) scanf("%lf",&f[i].x);
    for(int i=0;i<=m;i++) scanf("%lf",&f[i].y);
    int sm=n+m,h=-1;
    n=max(n,m);
    while(lim<=(n<<1)) lim<<=1,h++;
    for(re int i=1;i<=lim;i++)
        r[i]=(r[i>>1]>>1)|((i&1)<<h);
    FFT(f,1);
    for(re int i=0;i<=lim;i++) f[i]=f[i]*f[i];
    FFT(f,-1);
    for(re int i=0;i<=sm;i++)
        printf("%d%c",int(f[i].y/2+0.5),kk(i,sm));
    return 0;
}
上一篇:P4173 残缺的字符串 [FFT]


下一篇:使用Python 2.6快速分析音频文件的频谱?