超大數乘法---FFT

思路:

   算法導論第30章有詳細說明。此處只是簡略說明其主要的步驟。

一個知識點是：

  A(x)=a₀+a₁x+a₂x²+a₃x³+……+a_n-1x^n-1

 A^[0](x)=a₀+a₂x+a₄x²+……+a_n-2x^n/2-1

 A^[1](x)=a₁+a₃x+a₅x²+……+a_n-1x^n/2-1

 

 A^[0](x²)+x*A^[1](x²)=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)進行插值，求出其系數表示。

---

 代碼如下：

 

#include<iostream>
#include<cstdio>
#include<cmath>
#define Pi acos(-1.0)//定義Pi的值
#define N 200000
using namespace std;
struct complex //定義復數結構體
{
  double re,im;
  complex(double r=0.0,double i=0.0)
  {  re=r,im=i; }  //初始化
   //定義三種運算 
  complex operator +(complex o)
  { return complex(re+o.re,im+o.im);}
  complex operator -(complex o)
  { return complex(re-o.re,im-o.im);}
  complex operator *(complex o)
  { return complex(re*o.re-im*o.im,re*o.im+im*o.re);}
}x1[N],x2[N];
char a[N/2],b[N/2];
int sum[N];    //存儲最后的結果

void BRC(complex *y,int len)//二進制反轉倒置
{
  int i,j,k;
  for(i=1,j=len/2;i<len-1;i++)
  {
     if(i<j)swap(y[i],y[j]);//i<j保證只交換一次
     k=len/2;
     while(j>=k)
     {
       j-=k;k=k/2;
     }
   if(j<k)j+=k;
  }
}
void FFT(complex *y,int len ,double on)//on=1表示順，-1表示逆
{
  int i,j,k,h;
  complex u,t;
  BRC(y,len);
  for(h=2;h<=len;h<<=1)//控制層數
  {
      //初始化單位復根
    complex wn(cos(on*2*Pi/h),sin(on*2*Pi/h));
    for(j=0;j<len;j+=h)//控制起始下標
    {
      //初始化螺旋因子
        complex w(1,0);
        for(k=j;k<j+h/2;k++)
        {
            u=y[k];
            t=w*y[k+h/2];
            y[k]=u+t;
            y[k+h/2]=u-t;
            w=w*wn;//更新螺旋因子
        }    
    }
  }
if(on==-1)
   for(i=0;i<len;i++) //逆FFT(IDFT)
       y[i].re/=len;

}
int main()
{
    int len1,len2,len,i;
  while(scanf("%s%s",a,b)!=EOF)
  { 
     len1=strlen(a);
     len2=strlen(b);
     len=1;
//擴充次數界至2*n
     while(len<2*len1||len<2*len2)len<<=1;//右移相當於len=len*2
//倒置存儲
     for(i=0;i<len1;i++)
     { x1[i].re=a[len1-1-i]-'0';x1[i].im=0.0;}
     for(;i<len1;i++)  //多余次數界初始化為0
     {x1[i].re=x1[i].im=0.0;}
     for(i=0;i<len2;i++)
     { x2[i].re=b[len2-1-i]-'0';x2[i].im=0.0;}
     for(;i<len2;i++)  //多余次數界初始化為0
     {x2[i].re=x2[i].im=0.0;}
//FFT求值
     FFT(x1,len,1);//FFT(a) 1表示順 -1表示逆
     FFT(x2,len,1);//FFT(b)
//點乘，結果存入x1
     for(i=0;i<len;i++)
         x1[i]=x1[i]*x2[i];
//插值，逆FFT（IDTF）
     FFT(x1,len,-1);

//細節處理
     for(i=0;i<len;i++)
         sum[i]=x1[i].re+0.5;//四舍五入
     for(i=0;i<len;i++)     //進位
     {
       sum[i+1]+=sum[i]/10;
       sum[i]%=10;
     }
//輸出
     len=len1+len2-1;
     while(sum[len]<=0&&len>0)len--;//檢索最高位
     for(i=len;i>=0;i--)   //倒序輸出
      cout<<sum[i];
      cout<<endl;
  }
return 0;
}

 

值得注意的細節：a,b輸入的是字符串，x1,x2的re中存儲的是double類型的數據。