1 /*  2     algorithm : High-Precision FFT  3   4 */  5 #include <cstdio>  6 #include <cstring>  7 #include <cmath>  8 #include <algorithm>  9 #define N 200005 10 #define pi acos(-1.0) // PI值 11 using namespace std; 12 struct complex 13 { 14     double r,i; 15     complex(double real=0.0,double image=0.0){ 16         r=real; i=image; 17     } 18     // 以下为三种虚数运算的定义 19     complex operator + (const complex o){ 20         return complex(r+o.r,i+o.i); 21     } 22     complex operator - (const complex o){ 23         return complex(r-o.r,i-o.i); 24     } 25     complex operator * (const complex o){ 26         return complex(r*o.r-i*o.i,r*o.i+i*o.r); 27     } 28 }x1[N],x2[N]; 29 char a[N/2],b[N/2]; 30 int sum[N]; // 结果存在sum里 31 void brc(complex *y,int l) // 二进制平摊反转置换 O(logn) 32 { 33     register int i,j,k; 34     for(i=1,j=l/2;i<l-1;i++) 35     { 36         if(i<j)  swap(y[i],y[j]); // 交换互为下标反转的元素 37                                 // i<j保证只交换一次 38         k=l/2; 39         while(j>=k) // 由最高位检索，遇1变0，遇0变1，跳出 40         { 41             j-=k; 42             k/=2; 43         } 44         if(j<k)  j+=k; 45     } 46 } 47 void fft(complex *y,int l,double on) // FFT O(nlogn) 48                             // 其中on==1时为DFT，on==-1为IDFT 49 { 50     register int h,i,j,k; 51     complex u,t; 52     brc(y,l); // 调用反转置换 53     for(h=2;h<=l;h<<=1) // 控制层数 54     { 55         // 初始化单位复根 56         complex wn(cos(on*2*pi/h),sin(on*2*pi/h)); 57         for(j=0;j<l;j+=h) // 控制起始下标 58         { 59             complex w(1,0); // 初始化螺旋因子 60             for(k=j;k<j+h/2;k++) // 配对 61             { 62                 u=y[k]; 63                 t=w*y[k+h/2]; 64                 y[k]=u+t; 65                 y[k+h/2]=u-t; 66                 w=w*wn; // 更新螺旋因子 67             } // 据说上面的操作叫蝴蝶操作… 68         } 69     } 70     if(on==-1)  for(i=0;i<l;i++) y[i].r/=l; // IDFT 71 } 72 int main(void){ 73     int l1,l2,l; 74     register int i; 75     while(scanf("%s%s",a,b)!=EOF){ 76         l1 = strlen(a),l2 = strlen(b); 77         l = 1; while(l < l1 * 2 || l < l2 * 2) l <<= 1; // 将次数界变成2^n 配合二分与反转置换 78         for(i = 0 ; i < l1 ; i++){ // 倒置存入 79             x1[i].r = a[l1 - i - 1] - ‘0‘; 80             x1[i].i = 0.0; 81         } 82         for( ; i < l ; i++) x1[i].r = x1[i].i = 0.0; 83         // 将多余次数界初始化为0 84         for(i = 0 ; i < l2 ; i++){ // same 85             x2[i].r = b[l2 - i - 1] - ‘0‘; 86             x2[i].i = 0.0; 87         } 88         for( ; i < l ; i++) x2[i].r = x2[i].i = 0.0; 89         fft(x1,l,1); // DFT(a) 90         fft(x2,l,1); // DFT(b) 91         for(i = 0 ; i < l ; i++) x1[i] = x1[i] * x2[i]; // 点乘结果存入a 92         fft(x1,l,-1); // IDFT(a*b) 93         for(i = 0 ; i < l ; i++) sum[i] = x1[i].r + 0.5; // 四舍五入 94         for(i = 0 ; i < l ; i++){ // 进位 95             sum[i + 1] += sum[i] / 10; 96             sum[i] %= 10; 97         } 98         l = l1 + l2 - 1; 99         while(sum[l] <= 0 && l > 0) l--; // 检索最高位100         for(i = l ; i >= 0 ; i--) putchar(sum[i] + ‘0‘); // 倒序输出101         putchar(‘\n‘);102     }103     return 0;104 }