FFT各种模板

丑陋敬请谅解:

求两列数的卷积:

递归版:

#include <stdio.h>
#include <algorithm>
#include <math.h>
using namespace std;
const double pi=acos(-1);
int p[5000010],q[5000010];
struct complex
{
    double x,y;
    complex(double xx=0,double yy=0)
    {
        x=xx;
        y=yy;
    }
}a[5000100],b[5000100];
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==1)
        return;
    int mid=limit>>1;
    complex a1[mid],a2[mid];
    for(int i=0;i<=limit;i+=2)
    {
        a1[i>>1]=a[i];
        a2[i>>1]=a[i+1];
    }
    DFT(mid,a1);
    DFT(mid,a2);
    complex wn=complex(cos(2.0*pi/limit),sin(2.0*pi/limit));
    complex w=complex(1,0);
    for(int i=0;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==0)
        return 1;
    int t=pow(x,y/2);
    if(y%2==0)
        return t*t;
    return t*t*x;
}
void IDFT(int limit,complex *a)
{
    if(limit==1)
        return;
    int mid=limit>>1;
    complex a1[mid],a2[mid];
    for(int i=0;i<=limit;i+=2)
    {
        a1[i>>1]=a[i];
        a2[i>>1]=a[i+1];
    }
    IDFT(mid,a1);
    IDFT(mid,a2);
    complex wn=complex(cos(2.0*pi/limit),-sin(2.0*pi/limit));
    complex w=complex(1,0);
    for(int i=0;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(2,k);
    for(int i=0;i<=n-1;i++)
        scanf("%d",&p[i]);
    for(int i=0;i<=n-1;i++)
        scanf("%d",&q[i]);
    for(int i=0;i<=n-1;i++)
        a[i]=p[i];
    for(int i=0;i<=n-1;i++)
        b[i]=q[i];
    int limit=2*n;
    DFT(limit,a);
    DFT(limit,b);
    for(int i=0;i<=limit;i++)
        a[i]=a[i]*b[i];
    IDFT(limit,a);
    for(int i=0;i<=n-2;i++)
        printf("%d
",(int)(a[i].x/limit+0.5));
    printf("%d",(int)(a[n-1].x/limit+0.5));
}

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

#include<stdio.h>
#include <algorithm>
#include<math.h>
using namespace std;
const int MAXN=1e7+10;
const double Pi=acos(-1.0);
int p[MAXN],q[MAXN];
struct complex
{
        double x,y;
        complex (double xx=0,double yy=0){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=1;
void DFT(complex *A)
{
        for(int i=0;i<limit;i++) 
                if(i<r[i]) swap(A[i],A[r[i]]);
        for(int mid=1;mid<limit;mid<<=1)
        {
                complex Wn( cos(Pi/mid) , sin(Pi/mid) );
                for(int R=mid<<1,j=0;j<limit;j+=R)
                {
                        complex w(1,0);
                        for(int k=0;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=0;i<limit;i++) 
                if(i<r[i]) swap(A[i],A[r[i]]);
        for(int mid=1;mid<limit;mid<<=1)
        {
                complex Wn( cos(Pi/mid) , -1.0*sin(Pi/mid) );
                for(int R=mid<<1,j=0;j<limit;j+=R)
                {
                        complex w(1,0);
                        for(int k=0;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=0;i<=N;i++) scanf("%d",&p[i]);
        for(int i=0;i<=M;i++) scanf("%d",&q[i]);
          for(int i=0;i<=N;i++) a[i]=p[i];
        for(int i=0;i<=M;i++)  b[i]=q[i];
        while(limit<=N+M) limit<<=1,l++;
        for(int i=0;i<limit;i++)
                r[i]= ( r[i>>1]>>1 )| ( (i&1)<<(l-1) ) ;
        DFT(a);
        DFT(b);
        for(int i=0;i<=limit;i++) a[i]=a[i]*b[i];
        IDFT(a);
        for(int i=0;i<=N+M;i++)
                printf("%d ",(int)(a[i].x/limit+0.5));
        return 0;
}

FFT版高精度乘法:

#include<stdio.h>
#include <algorithm>
#include<math.h>
#include <string.h>
using namespace std;
const double Pi=acos(-1.0);
char s1[1200010],s2[1200010];
int sum[1200010];
struct complex
{
    double x,y;
    complex (double xx=0,double yy=0){x=xx,y=yy;}
}a[1200010],b[1200010];
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=1;
int l,r[1200010];
int limit=1;
void  fft(complex *A,int limit,int type)
{
    for(int i=0;i<limit;i++)
        if(i<r[i]) swap(A[i],A[r[i]]);
    for(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(int k=0;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==-1)for(int i=0;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=1;len<lenx;len<<=1)
            L++;
        for(int i=0;i<len;i++)
            r[i]=((r[i>>1])>>1)|((i&1)<<(L-1));
        for(int i=0;i<len1;i++)
            a[i]=complex(s1[len1-i-1]-'0',0);
        for(int i=len1;i<len;i++)
            a[i]=complex(0,0);
        for(int i=0;i<len2;i++)
            b[i]=complex(s2[len2-i-1]-'0',0);
        for(int i=len2;i<len;i++)
            b[i]=complex(0,0);
         fft(a,len,1);fft(b,len,1);
        for(int i=0;i<len;i++)
            a[i]=a[i]*b[i];
        fft(a,len,-1);
        for(int i=0;i<len;i++)
            sum[i]=int(a[i].x+0.5);
        for(int i=0;i<len;i++)
        {
            sum[i+1]+=sum[i]/10;
            sum[i]%=10;
        }
    len=len1+len2-1;
        while(sum[len]<=0 && len>0)
            len--;
        for(int i=len;i>=0;i--)
            printf("%c",sum[i]+'0');
    return 0;
}