离散傅里叶变换计算乘法

利用傅立叶变换可以把大数乘法时间控制在：O（n*lgn)，首先把整数a1a2...an看作多项式f(x) = a1x + a2x^2 + .... an x^n(其中x取10），把整数乘法转换为多项式乘法，多项式一般乘法也要O(n*n)时间，而多项式点值表示法的乘法只需要O(n)的时间，于是就有了条折返路线，如下：

1.首先把f(x)扩展为2n次多项式，取2n个点求多项式值  (O(nlgn))

2.在对这两的多项式的点值表示进行乘积     （O（n))

3.在对所得的点值进行插值得到结果多项式的系数   (O(nlgn))

过程就是：系数  -> 点值  -> 系数

第一步如果随意取2n个点进行求值，也要O（n*n)，而如果利用单位复根exp(iu)的诸多性质就可以在O(nlgn)时间内求的点值对，也可以在O(nlgn)时间内插值，具体可以参考《计算方法引论》中的证明。

代码：

 

#include <cstdio>
#include <cstring>
#include <cmath>

using namespace std;

#define MAXN 500005
#define PI acos(-1.0)

struct vir    
{
    double r,i;
    
    vir(double r = 0.0,double i = 0.0)
    {
        this -> r = r;
        this -> i = i;
    }
    vir operator +(const vir & b)
    {
        return vir(r + b.r,i + b.i);
    }
    vir operator -(const vir & b)
    {
        return vir(r - b.r,i - b.i);
    }
    vir operator *(const vir & b)
    {
        return vir(r * b.r - i * b.i,r * b.i + i * b.r);
    }
};

char stra[MAXN],strb[MAXN];
vir vira[MAXN],virb[MAXN];
int ans[MAXN];

int init(char *stra,char *strb,vir *vira,vir *virb)
{
    int lena = strlen(stra);
    int lenb = strlen(strb);
    int len = 1;
    while(len < 2 * lena || len < 2 * lenb) len <<= 1;    //位数扩展2倍
    
    int i;
    for(i = 0;i < lena;i++)
    {
        vira[i].r = stra[lena-i-1] - '0';
        vira[i].i = 0.0;
    }
    for(;i < len;i++)
    {    
        vira[i].r = vira[i].i = 0.0;
    }
    for(i = 0;i < lenb;i++)
    {
        virb[i].r = strb[lenb-i-1] - '0';
        virb[i].i = 0.0;
    }
    for(;i < len;i++)
    {
        virb[i].r = virb[i].i = 0.0;
    }
    return len;
}

void bitReverseSort(vir * v,int len)        //按逆序数排序
{
    for(int i = 0;i < len;i++)
    {
        int n = i;
        int revn = 0;
        for(int j = 1;j < len;j <<= 1)
        {
            revn = (revn << 1) | (n & 1);
            n >>= 1;
        }
        if(revn > i)
        {
            vir t = v[i];
            v[i] = v[revn];
            v[revn] = t;
        }
    }
}

void FFT(vir *v,int len,int flag)        //flag为1求值，-1插值
{
    bitReverseSort(v,len);    

    for(int r = 1;r < len;r <<= 1)
    {
        vir wp(cos(PI / r),sin(-flag * PI / r));
        vir w(1,0);
        for(int s = 0;s < r;s++)
        {
            for(int k = s;k < len;k += 2 * r)
            {
                int l = k + r;
                vir d = w * v[l];
                v[l] = v[k] - d;
                v[k] = v[k] + d;
            }
            w = w * wp;
        }
    }
    
    if(flag == -1)
    {
        for(int i = 0;i < len;i++)
        {
            v[i].r /= len;
        }
    }
}

void mul(vir * vira,vir * virb,int len)
{
    for(int i = 0;i < len;i++)
    {
        vira[i] = vira[i] * virb[i];
    }
}

void output(vir *vira,int *ans,int len)
{
    for(int i = 0;i < len;i++)
    {
        ans[i] = vira[i].r + 0.5;
    }
    for(int i=0;i<len;i++)
        {
          ans[i+1] += ans[i]/10;
        ans[i] %= 10;
        }
    int pos = len - 1;
    while(!ans[pos] && pos > 0)pos--;
    for(int i = pos;i >= 0;i--)
    {
        printf("%d",ans[i]);
    }
    printf("\n");
}

int main()
{
    while(scanf("%s%s",stra,strb) != EOF)
    {
        int len = init(stra,strb,vira,virb);
        FFT(vira,len,1);
        FFT(virb,len,1);
        mul(vira,virb,len);
        FFT(vira,len,-1);
        output(vira,ans,len);
    }
    return 0;
}

//9416236    2013-10-26 10:13:27    Accepted    1402    265MS    16512K    2836 B    G++    超级旅行者

 

参考资料：

 

《计算方法引论》 第三版 徐萃薇

《算法导论》