【FFT-快速傅立叶变换】

题意：给定两个整数，整数的长度最多能达到50000位，输出两个整数的乘积。

分析：题意非常的明了，一个惊世骇俗的想法是使用两个数组将整数保留起来，然后模拟我们平常手算时的乘法，不过这样一来时间复杂度将是O(N^2)，由于N过大，因此此题因寻求更加快速的解法。

　　对于任何一个N位的整数都可以看作是An*10^(n-1) + An-1*10^(n-2) + ... + A2*10^2 + A1*10 + A0。如果把10看作是一个自变量，那么任何一个整数就可以视作为一个多项式，两个整数相乘也便可以看作是两个多项式相乘。对于一个多项式，我们平时所接触到的多是其系数表示法，普通的相乘也就是建立在两个整数均采用系数表示法的基础上进行的。那么要使得计算多项式相乘的复杂度下降的另一种方式就是寻找一种新的表示多项式的方法......

　　若一个多项式的最高阶位为N-1，那么取N个点对(xi, yi)就能够唯一确定这个多项式，可以想象成有N个系数需要N个方程去求解。那么在此就可以寻找点对来表示一个多项式，对于一个大的数，看作多项式之后，那么舍弃掉原来以10为自变量的取值，而选取其他值，再通过计算多项式An*xi^(n-1) + An-1*xi^(n-2) + ... + A2*xi^2 + A1*xi + A0来保存这个多项式的信息。需要选取N个xi形成N对(xi, yi)方可唯一确定原来各个项前的系数，通过选取1的N次单位复根即可，并且利用单位复根的性质，可以使得计算量下降。

　　通过点值法表示多项式后，计算乘法也就是O(N)的时间了，由于两个数相乘使得项数变多，因此需要在之前尽可能多取点。FFT算法能够在O(NlogN)时间内将系数法转化为点值法，相乘后再有点值法转为系数法，该题就是使用的这个方法。

　　顺便说下一FFT过程中，计算叶子DTF时采用的二进制平摊反转置换，其作用是为了避免算法的递归而实现自底向上的计算方式。回顾一下在计算原串DFT的时候，假设离散点数为0-7，那么有以下过程：

(0 1 2 3 4 5 7) = (0 2 4 6) + (1 3 5 7) 1
(0 2 4 6) = (0 4) + (2 6) 2
(1 3 5 7) = (1 3) + (5 7) 2
(0 4) = (0) + (4) 3
(2 6) = (2) + (6) 3
(1 3) = (1) + (3) 3
(5 7) = (5) + (7) 3

分析这些分组的二进制位会发现，第1次分组是根据第0位是否为1来划分的，即奇偶性；第2次分组是根据第1位是否为1来划分的；第三次分组是根据第2位是否为1来划分的。这个特性与与一般的按照大小划分的数很类似（首先按照最高为是否为1划分，然后是次高位...），因此就可以通过一个算法来使得00...00 - 11...11这样递增的序列中的每一个数实现高位和低位的翻转，二进制平摊反转置换就是用于达到这个目的。

算法从1开始到N-2（FFT算法要求N必须是2的幂，保证每次折半之后不会出现奇数），因此0和N-1翻转后还是本身，接着维护好一个下标 j ，这个数就是与递增中的第 i 个数翻转之后对应的数，初始化 j 的下标示 N/2，这个数要是从后往前来定义二进制数的话，就会是1，例如若N=8，那么4的二进制位为100，假定只有3位2进制位组成2进制数，从右往左看，其值为1。接下来就是要找 j 的下一个数了，这个数从右往左看应该是2才能够满足要求，于是从右往左寻找，遇到1变为0，知道遇到0就跳出，并且将该位赋值为1，这个伟大的过程的作用仅仅只是给在右往左定义的二进制数 j 加了一个1。

#include <cstdlib>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <iostream>
#include <algorithm>
using namespace std;

struct C {
    double r, i;
    C() {}
    C(double _r, double _i) : r(_r), i(_i) {}
    inline C operator + (const C & a) const {
        return C(r + a.r, i + a.i);
    }
    inline C operator - (const C & a) const {
        return C(r - a.r, i - a.i);
    }
    inline C operator * (const C & a) const {
        return C(r*a.r - i*a.i, r*a.i + i*a.r);
    }
};

typedef long long LL;
const double pi = acos(-1.0);
const int N = 50005;
C a[N<<2], b[N<<2];
char num1[N], num2[N];
LL ret[N<<2];

void brc(C *y, int L) {
    int i, j, k;
    for (i=1,j=L>>1; i<L-1; ++i) { // 二进制平摊反转置换 O(NlogN)
        if (i < j) swap(y[i], y[j]);
        k = L>>1;
        while (j >= k) {
            j -= k;
            k >>= 1;
        }
        j += k;
    }
}

void FFT(C *y, int L, int dir) {
    brc(y, L);
    for (int h = 2; h <= L; h <<= 1) { // 枚举所需计算的点数 
        C wn(cos(dir*2*pi/h), sin(dir*2*pi/h)); // h次单位复根 
        for (int j = 0; j < L; j += h) { // 原序列被分成了L/h段h长序列 
            C w(1, 0); // 旋转因子 
            for (int k = j; k < j+h/2; ++k) { // 因为折半定理，只需要计算枚举一半的长度即可 
                C u = y[k];
                C t = w*y[k+h/2];
                y[k] = u + t;
                y[k+h/2] = u - t;
                w = w * wn; // 更新旋转因子 
            }
        }
    }
    if (dir == 1) {
        for (int i = 0; i < L; ++i) {
            y[i] = y[i] * C(1.0/L, 0);
        }
    }
}

int main() {
    while (scanf("%s %s", num1, num2) != EOF) {
        memset(ret, 0, sizeof (ret));
        int len1 = strlen(num1), len2 = strlen(num2);
        int ML = len1+len2-1, L = 1;
        while (L < ML) L <<= 1;
        for (int i = len1-1, j = 0; i >= 0; --i, ++j) {
            a[j] = C(num1[i]-'0', 0);
        }
        for (int i = len2-1, j = 0; i >= 0; --i, ++j) {
            b[j] = C(num2[i]-'0', 0);
        }
        for (int i = len1; i < L; ++i) a[i] = C(0, 0);
        for (int i = len2; i < L; ++i) b[i] = C(0, 0);
        FFT(a, L, -1), FFT(b, L, -1);
        for (int i = 0; i < L; ++i) {
            a[i] = a[i] * b[i];
        }
        FFT(a, L, 1);
        for (int i = 0; i < L; ++i) {
            ret[i] = (LL)floor(a[i].r + 0.5);
        }
        for (int i = 0; i < L; ++i) {
            ret[i+1] += ret[i] / 10;
            ret[i] %= 10;
        }
        int p = L;
        while (!ret[p] && p) --p;
        while (p >= 0) printf("%d", (int)ret[p--]);
        puts("");
    }
    return 0;
}

View Code