基于FFT算法的超大数字乘法优化技术

  1 #include
  2 #include
  3 #include
  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;i1;i++)
 26   {
 27      if(i<j)swap(y[i],y[j]);//i 28      k=len/2;
 29      while(j>=k)
 30      {
 31        j-=k;k=k/2;
 32      }
 33    if(jk;
 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//控制起始下标
 46     {
 47       //初始化螺旋因子
 48         complex w(1,0);
 49         for(k=j;k2;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//逆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)
 76      { x1[i].re=a[len1-1-i]-'0';x1[i].im=0.0;}
 77      for(;i//多余次数界初始化为0
 78      {x1[i].re=x1[i].im=0.0;}
 79      for(i=0;i)
 80      { x2[i].re=b[len2-1-i]-'0';x2[i].im=0.0;}
 81      for(;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)
 88          x1[i]=x1[i]*x2[i];
 89 //插值，逆FFT（IDTF）
 90      FFT(x1,len,-1);
 91 
 92 //细节处理
 93      for(i=0;i)
 94          sum[i]=x1[i].re+0.5;//四舍五入
 95      for(i=0;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 }