FFT快速傅里叶变换以O(NlogN)的时间复杂度实现大数乘

任意一个整数均能表示成An*10^(n-1) + An-1*10^(n-2) + ... + A2*10^2 + A1*10 + A0的形式，视10为自变量X，则化为一个多项式。两数相乘转化为两多项式相乘。以系数表示法表示的多项式相乘其复杂度为N^2，若采用点值表示法，结合适当的点的选取，能实现O(NlogN)的算法。

若一个多项式的最高次为N-1，那么取N个点对(xi,yi)就能够唯一确定这个多项式，其中各xi相异，yi为对应自变量xi的多项式的值。此时选取N次单位复根即可。具体定理见算法导论第30章。

此时过程很简单，先用FFT分别将两个输入的多项式由系数表示法转为点值表示法，即得到系数向量a=(a 0,a 1,...,a n-1)对应的离散傅里叶变换（DFT）y=(y 0,y1,y n-1)，此复杂度为O(NlogN)，再相乘，此复杂度仅为O(N)，最后再用FFT将结果转化回系数表示法，输出即可。

另外算法用到了二进制平摊反转置换（亦叫位反转置换）以将原本的递归过程转为迭代。

#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;
#define pi acos(-1.0)
struct complex{
    double r,i;
    complex(double real=0.0,double image=0.0)
    {
        r=real;
        i=image;
    }
    inline complex operator + (const complex a)
    {return complex(r+a.r,i+a.i);}
    inline complex operator - (const complex a)
    {return complex(r-a.r,i-a.i);}
    inline complex operator * (const complex a)
    {return complex(r*a.r-i*a.i,r*a.i+i*a.r);}
};
void brc(complex *y,int L) {
    int i,j,k;
    for (i=1,j=L>>1; i<L-1; ++i) { // 二进制平摊反转置换 O(NlogN)
        if (i < j) swap(y[i],y[j]);
        k = L>>1;
        while (j >= k) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
}
void FFT(complex *y,int L,double isI)
{
    //isI为1是DFT，-1则IDFT
    register int h,i,k;
    complex u,t;
    brc(y,L);
    for(h=2;h<=L;h<<=1)//层数自底向上
    {
        //初始化单位复根
        complex Wn( cos(isI*2*pi/h),sin(isI*2*pi/h) );
        for(j=0;j<L;j+=h)  // 原序列被分成了L/h段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(isI==-1)
        for(i=0;i<L;i++)
            y[i].r/=L;
}
const int N = 50024;
int ans[N<<2];
complex a[N<<2],b[N<<2];
char num1[N],num2[N];

int main()
{
    int i,j;
    while(~scanf("%s%s",num1,num2))
    {
        memset(ans,sizeof(ans));
        int len1=strlen(num1),len2=strlen(num2),L=1;
        int ll=len1+len2-1;
        while(L<ll) L<<=1;
        for(i=len1-1,j=0;i>=0;--i,++j)
        {
            a[j]=complex(num1[i]-'0',0); //高位实为多项式的低位
        }
        for(i=len2-1,++j)
        {
            b[j]=complex(num2[i]-'0',0);
        }
        for(i=len1;i<L;++i) a[i]=complex(0,0);//补0
        for(i=len2;i<L;++i) b[i]=complex(0,0);
        FFT(a,L,1);
        FFT(b,1);
        for(i=0;i<L;++i)
        {
            a[i]=a[i]*b[i];
        }
        FFT(a,-1);
        for(i=0;i<L;++i) ans[i]=a[i].r+0.5;
        for(i=0;i<L;++i)
        {
            ans[i+1]+=ans[i]/10; //进位
            ans[i]%=10;
        }
        int p=L;
        while(!ans[p] && p) --p;
        while(p>=0) printf("%d",ans[p--]);
        puts("");
    }
    return 0;
}