• hdu


    题意:求A * B(A,B的位数不超过50000位)。

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

    ——>>在我刚上大学的时候,若遇到它,猜想我会输入2个int,然后直接用*输出;在学了大数高精度之后,若遇到它,猜想我会用大数去乘;结果还是把神题看成了水题。。。

    10场多校,FFT常常在解题报告中出现,实在不能不去看看这个传奇的FFT。开始看LJ《训练指南》里的FFT,知道了个大概,接着看张家琳的论文《多项式的乘法》,通了一根筋,在学长的指导下,去啃《算法导论》里的FFT,不愧是经典,超级详细~

    ——>>步骤(LJ):补0,求值,乘法,插值。(详细步骤就不写了)下面是根据《算法导论》的思路写出来的(迭代实现,变量名也如算导,没做太大改变)~

    #include <cstdio>
    #include <cmath>
    #include <complex>
    #include <cstring>
    #include <algorithm>
    
    using namespace std;
    
    typedef complex<double> Complex;
    const int maxn = 1 << 17;       //50000 * 2 = 10^5 < 2^17
    const double pi = acos(-1);
    
    char sa[maxn], sb[maxn];        //输入的2个乘数
    Complex a[maxn], b[maxn], c[maxn], A[maxn], B[maxn], C[maxn];       //a,b,c为系数序列,A,B为点值对的值序列,C为A,B的乘积
    int ans[maxn], n, lay;      //ans为结果,n为扩展后的次数界,lay为层数-1
    
    int rev(int x){     //位反转置换
        int ret = 0, i;
        for(i = 0; i < lay; i++) if(x & (1 << i)) ret += (1 << (lay-1-i));
        return ret;
    }
    
    void bit_reverse_copy(Complex *a, Complex *A){      //把原序列排成迭代所需的顺序
        for(int i = 0; i < n; i++) A[rev(i)] = a[i];
    }
    
    void FFT(Complex *a, Complex *A, bool IDFT = 0){        //FFT,系数序列a,值序列A,是DFT还是IDFT
        bit_reverse_copy(a, A);     //最底层的值序列
        for(int s = 1; s <= lay; s++){      //lay次迭代
            int m = 1 << s;     //间隔
            double uu = (IDFT ? -1 : 1) * 2.0 * pi / m;     //DFT用正,IDFT用负
            Complex wm(cos(uu), sin(uu));       //主m次单位根
            for(int k = 0; k < n; k += m){      //每一层的偶序列(左边序列)的第1个元素
                Complex w(1);       //旋转因子
                for(int j = 0; j < m/2; j++){       //遍历子序列(偶序列,奇序列)的每1个元素
                    Complex u = A[k + j];       //取子(偶)序列的第j个元素
                    Complex t = w * A[k + j + m/2];     //计算上一行中对应的奇序列的值
                    A[k + j] = u + t;       //迭代计算到上一层
                    A[k + j + m/2] = u - t;     //根据分治原理和相消引理推出
                    w *= wm;        //旋转因子更新
                }
            }
        }
        if(IDFT) for(int i = 0; i < n; i++) A[i] /= n;      //如果是IDFT,记得除以n
    }
    
    void init(){
        int alen = strlen(sa), blen = strlen(sb), i;
        for(i = alen-1; i >= 0; i--) a[alen-1-i] = Complex(sa[i] - '0', 0);     //得到第1个乘数的系数序列
        for(i = blen-1; i >= 0; i--) b[blen-1-i] = Complex(sb[i] - '0', 0);     //得到第2个乘数的系数序列
        n = 1;
        lay = 0;        //迭代的次数
        while(n < alen + blen){     //n为乘积的次数界,且n需要是2的次幂
            n <<= 1;
            lay++;
        }
        for(i = alen; i < n; i++) a[i] = Complex(0);     //扩展到次数界为n,即到n-1次
        for(i = blen; i < n; i++) b[i] = Complex(0);     //扩展到次数界为n,即到n-1次
    }
    
    void solve(){
        int i;
        FFT(a, A);
        FFT(b, B);
        for(i = 0; i < n; i++) C[i] = A[i] * B[i];      //求积
        FFT(C, c, 1);
        memset(ans, 0, sizeof(ans));
        for(i = 0; i < n; i++){     //系数不能 > 9, >= 10的进位
            int x = (int)round(c[i].real());
            ans[i] += x;
            ans[i+1] += ans[i] / 10;        //用=会WA
            ans[i] %= 10;
        }
        for(i = n-1; !ans[i] && i > 0; i--);
        for(; i >= 0; i--) printf("%d", ans[i]);
        puts("");
    }
    
    int main()
    {
        while(scanf("%s%s", sa, sb) == 2){
            init();
            solve();
        }
        return 0;
    }
    



  • 相关阅读:
    Asp.net MVC 3实例学习之ExtShop(三)——完成首页
    Asp.net MVC 3实例学习之ExtShop(一)————创建应用并设置开发环境
    JDK里的设计模式
    Linux系统下ssh的相关配置
    简单介绍asp模式与saas模式
    Linux系统中xorg.conf文件简介
    [C++] MurmurHash2的性能
    [C++] 在程序里调用DOS命令
    Android SDK 1.5中文版 (Application基础—1)
    linux系统单网卡绑定双IP的方法
  • 原文地址:https://www.cnblogs.com/riskyer/p/3279941.html
Copyright © 2020-2023  润新知