傅里叶大数相乘算法HDU A * B Problem Plus

HDU A * B Problem Plus

这道题需要利用傅里叶变换的卷积性质：对于两个离散的序列求卷积的过程事实上就是一个求不进位乘法的过程，傅里叶变换的卷积性质告诉我们，两个离散序列的卷积运算经一个离散傅里叶变换（DFT）后，就变成了这两个序列对应每个位上面的乘法。快速傅里叶变换（FFT）就是DFT的一个O（nlogn）的快速算法，因此两个长乘法可以经FFT简化至O(nlogn)的时间级了。

 

#include <stdio.h>
#include<stdio.h>
#include<stdlib.h>
#include<string.h>
#include<math.h>
#define Maxn 200000
#define PI 3.141592653
typedef struct Virt {
double r;
    double i;
}Virt;
Virt bitRevTemp[2*Maxn];
Virt va[Maxn*2],vb[2*Maxn];
int ans[Maxn*2];
Virt Virt_Add(Virt a,Virt b)
{
    Virt t;
    t.r = a.r+b.r;
    t.i = a.i+b.i;
    return t;
}
Virt Virt_Sub(Virt a,Virt b)
{
    Virt t;
    t.r = a.r-b.r;
    t.i = a.i-b.i;
    return t;
}
Virt Virt_Multi(Virt a,Virt b)
{
    Virt t;
    t.r = a.r*b.r - a.i*b.i;
    t.i = a.r*b.i + a.i*b.r;
    return t;
}
int rev(int x,int l)   //位反置：将长度l的那几位倒过来就是
{
    int i,sum;
    sum=0;
    while (l--) {
        sum |= (x&1);
        x>>=1;
        sum<<=1;
    }
    sum>>=1;
   return sum;
}
void bitRevCpy(Virt *a,int len)
{
    int k,bits,t;
    t = 1;
    bits = 0;
    while (t<len){       //计算位长度
        ++bits;
        t = t<<1;
    }
    for (k = 0;k < len; k++)
        bitRevTemp[rev(k,bits)] = a[k];
    for (k = 0; k < len; k++)
        a[k] = bitRevTemp[k];
    return;
}
void FFT(Virt *Fourier,int len,int on)
{
    int i,j,k,loop,m;
    Virt u,t;
    Virt wm,w;
    bitRevCpy(Fourier,len);                                 //利用位反置生成输入序列
    for (loop = 2; loop <= len; loop = loop<<1) {
        m = loop;                                           //分治后计算长度为m的DFT
        wm.r = cos(on*2*PI/m);wm.i = sin(on*2*PI/m);        //这次单位复根 e^(2*pi/m) 用欧拉公式展开
        for (k = 0; k < len; k+=m) {                        //FFT过程
            w.r = 1; w.i = 0;                               //旋转因子
            for (j = 0; j < m/2; j++) {
                t = Virt_Multi(w,Fourier[k+j+m/2]);
                u = Fourier[k+j];
                Fourier[k+j] = Virt_Add(u,t);
                Fourier[k+j+m/2] = Virt_Sub(u,t);           //蝴蝶操作合并
                w = Virt_Multi(w,wm);                       //更新旋转因子
            }
        }
    }
    if (on==-1) {
        for (i = 0; i < len; i++)
            Fourier[i].r = Fourier[i].r/len;   //IDFT
    }
}
void conv(Virt *a,Virt *b,int len)
{   //求卷积
    FFT(a,len,1);
    FFT(b,len,1);
    for (int i = 0; i < len; i++)
        a[i] = Virt_Multi(a[i],b[i]);
    FFT(a,len,-1);
}
int main()
{
    char str1[Maxn*2],str2[Maxn*2];
    int i,j,top;
    int lena,lenb,lenc;
    while (scanf("%s%s",str1,str2)!=EOF)
    {
        lena = strlen(str1);lenb = strlen(str2);
        lenc = 1;
        while (lenc < lena*2 || lenc < lenb*2 ) {lenc = lenc<<1;}
        for (i = 0; str1[i]!='\0';i++) {
            va[i].r = str1[lena-i-1]-'0';
            va[i].i = 0.0;
        }
        while (i < lenc) {
            va[i].r = va[i].i = 0.0;
            i++;
        }
        for (i = 0; str2[i]!='\0';i++) {
            vb[i].r = str2[lenb-i-1]-'0';
            vb[i].i = 0.0;
        }
        while (i < lenc) {
            vb[i].r = vb[i].i = 0;
            i++;
        }
        conv(va,vb,lenc);    //FFT卷积相当于不进位乘法
        for (i = 0; i < lenc; i++)
            ans[i] = va[i].r+0.5;
        for (i = 0; i < lenc; i++) {
            ans[i+1] += ans[i]/10;
            ans[i]%=10;
        }
        for (top = lena+lenb-1;ans[top]<=0 && top>0; top--);
        top+=1;
        while (top-- >= 1)
            printf("%d",ans[top]);
        printf("\n");
    }
    return 0;
}