超大数乘法---FFT
思路:
算法导论第30章有详细说明。此处只是简略说明其主要的步骤。
一个知识点是：
A(x)=a0+a1x+a2x2+a3x3+……+an-1xn-1
A[0](x)=a0+a2x+a4x2+……+an-2xn/2-1
A[1](x)=a1+a3x+a5x2+……+an-1xn/2-1
A[0](x2)+x*A[1](x2)=A(x)
以上是 二进制平摊反转置换跟求和的主要式子。
多项式有两种表示形式：点值表示，系数表示。
快速FFT主要有以下四点：
1. 使次数界（上界）增加一倍。A(x)、B(x)的长度扩充到2*n
2. 求值。主要是求点值表示A(x)、B(x)的点值表示
3. 点乘。C(x)=A(x)*B(x)
4. 插值。对C(x)进行插值，求出其系数表示。
代码如下：
1 #include<iostream> 2 #include<cstdio> 3 #include<cmath> 4 #define Pi acos(-1.0)//定义Pi的值 5 #define N 200000 6 using namespace std; 7 struct complex //定义复数结构体 8 { 9 double re,im; 10 complex(double r=0.0,double i=0.0) 11 { re=r,im=i; } //初始化
//定义三种运算 12 complex operator +(complex o) 13 { return complex(re+o.re,im+o.im);} 14 complex operator -(complex o) 15 { return complex(re-o.re,im-o.im);} 16 complex operator *(complex o) 17 { return complex(re*o.re-im*o.im,re*o.im+im*o.re);} 18 }x1[N],x2[N]; 19 char a[N/2],b[N/2]; 20 int sum[N]; //存储最后的结果 21 22 void BRC(complex *y,int len)//二进制反转倒置 23 { 24 int i,j,k; 25 for(i=1,j=len/2;i<len-1;i++) 26 { 27 if(i<j)swap(y[i],y[j]);//i<j保证只交换一次 28 k=len/2; 29 while(j>=k) 30 { 31 j-=k;k=k/2; 32 } 33 if(j<k)j+=k; 34 } 35 } 36 void FFT(complex *y,int len ,double on)//on=1表示顺，-1表示逆 37 { 38 int i,j,k,h; 39 complex u,t; 40 BRC(y,len); 41 for(h=2;h<=len;h<<=1)//控制层数 42 { 43 //初始化单位复根 44 complex wn(cos(on*2*Pi/h),sin(on*2*Pi/h)); 45 for(j=0;j<len;j+=h)//控制起始下标 46 { 47 //初始化螺旋因子 48 complex w(1,0); 49 for(k=j;k<j+h/2;k++) 50 { 51 u=y[k]; 52 t=w*y[k+h/2]; 53 y[k]=u+t; 54 y[k+h/2]=u-t; 55 w=w*wn;//更新螺旋因子 56 } 57 } 58 } 59 if(on==-1) 60 for(i=0;i<len;i++) //逆FFT(IDFT) 61 y[i].re/=len; 62 63 } 64 int main() 65 { 66 int len1,len2,len,i; 67 while(scanf("%s%s",a,b)!=EOF) 68 { 69 len1=strlen(a); 70 len2=strlen(b); 71 len=1; 72 //扩充次数界至2*n 73 while(len<2*len1||len<2*len2)len<<=1;//右移相当于len=len*2 74 //倒置存储 75 for(i=0;i<len1;i++) 76 { x1[i].re=a[len1-1-i]-'0';x1[i].im=0.0;} 77 for(;i<len1;i++) //多余次数界初始化为0 78 {x1[i].re=x1[i].im=0.0;} 79 for(i=0;i<len2;i++) 80 { x2[i].re=b[len2-1-i]-'0';x2[i].im=0.0;} 81 for(;i<len2;i++) //多余次数界初始化为0 82 {x2[i].re=x2[i].im=0.0;} 83 //FFT求值 84 FFT(x1,len,1);//FFT(a) 1表示顺 -1表示逆 85 FFT(x2,len,1);//FFT(b) 86 //点乘，结果存入x1 87 for(i=0;i<len;i++) 88 x1[i]=x1[i]*x2[i]; 89 //插值，逆FFT（IDTF） 90 FFT(x1,len,-1); 91 92 //细节处理 93 for(i=0;i<len;i++) 94 sum[i]=x1[i].re+0.5;//四舍五入 95 for(i=0;i<len;i++) //进位 96 { 97 sum[i+1]+=sum[i]/10; 98 sum[i]%=10; 99 } 100 //输出 101 len=len1+len2-1; 102 while(sum[len]<=0&&len>0)len--;//检索最高位 103 for(i=len;i>=0;i--) //倒序输出 104 cout<<sum[i]; 105 cout<<endl; 106 } 107 return 0; 108 }
值得注意的细节：a,b输入的是字符串，x1,x2的re中存储的是double类型的数据。
