FFT多项式乘法加速

FFT基本操作。。。讲解请自己看大学信号转置系列。。。

15-5-30更新：改成结构体的，跪烂王学长啊啊啊啊机智的static。。。

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<queue>
#include<cstring>
#define PAU putchar(' ')
#define ENT putchar('
')
#pragma comment(linker,"/STACK:10240000,10240000")
using namespace std;
const double PI=acos(-1.0);
const int maxn=800000+10;
struct FFT{
    struct cox{
        double r,i;cox(double _r = 0.0,double _i = 0.0){r=_r;i=_i;}
        cox operator +(const cox &b){return cox(r+b.r,i+b.i);}
        cox operator -(const cox &b){return cox(r-b.r,i-b.i);}
        cox operator *(const cox &b){return cox(r*b.r-i*b.i,r*b.i+i*b.r);}
    }f[maxn];int len,lenx;
    void init(int*s,int L,int Len){
        len=L;lenx=Len;
        for(int i=0;i<L;i++) f[i]=cox(s[L-1-i],0);
        for(int i=L;i<Len;i++) f[i]=cox(0,0);return;
    }
    void change(){
        for(int i=1,j=lenx>>1;i<lenx-1;i++){
            if(i<j)swap(f[i],f[j]);
            int k=lenx>>1;
            while(j>=k) j-=k,k>>=1;
            if(j<k) j+=k;
        } return;
    }
    void cal(int tp){
        change();
        for(int h=2;h<=lenx;h<<=1){
            double tr=-tp*2*PI/h;
            cox wn(cos(tr),sin(tr));
            for(int j=0;j<lenx;j+=h){
                cox w(1,0);
                for(int k=j;k<j+(h>>1);k++){
                    cox u=f[k],t=w*f[k+(h>>1)];
                    f[k]=u+t;f[k+(h>>1)]=u-t;w=w*wn;
                }
            }
        } if(tp==-1) for(int i=0;i<lenx;i++) f[i].r/=lenx;return;
    }
};
void mul(int*s1,int*s2,int L1,int L2,int&L,int*ans){
    L=1;while(L<L1<<1||L<L2<<1) L<<=1;
    static FFT a,b;a.init(s1,L1,L);b.init(s2,L2,L);a.cal(1);b.cal(1);
    for(int i=0;i<L;i++) a.f[i]=a.f[i]*b.f[i];a.cal(-1);
    for(int i=0;i<L;i++) ans[i]=(int){a.f[i].r+0.5};return;
}
int s1[maxn>>1],s2[maxn>>1],ans[maxn],L1,L2,L;
void init(){
    char ch;int tot;
    do ch=getchar(); while(!isdigit(ch));tot=0;
    while(isdigit(ch)){s1[tot++]=ch-'0';ch=getchar();}L1=tot;
    do ch=getchar(); while(!isdigit(ch));tot=0;
    while(isdigit(ch)){s2[tot++]=ch-'0';ch=getchar();}L2=tot;
    mul(s1,s2,L1,L2,L,ans);
    return;
}
void work(){
    for(int i=0;i<L;i++){
        ans[i+1]+=ans[i]/10;ans[i]%=10;
    } L=L1+L2-1;return;
}
void print(){
    while(ans[L]<=0&&L>0) L--;
    for(int i=L;i>=0;i--) putchar(ans[i]+'0');
    return;
}
int main(){init();work();print();return 0;}

原来写的丑哭了TAT

#include<iostream>
#include<cstdio>
#include<cmath>
#include<algorithm>
#include<queue>
#include<cstring>
#define PAU putchar(' ')
#define ENT putchar('
')
using namespace std;
const double PI=acos(-1.0);
struct complex{
    double r,i;
    complex(double _r = 0.0,double _i = 0.0){r=_r;i=_i;}
    complex operator +(const complex &b){return complex(r+b.r,i+b.i);}
    complex operator -(const complex &b){return complex(r-b.r,i-b.i);}
    complex operator *(const complex &b){return complex(r*b.r-i*b.i,r*b.i+i*b.r);}
};
void change(complex y[],int len){
    int i,j,k;
    for(i=1,j=len/2;i<len-1;i++){
        if(i<j)swap(y[i],y[j]);
        k=len/2;
        while(j>=k) j-=k,k/=2;
        if(j<k) j+=k;
    } return;
}
void fft(complex y[],int len,int on){
    change(y,len);
    for(int h=2;h<=len;h<<=1){
        complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j=0;j<len;j+=h){
            complex w(1,0);
            for(int k=j;k<j+h/2;k++){
                complex u=y[k],t=w*y[k+h/2];
                y[k]=u+t;
                y[k+h/2]=u-t;
                w=w*wn;
            }
        }
    }
    if(on==-1) for(int i=0;i<len;i++) y[i].r/=len;
    return;
}
const int MAXN=2000+10;
complex x1[MAXN],x2[MAXN];
char str1[MAXN>>1],str2[MAXN>>1];
int sum[MAXN];
inline int read(){
    int x=0,sig=1;char ch=getchar();
    while(!isdigit(ch)){if(ch=='-')sig=-1;ch=getchar();}
    while(isdigit(ch))x=10*x+ch-'0',ch=getchar();
    return x*=sig;
}
inline void write(int x){
    if(x==0){putchar('0');return;}if(x<0)putchar('-'),x=-x;
    int len=0,buf[15];while(x)buf[len++]=x%10,x/=10;
    for(int i=len-1;i>=0;i--)putchar(buf[i]+'0');return;
}
int len,len1,len2;
void init(){
    scanf("%s%s",str1,str2);
    len1=strlen(str1),len2=strlen(str2),len=1;
    return;
}
void work(){
    while(len<len1<<1||len<len2<<1) len<<=1;
    for(int i=0;i<len1;i++) x1[i]=complex(str1[len1-1-i]-'0',0);
    for(int i=len1;i<len;i++) x1[i]=complex(0,0);
    for(int i=0;i<len2;i++) x2[i]=complex(str2[len2-1-i]-'0',0);
    for(int i=len2;i<len;i++) x2[i]=complex(0,0);
    fft(x1,len,1);fft(x2,len,1);
    for(int i=0;i<len;i++) x1[i]=x1[i]*x2[i];
    fft(x1,len,-1);
    for(int i=0;i<len;i++) sum[i]=(int)(x1[i].r+0.5);
    for(int i=0;i<len;i++){
        sum[i+1]+=sum[i]/10;
        sum[i]%=10;
    } len=len1+len2-1;return;
}
void print(){
    while(sum[len]<=0&&len>0) len--;
    for(int i=len;i>=0;i--) putchar(sum[i]+'0');
    return;
}
int main(){init();work();print();return 0;}

搜索

复制