快速傅里叶变换 FFT

  1 #include<iostream>
  2 #include<cstdio>
  3 #include<cmath>
  4 #include<cstring>
  5 #include<vector>
  6 #include<algorithm>
  7 using namespace std;
  8 #define MAXN 550000+100
  9 #define PI 3.14159265358
 10 #define EPS 1e-2
 11 char a[MAXN],b[MAXN];
 12 int n,m,top=0,label=0;
 13 int ans[MAXN];
 14 struct complex
 15 {
 16     double R,I;
 17     complex(){}
 18     complex(double r0,double i0): R(r0),I(i0) {}
 19 };
 20 complex operator *(const complex& A,const complex &B)
 21 {
 22     return complex(A.R*B.R-A.I*B.I,A.R*B.I+B.R*A.I);
 23 }
 24 complex operator +(const complex& A,const complex &B)
 25 {
 26     return complex(A.R+B.R,A.I+B.I);
 27 }
 28 complex operator -(const complex& A,const complex &B)
 29 {
 30     return complex(A.R-B.R,A.I-B.I);
 31 }
 32 vector<complex> C[3*MAXN];
 33 complex p[MAXN];
 34 int FFT(int first,int floor,int multi)
 35 {
 36     int i;
 37     if(n==floor) 
 38     {
 39         C[++top].push_back(p[first]);
 40         return top;
 41     }
 42     complex wn=complex(cos(2*PI/n*floor),sin(2*PI/n*floor)*multi);
 43     complex w=complex(1,0);
 44     label++;
 45     int top0=FFT(first,floor*2,multi);
 46     int top1=FFT(first+floor,floor*2,multi);
 47     top++;
 48     for(i=0;i<n/floor/2;i++,w=w*wn)
 49         C[top].push_back(C[top0][i]+w*C[top1][i]);
 50 
 51     w=complex(1,0);
 52     for(i=0;i<n/floor/2;i++,w=w*wn)
 53         C[top].push_back(C[top0][i]-w*C[top1][i]);
 54     return top;
 55 }
 56 int make_int(double x)
 57 {
 58     int A=x;
 59     if (fabs(x-A)<EPS)
 60         return A;
 61     return A+1;
 62 }
 63 void solve()
 64 {
 65     int i=1;
 66     for(i=0;i<n;i++)
 67         p[i].R=a[n-i-1]-48;
 68     while(i<n) i*=2;
 69     n=i*2;
 70     
 71     i=1;
 72     while(i<m) i*=2;
 73     n=max(i*2,n);
 74 
 75     int top0=FFT(0,1,1);
 76 
 77     memset(p,0,sizeof(p));
 78     for(i=0;i<m;i++)
 79         p[i].R=b[m-i-1]-48;
 80     int top1=FFT(0,1,1);
 81 
 82     memset(p,0,sizeof(p));
 83     for(i=0;i<n;i++)
 84         p[i]=C[top0][i]*C[top1][i];
 85     int top2=FFT(0,1,-1);
 86 
 87     for(i=0;i<n;i++) 
 88     {
 89         ans[i]+=make_int(C[top2][i].R/n);
 90         ans[i+1]=ans[i]/10;
 91         ans[i]=ans[i]%10;
 92     }
 93 
 94     for(i=n-1;i>=0;i--) if(ans[i]!=0) break;
 95     for(i;i>=0;i--) printf("%d",ans[i]);
 96     printf("\n");
 97 }
 98 
 99 int main()
100 {
101     int i;
102     scanf("%s%s",a,b);
103     n=strlen(a); m=strlen(b);
104     solve();
105 }
106