A * B Problem Plus

A * B Problem Plus

题目链接：http://acm.split.hdu.edu.cn/showproblem.php?pid=1402

FFT

（FFT的详细证明参见算法导论第三十章）

一个多项式有两种表达方式：

1.系数表示法，系数表示的多项式相乘，时间复杂度为O（n^2）；

2.点值表示法，点值表示的多项式相乘，时间复杂度为O（n）.

简单的说，FFT能办到的就是将系数表示的多项式转化为点值表示，其时间复杂度为O（nlgn），而将点值表示的多项式转化为系数表示需要IFFT（FFT的逆运算），其形式与FFT相似，时间复杂度也为O（nlgn）.

这道题需要用FFT将两个大数转化为点值表示，相乘后再用IFFT将点值表示转化回系数表示，总时间复杂度为O（nlgn）.

代码如下：

#include<cstdio>
#include<cmath>
#include<algorithm>
#include<cstring>
#include<iostream>
#define N 200005
using namespace std;
const double pi=acos(-1.0);
struct Complex{
    double r,i;
    Complex(double r=0,double i=0):r(r),i(i){};
    Complex operator + (const Complex &rhs){
        return Complex(r+rhs.r,i+rhs.i);
    }
    Complex operator - (const Complex &rhs){
        return Complex(r-rhs.r,i-rhs.i);
    }
    Complex operator * (const Complex &rhs){
        return Complex(r*rhs.r-i*rhs.i,i*rhs.r+r*rhs.i);
    }
}a[N],b[N],c[N];
char s1[N],s2[N];
int ans[N],n1,n2,len;
inline void sincos(double theta,double &p0,double &p1){
    p0=sin(theta);
    p1=cos(theta);
}
void FFT(Complex P[], int n, int oper){
    for(int i=1,j=0;i<n-1;i++){
        for(int s=n;j^=s>>=1,~j&s;);
        if(i<j)swap(P[i],P[j]);
    }
    Complex unit_p0;
    for(int d=0;(1<<d)<n;d++){
        int m=1<<d,m2=m*2;
        double p0=pi/m*oper;
        sincos(p0,unit_p0.i,unit_p0.r);
        for(int i=0;i<n;i+=m2){
            Complex unit=1;
            for(int j=0;j<m;j++){
                Complex &P1=P[i+j+m],&P2=P[i+j];
                Complex t=unit*P1;
                P1=P2-t;
                P2=P2+t;
                unit=unit*unit_p0;
            }
        }
    }
    if(oper==-1)for(int i=0;i<len;i++)P[i].r/=len;
}
void Conv(Complex a[],Complex b[],int len){//求卷积
    FFT(a,len,1);//FFT
    FFT(b,len,1);//FFT
    for(int i=0;i<len;++i)c[i]=a[i]*b[i];
    FFT(c,len,-1);//IFFT
}
void init(char *s1,char *s2){
    len=1;
    n1=strlen(s1),n2=strlen(s2);
    while(len<2*n1||len<2*n2)len<<=1;
    int idx;
    for(idx=0;idx<n1;++idx){
        a[idx].r=s1[n1-1-idx]-'0';
        a[idx].i=0;
    }
    while(idx<len){
        a[idx].r=a[idx].i=0;
        idx++;
    }
    for(idx=0;idx<n2;++idx){
        b[idx].r=s2[n2-1-idx]-'0';
        b[idx].i=0;
    }
    while(idx<len){
        b[idx].r=b[idx].i=0;
        idx++;
    }
}
int main(void){
    while(scanf("%s%s",s1,s2)==2){
        init(s1,s2);
        Conv(a,b,len);
        for(int i=0;i<len+len;++i)ans[i]=0;//93ms
        //memset(ans,0,sizeof(ans));//140ms
        int index;
        for(index=0;index<len||ans[index];++index){
            ans[index]+=(c[index].r+0.5);
            ans[index+1]+=(ans[index]/10);
            ans[index]%=10;
        }
        while(index>0&&!ans[index])index--;
        for(;index>=0;--index)printf("%d",ans[index]);
        printf("\n");
    }
}