快速傅里叶变换

Posted 2020-11-19 mediocrekonjac

FFT

#include<bits/stdc++.h>
using namespace std;
#define N (1e7+10)
inline int read()
{
    char c=getchar();int x=0,f=1;
    while(c<‘0‘||c>‘9‘){if(c==‘-‘)f=-1;c=getchar();}
    while(c>=‘0‘&&c<=‘9‘){x=x*10+c-‘0‘;c=getchar();}
    return x*f;
}
const double Pi=acos(-1.0);
struct complex
{      
    double x,y;
    complex (double xx=0,double yy=0):x(xx),y(yy){}
};
complex operator + (complex a,complex b){ return complex(a.x+b.x , a.y+b.y);}
complex operator - (complex a,complex b){ return complex(a.x-b.x , a.y-b.y);}
complex operator * (complex a,complex b){ return complex(a.x*b.x-a.y*b.y , a.x*b.y+a.y*b.x);}
namespace fft
{    
    int l,r[N];
    int limit=1;
    void FFT(complex A[],int type)
    {
        for(RG int i=0;i<limit;++i) if(i<r[i]) swap(A[i],A[r[i]]);
        for(RG int mid=1;mid<limit;mid<<=1)
        {
            complex Wn(cos(Pi/mid),type*sin(Pi/mid));
            for(int R=mid<<1,j=0;j<limit;j+=R)
            {
                complex w(1,0);
                for(RG int k=0;k<mid;++k,w=w*Wn)
                {
                    complex x=A[j+k],y=w*A[j+k+mid];
                    A[j+k]=x+y;A[j+mid+k]=x-y;
                }
            }
        }
    }
    void multi(int & n,int & m,complex *a,complex *b)
    {
        limit=1,l=0;memset(r,0,sizeof(r))； 
        while(limit<=n+m) limit<<=1,++l;
        for(RG int i=0;i<limit;++i) r[i]=(r[i>>1]>>1)|((i&1)<<(l-1));
        FFT(a,1);FFT(b,1);
        for(RG int i=0;i<=limit;++i) a[i]=a[i]*b[i];
        FFT(a,-1);
        for(RG int i=0;i<=n+m;++i)    a[i].x=a[i].x/limit+0.5));
        n+=m;
    }
};

int main()
{
    
    return 0;
}