hdu1402 A*B 快速傅里叶变换(FFT)

      两个不超过50000位的数，求乘积。高精度乘法，但是写暴力高精的话复杂度是O(n^2)，不压位估计是要T掉的=,这里介绍一种新的方法，可以在O(nlogn)的复杂度内求出答案。

先来说一下我对fft的理解吧，fft其实就是一个求多项式乘法的快速算法，两个n阶多项式相乘，传统的方法是循环相乘再累加，复杂度是O（n^2），而用fft去实现的话，复杂度可以降到O(nlogn)。对于n阶多项式，我们可以用两种方法去表示，一种是常用的稀疏表示法，也就是a0+a1*x^1+a2*x^2+...+an*x^n，另一种是点集表示法，在平面上取x坐标的集合X，对于每个x[i]求出y[i]=A（x)，其中A（）为一个多项式，这样可以用这个点集来描述这个多项式（可以看做是集合表示吧）。对于用系数表示两个多项式X,Y，求Z=X*Y的话，可以通过一下四步实现：

1：使次数界增加一倍：把X，Y的系数项扩大一倍，扩充成2n阶的多项式。

2: 求值：把X和Y分别转化成用的点集表示法表示。

3：点乘：Z[i]=X[i]*Y[i],乘法可以看做复数乘法。

4：插值：把Z由点集表示法转化成稀疏表示法。

其中过程2即离散傅里叶变化（DFT），过程4为逆离散傅里叶变换（IDFT），而如果有选择性的选取某些点来构成点集表示发的集合（即选取复平面上的“单位复根”），就可以将DFT优化成O(nlogn)的FFT从而降低整体的时间复杂度，最终过程2,4的时间复杂度为O(nlogn),过程1,3的复杂度为O（n）,整体的时间复杂度为O（nlogn）.FFT具体的具体实现可以参考算导上的讲解，当然看不明白的话就只好记模板了= -...

然后来看这一题，两个数相乘的话，可以看做x[i]=10^i的两个多项式的乘法，那么可以直接用FFT来实现，要注意的一点就是求出乘积的系数后，转化成整数时注意四舍五入，然后注意处理进位的情况。

#include <iostream>
#include <cstdio>
#include <algorithm>
#include <cmath>
#include <cstring>
using namespace std;
typedef long long ll;
const double PI = acos(-1.0);
struct comp
{
    double r,i;
    comp(double rt=0,double it=0)
    {
        r=rt;
        i=it;
    }
    comp operator +(const comp& b)
    {
        return comp(r+b.r,i+b.i);
    }
    comp operator -(const comp &b)
    {
        return comp(r-b.r,i-b.i);
    }
    comp operator *(const comp &b)
    {
        return comp(r*b.r-i*b.i,r*b.i+i*b.r);
    }
};

void change(comp 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]);
        k = len/2;
        while( j >= k)
        {
            j -= k;
            k /= 2;
        }
        if(j < k)j += k;
    }
}

void fft(comp y[],int len,int on)
/* on=1 DFT 把一个多项式的系数向量转化为点集表示；
on=-1,IDFT 把一个点集转化成多项式的系数向量*/
{
    change(y,len);
    for(int h = 2;h <= len;h <<= 1)
    {
        comp wn(cos(-on*2*PI/h),sin(-on*2*PI/h));
        for(int j = 0;j < len;j += h)
        {
            comp w(1,0);
            for(int k = j;k < j+h/2;k++)
            {
                comp u = y[k];
                comp t = w*y[k+h/2];
                y[k] = u+t;
                y[k+h/2] = u-t;
                w = w*wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0;i < len;i++)
            y[i].r /= len;
}

void conv(comp f[],int len)//求f的卷积
{
    fft(f,len,1);
    for (int i=0; i<len; i++)
    f[i]=f[i]*f[i];
    fft(f,len,-1);
}
int n,m,len1,len2,len;
char a[50500],b[50500];
comp x[450500],y[450500];
int ans[450500];
int main()
{
//    freopen("in.txt","r",stdin);
    while(~scanf("%s",a))
    {
        scanf("%s",b);
        int l1=strlen(a);
        for (int i=0; i<l1; i++)
        x[l1-i-1]=comp(a[i]-'0',0);

        int l2=strlen(b);
        for (int i=0; i<l2; i++)
        y[l2-i-1]=comp(b[i]-'0',0);

        len=1;
        while(len<(l1+l2)*2) len<<=1;
        for (int i=l1; i<len; i++) x[i]=comp(0,0);
        for (int i=l2; i<len; i++) y[i]=comp(0,0);

        fft(x,len,1);
        fft(y,len,1);
        for (int i=0; i<len; i++)
        x[i]=x[i]*y[i];
        fft(x,len,-1);

        for (int i=0; i<len; i++)
        ans[i]=(int)(x[i].r+0.5);
        for (int i=0; i<len; i++)
        if (ans[i]>9)
        {
            ans[i+1]+=ans[i]/10;
            ans[i]%=10;
        }

        bool ok=false;
        for (int i=len-1; i>=0; i--)
        {
            if (ok) printf("%d",ans[i]);
            else if (ans[i])
            {
                ok=true;
                printf("%d",ans[i]);
            }
        }
        if (!ok) puts("0");
        else puts("");


    }
}