FFT多项式乘法加速

 1 #include<iostream>
 2 #include<cstdio>
 3 #include<cmath>
 4 #include<algorithm>
 5 #include<queue>
 6 #include<cstring>
 7 #define PAU putchar(‘ ‘)
 8 #define ENT putchar(‘\n‘)
 9 using namespace std;
10 const double PI=acos(-1.0);
11 struct complex{
12     double r,i;
13     complex(double _r = 0.0,double _i = 0.0){r=_r;i=_i;}
14     complex operator +(const complex &b){return complex(r+b.r,i+b.i);}
15     complex operator -(const complex &b){return complex(r-b.r,i-b.i);}
16     complex operator *(const complex &b){return complex(r*b.r-i*b.i,r*b.i+i*b.r);}
17 };
18 void change(complex y[],int len){
19     int i,j,k;
20     for(i=1,j=len/2;i<len-1;i++){
21         if(i<j)swap(y[i],y[j]);
22         k=len/2;
23         while(j>=k) j-=k,k/=2;
24         if(j<k) j+=k;
25     } return;
26 }
27 void fft(complex y[],int len,int on){
28     change(y,len);
29     for(int h=2;h<=len;h<<=1){
30         complex wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
31         for(int j=0;j<len;j+=h){
32             complex w(1,0);
33             for(int k=j;k<j+h/2;k++){
34                 complex u=y[k],t=w*y[k+h/2];
35                 y[k]=u+t;
36                 y[k+h/2]=u-t;
37                 w=w*wn;
38             }
39         }
40     }
41     if(on==-1) for(int i=0;i<len;i++) y[i].r/=len;
42     return;
43 }
44 const int MAXN=2000+10;
45 complex x1[MAXN],x2[MAXN];
46 char str1[MAXN>>1],str2[MAXN>>1];
47 int sum[MAXN];
48 inline int read(){
49     int x=0,sig=1;char ch=getchar();
50     while(!isdigit(ch)){if(ch==‘-‘)sig=-1;ch=getchar();}
51     while(isdigit(ch))x=10*x+ch-‘0‘,ch=getchar();
52     return x*=sig;
53 }
54 inline void write(int x){
55     if(x==0){putchar(‘0‘);return;}if(x<0)putchar(‘-‘),x=-x;
56     int len=0,buf[15];while(x)buf[len++]=x%10,x/=10;
57     for(int i=len-1;i>=0;i--)putchar(buf[i]+‘0‘);return;
58 }
59 int len,len1,len2;
60 void init(){
61     scanf("%s%s",str1,str2);
62     len1=strlen(str1),len2=strlen(str2),len=1;
63     return;
64 }
65 void work(){
66     while(len<len1<<1||len<len2<<1) len<<=1;
67     for(int i=0;i<len1;i++) x1[i]=complex(str1[len1-1-i]-‘0‘,0);
68     for(int i=len1;i<len;i++) x1[i]=complex(0,0);
69     for(int i=0;i<len2;i++) x2[i]=complex(str2[len2-1-i]-‘0‘,0);
70     for(int i=len2;i<len;i++) x2[i]=complex(0,0);
71     fft(x1,len,1);fft(x2,len,1);
72     for(int i=0;i<len;i++) x1[i]=x1[i]*x2[i];
73     fft(x1,len,-1);
74     for(int i=0;i<len;i++) sum[i]=(int)(x1[i].r+0.5);
75     for(int i=0;i<len;i++){
76         sum[i+1]+=sum[i]/10;
77         sum[i]%=10;
78     } len=len1+len2-1;return;
79 }
80 void print(){
81     while(sum[len]<=0&&len>0) len--;
82     for(int i=len;i>=0;i--) putchar(sum[i]+‘0‘);
83     return;
84 }
85 int main(){init();work();print();return 0;}