hdu 1402 A * B Problem Plus (FFT&DFT)

http://acm.hdu.edu.cn/showproblem.php?pid=1402

　　大数乘法，将乘法转换成多项式来解，用DFT加上分治法，将O(n^2)的复杂度降至O(n logn).

　　在FFT（快速傅里叶变换）中，构造旋转因子并且利用蝴蝶操作将原来的每个系数离散化，然后将两个多项式对应的系数乘起来。因为傅里叶变换是一个可逆的操作，所以最后IDFT，将答案的每个系数还原，最后输出结果。

　　目前我对FFT的机理还没完全理解，只能明白是利用DFT（离散傅里叶变换）将一个多项式离散化，构造出傅里叶级数，然后再利用其周期性的性质，对应系数相乘，这时相当于每一位都和相对的每一位相乘，然后将结果储存到应在的位置。前面的FFT复杂度是O(n logn)，最后才进行进位的操作，复杂度是O(n)。因此，总的复杂度就从O(n^2)降到了O(n logn)了。

#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>

const int maxn = 1 << 16;
const double pi = acos((double)-1);

struct virt{ // 定义虚数，并重载虚数的运算符
    double r;
    double i;
    void ins(double a = 0.0, double b = 0.0){
        r = a;
        i = b;
    }
    virt operator + (const virt &x){
        virt ret;
        ret.ins(r + x.r, i + x.i);
        return ret;
    }
    virt operator - (const virt &x){
        virt ret;
        ret.ins(r - x.r, i - x.i);
        return ret;
    }
    virt operator * (const virt &x){
        virt ret;
        ret.ins(r * x.r - i * x.i, r * x.i + i * x.r);
        return ret;
    }
}revtmp[maxn << 1], val1[maxn << 1], val2[maxn << 1];
int ep[21], ans[maxn << 1];
double base[19];

void pre(){
    ep[0] = 1;
    for (int i = 1; i <= 20; i++){
        ep[i] = ep[i - 1] << 1;
    }
    for (int i = 0; i < 19; i++){
        base[i] = 2 * pi / ep[i + 1];
    }
}

int rev(int x, int l){
    int ret = 0;

    while (l--){
        ret |= (x & 1);
        ret <<= 1;
        x >>= 1;
    }

    return ret >> 1;
} // 将数x在二进制形式的最后l位倒置

void reverse(virt *a, int len){
    int bits;

    bits = (int)ceil(log((double)len) / log(2.0));
    for (int i = 0; i < len; i++){
        revtmp[rev(i, bits)] = a[i];
    }
    for (int i = 0; i < len; i++){
        a[i] = revtmp[i];
    }
} // 重组虚数数组，以便之后分治处理DFT，假设长度为8，重置后按自然序排列的倒位序相应为0 4 2 6 1 5 3 7

void fft(virt *a, int len, bool idft){
    virt u, t;
    double bits = (int)ceil(log((double)len) / log(2.0));

    reverse(a, len); // 先将数倒位排序
    for (int i = 0; i < bits; i++){
        virt wi, w; // wi是单位旋转因子

        wi.ins(cos((idft ? -1 : 1) * base[i]), sin((idft ? -1 : 1) * base[i])); // 判断是dft还是idft，两种操作相逆
        for (int k = 0; k < len; k += ep[i + 1]){
            w.ins(1.0, 0.0); // 重置旋转因子
            for (int j = 0; j < ep[i]; j++){
/************下面几行是蝴蝶操作的几个步骤************/
                t = w * a[k + j + ep[i]];
                u = a[k + j];
                a[k + j] = u + t;
                a[k + j + ep[i]] = u - t;
                w = w * wi;
            }
        }
    }
    if (idft){ // 如果是idft，要将每一位除以多项式的长度
        for (int i = 0; i < len; i++){
            a[i].r /= len;
        }
    }
}

void cal(virt *a, virt *b, int len){
    fft(a, len, false);
    fft(b, len, false);
    for (int i = 0; i < len; i++){
        a[i] = a[i] * b[i];
    }
    fft(a, len, true);
#ifndef ONLINE_JUDGE
    puts("answer:");
    for (int i = 0; i < len; i++){
        printf("%2d : %10.5f %10.5f\n", i, a[i].r, a[i].i);
    }
    puts("");
#endif
}

void deal(char *a, char *b){
    int len_a = strlen(a);
    int len_b = strlen(b);
    int len = 1;

    if (len_a > len_b) len = (int)ceil(log((double)len_a) / log(2.0));
    else len = (int)ceil(log((double)len_b) / log(2.0));
    len = ep[len + 1];
#ifndef ONLINE_JUDGE
    printf("len %d\n", len);
#endif
    for (int i = 0; i < len_a; i++)
        val1[i].ins((double)a[len_a - i - 1] - '0');
    for (int i = len_a; i < len; i++)
        val1[i].ins();
    for (int i = 0; i < len_b; i++)
        val2[i].ins((double)b[len_b - i - 1] - '0');
    for (int i = len_b; i < len; i++)
        val2[i].ins();

    cal(val1, val2, len);
    for (int i = 0; i < len; i++){
        ans[i] = val1[i].r + 0.5;
    }
    for (int i = 1; i < len; i++){
        ans[i] += ans[i - 1] / 10;
        ans[i - 1] %= 10;
    }

    len = len_a + len_b;
    while (ans[len] <= 0) len--;
    while (len >= 0){
        printf("%d", ans[len]);
        len--;
    }
    puts("");
}

char in_a[maxn], in_b[maxn];

int main(){
    pre();
    while (~scanf("%s%s", in_a, in_b)){
        if (in_a[0] == '0' || in_b[0] == '0'){
            printf("0\n");
            continue;
        }
        deal(in_a, in_b);
    }

    return 0;
}

——written by Lyon