FFT快速傅立叶

 1 #include<cstdio>
 2 #include<iostream>
 3 #include<cmath>
 4 #include<cstring>
 5 #include<algorithm>
 6 #define maxn 131072
 7 #define pi 3.14159265358979323846
 8 using namespace std;
 9 char ch;
10 int m,n,ans[maxn];
11 bool ok;
12 struct comp{
13     double rea,ima;
14     void clear(){rea=ima=0;}
15     comp operator +(const comp &x){return (comp){rea+x.rea,ima+x.ima};}
16     comp operator -(const comp &x){return (comp){rea-x.rea,ima-x.ima};}
17     comp operator *(const comp &x){return (comp){rea*x.rea-ima*x.ima,rea*x.ima+ima*x.rea};}
18 }a[maxn],b[maxn],c[maxn],tmp[maxn],w,wn;
19 void read(int &x){
20     for (ok=0,ch=getchar();!isdigit(ch);ch=getchar()) if (ch==‘-‘) ok=1;
21     for (x=0;isdigit(ch);x=x*10+ch-‘0‘,ch=getchar());
22     if (ok) x=-x;
23 }
24 void read(comp *a){
25     for (int i=m-1;i>=0;i--){
26         for (ch=getchar();!isdigit(ch);ch=getchar());
27         a[i].rea=ch-‘0‘;
28     }
29 }
30 void write(int *a){
31     int i;
32     for (i=n-1;i>=0&&!a[i];i--);
33     if (i<0){puts("0");return;}
34     for (;i>=0;i--) printf("%d",(int)a[i]);
35     puts("");
36 }
37 void fft(comp *a,int st,int siz,int step,int op){
38     if (siz==1) return;
39     fft(a,st,siz>>1,step<<1,op),fft(a,st+step,siz>>1,step<<1,op);
40     int x=st,x1=st,x2=st+step;
41     w=(comp){1,0},wn=(comp){cos(op*2*pi/siz),sin(op*2*pi/siz)};
42     for (int i=0;i<(siz>>1);i++,x+=step,x1+=(step<<1),x2+=(step<<1),w=w*wn)
43         tmp[x]=a[x1]+(w*a[x2]),tmp[x+(siz>>1)*step]=a[x1]-(w*a[x2]);
44     for (int i=st;siz;i+=step,siz--) a[i]=tmp[i];
45 }
46 int main(){
47     read(m),n=1;
48     while (n<(m<<1)) n<<=1;
49     read(a),read(b);
50     fft(a,0,n,1,1),fft(b,0,n,1,1);
51     for (int i=0;i<n;i++) c[i]=a[i]*b[i];
52     fft(c,0,n,1,-1);
53     for (int i=0;i<n;i++) ans[i]=(int)round(c[i].rea/(1.0*n));
54     for (int i=0;i<n;i++) ans[i+1]+=ans[i]/10,ans[i]%=10;
55     write(ans);
56     system("pause");
57     return 0;
58 }