HDU 1402 A * B Problem Plus FFT

 

大数乘法,暴力逐项乘后相加 O(n^2) 或者 用JAVA的BigInteger都超时

这时候就要用到FFT了...多项式乘法,大数乘法都能在O(nlogn)时间内解决

因为之前没有复变和数字信号课程的基础...所以今天看FFT异常的吃力...

一个上午都在研究复数的概念看FFT的定义,下午的时候才开始看FFT的算法...晚上就学了一个别人的板...过了这个水题

发现FFT完全就是模板,直接会用就行了...

但是这里还是讲解一下原理: 

首先这个课件应该是最能解决概念性问题的: http://wenku.baidu.com/view/14457119de80d4d8d15a4f34.html

然后就是kuangbin神推荐的这个课件: http://wenku.baidu.com/view/8bfb0bd476a20029bd642d85.html

这题对于FFT,我们一共就要做三件事 -> 对a,b求值(做DFT) , 两个值点乘 , 对结果插值(IDFT) 

其中值得一提的也就是求值插值的过程了...(插值是求值的逆过程,也就是对Vandermonde矩阵求逆后的DFT,具体的算导上写的很详细...结论在在515页上面的30.11式,和30.08式很像对吧 :>

递归法的求值在算导513页,迭代法是在517页,其中递归法我暂时还没去敲这个板...我先讲解一下楼下要贴的迭代的板

首先一个难点是反转置换,这个太神奇了...我无法解释为什么反转下标的二进制能够达到这个效果...举个图中的例子 4是100 在图中的位置是001 , 6是110 在图中的位置是011

通过这个反转操作, 我们就可以自底向上进行算式合并了....下面这个图很直观的说明了反转的重要性

这个板很逆天的给了一个我看不懂的过程...具体实现如下

void rev(complex *y,int len) // logn
{
    // 对向量反转置换 见算导517页
    register 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;
    }
}

然后就是求值的过程了...

也就两个算式进行一下蝴蝶操作 ... (名字很屌...其实很好理解就是交叉加减

前n/2个项的算式:

后n/2个项的算式: 

你可以发现前后两个算式除了符号不同以外其他都一样...

我们将其中共同的算式提取出来,这也就是俗称的蝴蝶操作...

关于单位复根和Vandermonde矩阵我就不做多余的解释了...可以参见算导或者上面的课件

以下是代码

void fft(complex *y , int len , double op) // nlogn
{
    // a为向量, h为向量长度, OP=1时求值 OP=-1时插值
    register int i,j,k,h;
    complex u,t;
    rev(y,len); // 对向量反转置换
    for(h = 2;h <= len; h <<= 1) // 控制层数
    {
        complex wn(cos(op*2*PI/h),sin(op*2*PI/h));//该层单位复根,h为层数
        for(j = 0 ; j < len ; j += 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(op == -1)
        for(i = 0 ; i < len ; i++)
            y[i].real /= len; // IDFT
}

不是自己的代码看起来还是不爽,今天再推掉重写一遍好了...

HDU1402:

struct complex
{
    double real,imag;
    bool vis;
    complex(double a = 0.,double b = 0.){
        real = a;
        imag = b;
    }
    complex operator + (complex e){
        return complex(real+e.real,imag+e.imag);
    }
    complex operator - (complex e){
        return complex(real-e.real,imag-e.imag);
    }
    complex operator * (complex e){
        return complex(real*e.real-imag*e.imag,real*e.imag+imag*e.real);
    }
}x1[MAXN],x2[MAXN];
char a[MAXN],b[MAXN];
int sum[MAXN];
void rev(complex *y,int len) // logn
{
    // 对向量反转置换 见算导517页
    register 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(complex *y , int len , double op) // nlogn
{
    // a为向量, h为向量长度, OP=1时求值 OP=-1时插值
    register int i,j,k,h;
    complex u,t;
    rev(y,len); // 对向量反转置换
    for(h = 2;h <= len; h <<= 1) // 控制层数
    {
        complex wn(cos(op*2*PI/h),sin(op*2*PI/h));//该层单位复根,h为层数
        for(j = 0 ; j < len ; j += 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(op == -1)
        for(i = 0 ; i < len ; i++)
            y[i].real /= len; // IDFT
}

int main()
{
    //freopen("in.txt","r",stdin);
    //CHEAT;
    while(~scanf("%s%s",a,b))
    {
        memset(sum,0,sizeof(sum));
        int i ;
        int len1 = strlen(a);
        int len2 = strlen(b);
        int len = 1 ;
        while(len < len1 * 2 || len < len2 * 2) len <<= 1; //求最大次界
        //cout<<len<<endl;
        //倒置 将多余次数界初始化为0
        for(i = 0 ; i < len1 ; i++)
        {
            x1[i].real = a[len1-i-1] - '0';
            x1[i].imag = 0.0;
        }
        for(; i < len ; i++)
        {
            x1[i].real = 0.0;
            x1[i].imag = 0.0;
        }
        for(i = 0 ; i < len2 ; i++)
        {
            x2[i].real = b[len2-i-1] - '0';
            x2[i].imag = 0.0;
        }
        for(; i < len ; i++)
        {
            x2[i].real = 0.0;
            x2[i].imag = 0.0;
        }
        //对x1,x2求值
        fft(x1,len,1);
        fft(x2,len,1);
        for(i = 0 ; i < len ; i++)    x1[i] = x1[i] * x2[i]; // 点乘
        //对乘积插值
        fft(x1,len,-1);
        for(i = 0 ; i < len ; i++)
            sum[i] = x1[i].real + 0.5; // 四舍五入
        for(i = 0 ; i < len ; i++) // 进位
        {
            sum[i+1] += sum[i] / 10;
            sum[i] %= 10;
        }
        len = len1 + len2 - 1;
        while(sum[len] <= 0 && len > 0)    len--; // 除前导0
        for(i = len ; i >= 0 ; i--) printf("%c",sum[i]+'0');
        puts("");
    }
    return 0;
}