快速傅里叶变换(FFT)

快速傅里叶变换大学的时候就学了，可现在想起来都还回去了已经。前两天做到一道大数乘法题：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1028

需要用到傅里叶变换，就又从新拿出算法导论看了一遍，记下来加深一下记忆！！

傅里叶变换多用于信号处理，在频率域和时间域之间进行变换，我只知道这么多了！！

对于大数乘法，对每个大数，可以表示成一个多项式，如果直接进行多项式相乘，那么时间复杂度是O(n^2)的。那为什么用了傅里叶变换就可以做到O(n*lg(n))了呢？听我慢慢道来

一个n次多项式可以有两种表示方法：

1.系数表示法：即用系数向量来表示一个n次多项式(a0, a1, a2,......an-1)

2.点值法：用n个点值对来表示，{(x0, y0), (x1, y1),......(xn-1, yn-1)};要求这n个x不能相同，为什么呢?其实说白了，给你n个点值对，而且你知道多项式是一个n次多项式，那么你可以列出n个方程，解这个方程组，你就可以得出这个n次多项式的系数，之所以要求这n个x不能相同，就是为了保证你能够得到n个有效的方程，这里面会有一个范德蒙德矩阵，线性代数里学过，也忘记了，就记得这个矩阵的行列式可以直接套公式~~~

闲言少叙，点值法和系数法其实是一一对应的，如果给你两个点值法表示的多项式[两个多项式用的是同一组(x0, x1, ...xn)]，让你求他们点值法表示的乘积：很明显，直接对应相乘就好了,时间复杂度为O(n)；FFT之所以能将时间复杂度从O(n^2)降到O(n*lgn)，就是因为它找到了将点值法表示的多项式在O(n*lgn)时间内转化成对应的系数表示法！

整个过程是这样的：

1.给两个系数表示法表示的n次多项式 ===》》 变换成两个点值表示法表示的多项式：这一步涉及到了扩充，从n次扩充到了2n次，称之为求值

2.给两个点值表示法表示的多项式对应项相乘得到对应结果的点值法表示

3.将结果的点值法表示转化成系数法表示

那用什么方法可以在O(n*lgn)时间内完成点值法和系数法之间的转换呢？

这就是FFT

其中的思想主要是二分，具体方法我也没有弄明白，递归实现的FFT还可以看懂，迭代实现的FFT真心没有弄明白，仅仅那个位置换操作的写法我就看不懂了~~~~

先说从系数法到点值法的转换：

递归实现的FFT里是利用了单位复根的特殊性质：n个n次单位复根,把这n个复数平方,你会得到n/2个不同的数，每个出现两次；所以原问题就是这样被分成了更小的子问题~~~~T(n) = 2T(n/2) + O(n),其实最后一项不是O，而是西塔，你懂的，不好打出来的字符！所以整个时间复杂度就成了O(n*lgn)

具体是这样的：给你一个n次多项式的系数表示A(x) = a0 + a1*x^1 + a2*x^2......an-1*x^(n-1);

让你求他的点值表示，你需要求n个点值对，那么这时候你选择求这n个点对应的值，怎么求呢？

你把这个n次多项式的系数表示拆成两组：

一组是偶数项系数：A[0](x) = a0 + a2*x + a4*x^2 + .......a(n-2) * x^(n/2 - 1)

一组是奇数项系数：A[1](x) = a1 + a3*x + a5*x^2 + .......a(n-1) * x^(n/2 - 1)

这样原来的A(x) = A[0](x^2) + x * A[1](x^2)

所以原来你需要求A(x)在 这n个点的值,现在你只要求这n个点中的一半就可以了,因为这n个点的平方结果中只有n/2个不同的值.也就是说你只要求出A[0](x)和A[1](x)在点这n/2点处的值就可以了，具体为什么就不说了，一算就明白了；最后贴一个FFT的迭代实现，具体我也不是很明白，希望以后能弄懂~~~~~

// Copyright (c) 2013
// Author: Muye (muyepiaozhou@gmail.com)

#include <iostream>
#include <complex>
#include <string>
#include <cmath>
#include <algorithm>
using namespace std;

const double Pi = acos(-1);
const int LEN = 262145;
int ans[LEN] = {0};
complex<double> A[LEN], B[LEN], C[LEN];

int Init(string& sa, string& sb) {
    int la = sa.length();
    int lb = sb.length();
    int lc = 1;
    int l = max(la, lb) * 2 - 1;
    while(lc < l) lc <<= 1;
    for(int i = 0; i < la; ++i) A[i] = sa[la - i - 1] - '0';
    for(int i = la; i < lc; ++i) A[i] = 0;
    for(int i = 0; i < lb; ++i) B[i] = sb[lb - i - 1] - '0';
    for(int i = lb; i < lc; ++i) B[i] = 0;
    return lc;
}

void Bit_Reverse_Copy(complex<double> *p, int n) {
    for(int i = 1, j = n >> 1, k; i < n - 1; ++i) {
        if(i < j) swap(p[i], p[j]);
        k = n >> 1;
        while(j >= k) {
            j -= k;
            k >>= 1;
        }
        if(j < k) j += k;
    }
}

void Iterative_FFT(complex<double> *p, int n, int inverse = 1) {
    Bit_Reverse_Copy(p, n);
    for(int m = 2; m <= n; m <<= 1) {
        complex<double> wn(cos(inverse * 2 * Pi / m), sin(inverse * 2 * Pi / m));
        for(int k = 0; k < n; k += m) {
            complex<double> w(1, 0);
            for(int j = 0; j < m >> 1; ++j) {
                complex<double> t = w * p[k + j + m / 2];
                complex<double> u = p[k + j];
                p[k + j] = u + t;
                p[k + j + m / 2] = u - t;
                w = w * wn;
            }

        }
    }
    if(inverse == -1) {
        for(int i = 0; i < n; ++i) p[i] /= n;
    }
}

int main() {
    string sa, sb;
    cin >> sa >> sb;
    int n = Init(sa, sb);
    Iterative_FFT(A, n);
    Iterative_FFT(B, n);
    for(int i = 0; i < n; ++i) C[i] = A[i] * B[i];
    Iterative_FFT(C, n, -1);
    for(int i = 0; i < n; ++i) ans[i] = C[i].real() + 0.5;
    for(int i = 0; i < n; ++i) {
        ans[i+1] += ans[i] / 10;
        ans[i] %= 10;
    }
    while(!ans[--n] && n);
    while(n >= 0) cout << ans[n--];
    cout << endl;
    return 0;
}