HDU - 1402 A * B Problem Plus (FFT实现高精度乘法)
题意：计算A*B，A,B均为长度小于50000的整数。
这是FFT在大整数相乘中的一个应用，我本来想用NTT做的，但NTT由于取模很可能取炸，所以base必须设得很小，而且效率也比不上FFT。
A和B的存储均用long long，在计算乘积的时候转化成double，计算完成后再转回来即可。
测得base在精度允许范围内最多能开到10000。
把平方和快速幂的函数也写上了，可以当模板用~
1 #include<bits/stdc++.h> 2 using namespace std; 3 typedef long long ll; 4 typedef double db; 5 const int N=1e5+10,inf=0x3f3f3f3f; 6 const db pi=acos(-1); 7 const ll P[]= {1,10,100,1000,10000}; 8 struct cpl { 9 db x,y; 10 cpl operator-(cpl& b) {return {x-b.x,y-b.y};} 11 cpl operator+(cpl& b) {return {x+b.x,y+b.y};} 12 cpl operator*(cpl& b) {return {x*b.x-y*b.y,x*b.y+y*b.x};} 13 } A[N],B[N]; 14 void change(cpl* a,int n) { 15 for(int i=1,j=n>>1,k; i<n-1; ++i) { 16 if(i<j)swap(a[i],a[j]); 17 k=n>>1; 18 while(j>=k)j-=k,k>>=1; 19 j+=k; 20 } 21 } 22 void FFT(cpl* a,int n,int f) { 23 change(a,n); 24 for(int k=1; k<n; k<<=1) { 25 cpl wn= {cos(pi*f/k),sin(pi*f/k)}; 26 for(int i=0; i<n; i+=k<<1) { 27 cpl w{1,0}; 28 for(int j=i; j<i+k; ++j) { 29 cpl x=a[j],y=w*a[j+k]; 30 a[j]=x+y,a[j+k]=x-y; 31 w=w*wn; 32 } 33 } 34 } 35 if(!~f)for(int i=0; i<n; ++i)a[i].x/=n,a[i].y/=n; 36 } 37 struct Bigint { 38 static const int N=1e5+10; 39 static const int bit=4; 40 static const ll base=pow(10,bit); 41 int n; 42 ll a[N]; 43 ll& operator[](int x) {return a[x];} 44 void init(ll x) { 45 memset(a,0,sizeof a); 46 a[0]=x,n=0; 47 } 48 void init(char* s) { 49 memset(a,0,sizeof a); 50 int m=strlen(s); 51 for(int i=0; i<m; ++i)a[(m-1-i)/bit]+=(s[i]-'0')*P[(m-1-i)%bit]; 52 n=(m-1)/bit; 53 } 54 void Sqr() { 55 int m; 56 for(m=1; m<=n*2; m<<=1); 57 for(int i=0; i<m; ++i)A[i]= {a[i],0}; 58 FFT(A,m,1); 59 for(int i=0; i<m; ++i)A[i]=A[i]*A[i]; 60 FFT(A,m,-1); 61 for(int i=0; i<m; ++i)a[i]=A[i].x+0.5; 62 for(int i=0; i<m; ++i) { 63 if(a[i]>=base) { 64 a[i+1]+=a[i]/base; 65 a[i]%=base; 66 if(i==m-1)++m; 67 } 68 } 69 for(n=m-1; n>1&&!a[n]; --n); 70 } 71 void Mul(Bigint& b) { 72 int m; 73 for(m=1; m<=n*2||m<=b.n*2; m<<=1); 74 for(int i=0; i<m; ++i)A[i]= {a[i],0},B[i]= {b[i],0}; 75 FFT(A,m,1),FFT(B,m,1); 76 for(int i=0; i<m; ++i)A[i]=A[i]*B[i]; 77 FFT(A,m,-1); 78 for(int i=0; i<m; ++i)a[i]=A[i].x+0.5; 79 for(int i=0; i<m; ++i) { 80 if(a[i]>=base) { 81 a[i+1]+=a[i]/base; 82 a[i]%=base; 83 if(i==m-1)++m; 84 } 85 } 86 for(n=m-1; n>1&&!a[n]; --n); 87 } 88 void Pow(int p) { 89 Bigint b=*this; 90 this->init(1); 91 for(; p; p>>=1,b.Sqr())if(p&1)this->Mul(b); 92 } 93 void pr() { 94 for(int i=n; i>=0; --i) { 95 if(i==n)printf("%lld",a[i]); 96 else printf("%04lld",a[i]); 97 } 98 puts(""); 99 } 100 } a,b; 101 char s1[N],s2[N]; 102 103 int main() { 104 while(scanf("%s%s",s1,s2)==2) { 105 a.init(s1),b.init(s2); 106 a.Mul(b); 107 a.pr(); 108 } 109 return 0; 110 }
