FFT各种模板

丑陋敬请谅解:

求两列数的卷积:

递归版:

#include <stdio.h>

#include <algorithm>

#include <math.h>

using namespace std;

const double pi=acos(-);

int p[],q[];

struct complex

{

    double x,y;

    complex(double xx=,double yy=)

    {

        x=xx;

        y=yy;

    }

}a[],b[];

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);}

void DFT(int limit,complex *a)

{

    if(limit==)

        return;

    int mid=limit>>;

    complex a1[mid],a2[mid];

    for(int i=;i<=limit;i+=)

    {

        a1[i>>]=a[i];

        a2[i>>]=a[i+];

    }

    DFT(mid,a1);

    DFT(mid,a2);

    complex wn=complex(cos(2.0*pi/limit),sin(2.0*pi/limit));

    complex w=complex(,);

    for(int i=;i<mid;i++,w=w*wn)

    {

        a[i]=a1[i]+w*a2[i];

        a[i+mid]=a1[i]-w*a2[i];

    }

}

int pow(int x,int y)

{

    if(y==)

        return ;

    int t=pow(x,y/);

    if(y%==)

        return t*t;

    return t*t*x;

}

void IDFT(int limit,complex *a)

{

    if(limit==)

        return;

    int mid=limit>>;

    complex a1[mid],a2[mid];

    for(int i=;i<=limit;i+=)

    {

        a1[i>>]=a[i];

        a2[i>>]=a[i+];

    }

    IDFT(mid,a1);

    IDFT(mid,a2);

    complex wn=complex(cos(2.0*pi/limit),-sin(2.0*pi/limit));

    complex w=complex(,);

    for(int i=;i<mid;i++,w=w*wn)

    {

        a[i]=a1[i]+w*a2[i];

        a[i+mid]=a1[i]-w*a2[i];

    }

}

int main()

{

    int k;

    scanf("%d",&k);

    int n=pow(,k);

    for(int i=;i<=n-;i++)

        scanf("%d",&p[i]);

    for(int i=;i<=n-;i++)

        scanf("%d",&q[i]);

    for(int i=;i<=n-;i++)

        a[i]=p[i];

    for(int i=;i<=n-;i++)

        b[i]=q[i];

    int limit=*n;

    DFT(limit,a);

    DFT(limit,b);

    for(int i=;i<=limit;i++)

        a[i]=a[i]*b[i];

    IDFT(limit,a);

    for(int i=;i<=n-;i++)

        printf("%d\n",(int)(a[i].x/limit+0.5));

    printf("%d",(int)(a[n-].x/limit+0.5));

}

非递归版+蝶形算法优化:

#include<stdio.h>

#include <algorithm>

#include<math.h>

using namespace std;

const int MAXN=1e7+;

const double Pi=acos(-1.0);

int p[MAXN],q[MAXN];

struct complex

{

        double x,y;

        complex (double xx=,double yy=){x=xx,y=yy;}

}a[MAXN],b[MAXN];

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);}

int N,M;

int l,r[MAXN];

int limit=;

void DFT(complex *A)

{

        for(int i=;i<limit;i++)

                if(i<r[i]) swap(A[i],A[r[i]]);

        for(int mid=;mid<limit;mid<<=)

        {

                complex Wn( cos(Pi/mid) , sin(Pi/mid) );

                for(int R=mid<<,j=;j<limit;j+=R)

                {

                        complex w(,);

                        for(int k=;k<mid;k++,w=w*Wn)

                        {

                                 complex x=A[j+k],y=w*A[j+mid+k];

                                A[j+k]=x+y;

                                A[j+mid+k]=x-y;

                        }

                }

        }

}

void IDFT(complex *A)

{

        for(int i=;i<limit;i++)

                if(i<r[i]) swap(A[i],A[r[i]]);

        for(int mid=;mid<limit;mid<<=)

        {

                complex Wn( cos(Pi/mid) , -1.0*sin(Pi/mid) );

                for(int R=mid<<,j=;j<limit;j+=R)

                {

                        complex w(,);

                        for(int k=;k<mid;k++,w=w*Wn)

                        {

                                 complex x=A[j+k],y=w*A[j+mid+k];

                                A[j+k]=x+y;

                                A[j+mid+k]=x-y;

                        }

                }

        }

}

int main()

{

        scanf("%d%d",&N,&M);

        for(int i=;i<=N;i++) scanf("%d",&p[i]);

        for(int i=;i<=M;i++) scanf("%d",&q[i]);

          for(int i=;i<=N;i++) a[i]=p[i];

        for(int i=;i<=M;i++)  b[i]=q[i];

        while(limit<=N+M) limit<<=,l++;

        for(int i=;i<limit;i++)

                r[i]= ( r[i>>]>> )| ( (i&)<<(l-) ) ;

        DFT(a);

        DFT(b);

        for(int i=;i<=limit;i++) a[i]=a[i]*b[i];

        IDFT(a);

        for(int i=;i<=N+M;i++)

                printf("%d ",(int)(a[i].x/limit+0.5));

        return ;

}

FFT版高精度乘法:

#include<stdio.h>

#include <algorithm>

#include<math.h>

#include <string.h>

using namespace std;

const double Pi=acos(-1.0);

char s1[],s2[];

int sum[];

struct complex

{

    double x,y;

    complex (double xx=,double yy=){x=xx,y=yy;}

}a[],b[];

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);}

int N,M,L;

int len=;

int l,r[];

int limit=;

void  fft(complex *A,int limit,int type)

{

    for(int i=;i<limit;i++)

        if(i<r[i]) swap(A[i],A[r[i]]);

    for(int mid=;mid<limit;mid<<=)

    {

        complex Wn( cos(Pi/mid) , type*sin(Pi/mid) );

        for(int R=mid<<,j=;j<limit;j+=R)

        {

            complex w(,);

            for(int k=;k<mid;k++,w=w*Wn)

            {

                 complex x=A[j+k],y=w*A[j+mid+k];

                A[j+k]=x+y;

                A[j+mid+k]=x-y;

            }

        }

    }

    if(type==-)for(int i=;i<limit;i++)A[i].x/=limit;

}

int main()

{

    int n;

    scanf("%d",&n);

    scanf("%s",s1);

    scanf("%s",s2);

    int len1=strlen(s1);

    int len2=strlen(s2);

        int lenx=len1+len2;

        for(len=;len<lenx;len<<=)

            L++;

        for(int i=;i<len;i++)

            r[i]=((r[i>>])>>)|((i&)<<(L-));

        for(int i=;i<len1;i++)

            a[i]=complex(s1[len1-i-]-'',);

        for(int i=len1;i<len;i++)

            a[i]=complex(,);

        for(int i=;i<len2;i++)

            b[i]=complex(s2[len2-i-]-'',);

        for(int i=len2;i<len;i++)

            b[i]=complex(,);

         fft(a,len,);fft(b,len,);

        for(int i=;i<len;i++)

            a[i]=a[i]*b[i];

        fft(a,len,-);

        for(int i=;i<len;i++)

            sum[i]=int(a[i].x+0.5);

        for(int i=;i<len;i++)

        {

            sum[i+]+=sum[i]/;

            sum[i]%=;

        }

    len=len1+len2-;

        while(sum[len]<= && len>)

            len--;

        for(int i=len;i>=;i--)

            printf("%c",sum[i]+'');

    return ;

}

巴特西

FFT各种模板

最新文章

热门文章