Time Limit: 2000/1000 MS (Java/Others) Memory Limit: 65536/32768 K (Java/Others)
Total Submission(s): 9413 Accepted Submission(s): 1468
Problem Description
Calculate A * B.
Input
Each line will contain two integers A and B. Process to end of file.
Note: the length of each integer will not exceed 50000.
Note: the length of each integer will not exceed 50000.
Output
For each case, output A * B in one line.
Sample Input
1 2 1000 2
Sample Output
2 2000
Author
DOOM III
Recommend
DOOM III
就一个高精度乘法 FFT加速。
最近正好要捡起fft,就顺便整理了模板。
FFT的原理还是算法导论靠谱,没有那么艰深难懂,就涉及怎么进行FFT和FFT需要的原理和定理。
看看算法导论里FFT的部分,一定要读到迭代实现那部分!!
看了好久求和引理,才发觉他是为了保证$w_n^k$与$w_n^{k+2/h}$的对称性(即$w_n^{k+2/h}=-w_n^k$)的,这个引理是必要的。
对于多项式序列,我们可以用两个O(nlgn)(n>max(len1,len2)*2)的FFT将其系数表示转化为点值表示(DFT),然后用O(n) 相乘,接着用FFT把结果的点值表示变为系数表示(IDFT),总体算起来是3O(nlgn)+O(n),即O(nlgn)的时间复杂度。比O(n^2)好多了。
以下是学习的两个版本。
1 #include<bits/stdc++.h> 2 #define clr(x) memset(x,0,sizeof(x)) 3 #define clr_1(x) memset(x,-1,sizeof(x)) 4 #define clrmax(x) memset(x,0x3f3f3f3f,sizeof(x)) 5 #define LL long long 6 #define mod 1000000007 7 #define PI 3.1415926535 8 using namespace std; 9 char s1[200010],s2[200010]; 10 int a[200010],b[200010]; 11 //复数序列结构体 12 struct complexed 13 { 14 double r,i; 15 complexed(double _r=0.0,double _i=0.0) 16 { 17 r=_r; 18 i=_i; 19 } 20 complexed operator +(complexed b) 21 { 22 return complexed(r+b.r,i+b.i); 23 } 24 complexed operator -(complexed b) 25 { 26 return complexed(r-b.r,i-b.i); 27 } 28 complexed operator *(complexed b) 29 { 30 return complexed(r*b.r-i*b.i,i*b.r+r*b.i); 31 } 32 }num1,num2; 33 vector<complexed> multi1,multi2; 34 inline int max(int a,int b) 35 { 36 return a>b?a:b; 37 } 38 //并将长度变为2…^(k+1) 39 void changelen(int &len) 40 { 41 int mul=1; 42 while(mul<len) 43 mul<<=1; 44 mul<<=1; 45 len=mul; 46 return ; 47 } 48 //将整数序列复制到复数序列中 49 void copyed(int *a,vector<complexed> &multi,int len) 50 { 51 multi.resize(len); 52 for(int i=0;i<len;i++) 53 multi[i]=(complexed){a[i],0}; 54 return; 55 } 56 //DFT的话on=1,IDFT on=-1; 57 void fft(vector<complexed> &multi,int len,int on) 58 { 59 complexed wn,w,u,t; 60 //wn,w,u,t如算法导论中所示 61 vector<complexed> ans; 62 ans.resize(len); 63 //ans存每次操作计算后的y,最后再作为下次的multi。 64 for(int h=len/2;h>=1;h>>=1) 65 { 66 wn=(complexed){cos(2*on*PI/(len/h)),sin(2*on*PI/(len/h))}; 67 for(int i=0;i<h;i++) 68 { 69 w=(complexed){1,0}; 70 for(int j=0;j<len/h/2;j++) 71 { 72 //蝴蝶操作 73 u=multi[i+2*h*j]; 74 t=multi[i+2*h*j+h]*w; 75 ans[i+h*j]=u+t; 76 ans[i+h*j+len/2]=u-t; 77 w=w*wn; 78 } 79 } 80 //ans作为下次计算的multi 81 multi=ans; 82 } 83 //IDFT每个元素都得除以n 84 if(on==-1) 85 for(int i=0;i<len;i++) 86 multi[i].r/=len; 87 return ; 88 } 89 int main() 90 { 91 int len1,len2,len; 92 while(scanf("%s%s",s1,s2)!=EOF) 93 { 94 len1=strlen(s1); 95 len2=strlen(s2); 96 clr(a); 97 clr(b); 98 for(int i=0;i<len1;i++) 99 { 100 a[len1-i-1]=s1[i]-'0'; 101 } 102 for(int i=0;i<len2;i++) 103 { 104 b[len2-i-1]=s2[i]-'0'; 105 } 106 len=max(len1,len2); 107 //取长度较长者作为长度,并将长度变为2…^(k+1) 108 changelen(len); 109 //将两个整数序列复制到复数序列中 110 copyed(a,multi1,len); 111 copyed(b,multi2,len); 112 //对两个复数序列进行DFT,变为点值表示 113 fft(multi1,len,1); 114 fft(multi2,len,1); 115 //对应点点值相乘 116 for(int i=0;i<len;i++) 117 multi1[i]=multi1[i]*multi2[i]; 118 //将的出来的点值表示进行IDFT变为系数表示 119 fft(multi1,len,-1); 120 //四舍五入减小损失精度 121 for(int i=0;i<len;i++) 122 { 123 a[i]=(int)(multi1[i].r+0.5); 124 } 125 //进位 126 for(int i=0;i<len;i++) 127 { 128 a[i+1]=a[i+1]+a[i]/10; 129 a[i]%=10; 130 } 131 len=len1+len2-1; 132 //去掉前导0 133 while(a[len]<=0 && len>0) len--; 134 for(int i=len;i>=0;i--) 135 printf("%d",a[i]); 136 printf("\n"); 137 } 138 return 0; 139 }