快速傅里叶变换

今天学习了一下算法导论上的快速傅里叶变换，并写了个模板，留作今后使用吧~~~

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <complex>
#include <iostream>
#include <algorithm>
using namespace std;
typedef long long Int64;
const int MAXN = 200010;
const double PI = acos(-1.0);
typedef complex<double> Complex;
char str1[MAXN], str2[MAXN];
Complex A[MAXN], B[MAXN];
int str_len1, str_len2, len;
int ans[MAXN];

void bulid(Complex _P[], Complex P[], int n, int m, int curr, int &cnt)
{
    if(m == n)
    {
        _P[curr] = P[cnt++];
    }
    else
    {
        bulid(_P, P, n, m * 2, curr, cnt);
        bulid(_P, P, n, m * 2, curr + m, cnt);
    }
    return ;
}

/*
*P 需要进行DFT变换的数据
*n 数据长度(2的整数次幂)
*oper 1(正变换), -1(负变换)
*/
void FFT(Complex P[], int n, int oper)
{
    static Complex _P[MAXN];
    int cnt = 0;
    bulid(_P, P, n, 1, 0, cnt);
    copy(_P, _P + n, P);
    for(int d = 0; (1 << d) < n; d++)
    {
        int m = 1 << d;
        int m2 = m * 2;
        double p0 = PI / m * oper;
        Complex unit_p0 = Complex(cos(p0), sin(p0));
        for(int i = 0; i < n; i += m2)
        {
            Complex unit = 1;
            for(int j = 0; j < m; ++j)
            {
                Complex &P1 = P[i + j + m], &P2 = P[i + j];
                Complex t = unit * P1;
                P1 = P2 - t;
                P2 = P2 + t;
                unit = unit * unit_p0;
            }
        }
    }
    if(oper == -1)
        for(int i = 0; i < n; ++i)
            P[i] = Complex(P[i].real()/n, P[i].imag());
    return ;
}


int main()
{
    //freopen("aa.in", "r", stdin);
    //freopen("bb.out", "w", stdout);

    while(scanf("%s %s", str1, str2) != EOF)
    {
        memset(ans, 0, sizeof(ans));
        str_len1 = strlen(str1);
        str_len2 = strlen(str2);
        int len = 1;
        while(len < (2*str_len1) || len < (2*str_len2))
            len <<= 1;
        for(int i = 0; i < str_len1; ++i)
        {
            A[i] = Complex(static_cast<double>(str1[str_len1-i-1]-'0'), 0.0);
        }
        for(int i = str_len1; i < len; ++i)
        {
            A[i] = Complex(0.0, 0.0);
        }
        for(int i = 0; i < str_len2; ++i)
        {
            B[i] = Complex(static_cast<double>(str2[str_len2-i-1]-'0'), 0.0);
        }
        for(int i = str_len2; i < len; ++i)
        {
            B[i] = Complex(0.0, 0.0);
        }
        FFT(A, len, 1); FFT(B, len, 1);
        for(int i = 0; i < len; ++i)
        {
            A[i] = A[i] * B[i];
        }
        FFT(A, len, -1);
        for(int i = 0; i < len; ++i)
        {
            ans[i] = static_cast<int>(A[i].real() + 0.5);
        }
        for(int i = 0; i < len; ++i)
        {
            ans[i + 1] += ans[i] / 10;
            ans[i] %= 10;
        }
        len = str_len1 + str_len2;
        while(!ans[len] && len > 0)
            len--;
        for(int i = len; i >= 0; --i)
            putchar((ans[i]+'0'));
        putchar('\n');

    }
    return 0;
}

/*自定义的Complex类，只需要重载相应的运算符即可*/

#include <cmath>
#include <cstdio>
#include <cstdlib>
#include <cstring>
#include <complex>
#include <iostream>
#include <algorithm>
using namespace std;
const int MAXN = 200010;
const double PI = acos(-1.0);

struct Complex
{
    double r, i;
    Complex(double _r = 0.0, double _i = 0.0)
    {
        this->r = _r;
        this->i = _i;
    }
    Complex operator+ (const Complex &b)
    {
        return Complex(r + b.r, i + b.i);
    }
    Complex operator- (const Complex &b)
    {
        return Complex(r - b.r, i - b.i);
    }
    Complex operator* (const Complex &b)
    {
        return Complex(r * b.r - i * b.i, r * b.i + i * b.r);
    }
};


/*
 * 进行FFT和IFFT前的反转变换。
 * 位置i和 （i二进制反转后位置）互换
 * len必须去2的幂
*/
void change(Complex y[], int len)
{
    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;
    }
}


/*
    *FTT求解
    *len必须为2^k形式
    *on==1时DFT，on==-1时是DFT
*/
void FFT(Complex y[], int len, int on)
{
    change(y, len);
    for(int h = 2; h <= len; h <<= 1)
    {
        Complex wn(cos(-on*2*PI/h), sin(-on*2*PI/h));
        for(int j = 0; j < len; j += h)
        {
            Complex w(1, 0);
            for(int k = j; k < j + h/2; ++k)
            {
                Complex u = y[k];
                Complex t = w * y[k+h/2];
                y[k] = u + t;
                y[k+h/2] = u - t;
                w = w * wn;
            }
        }
    }
    if(on == -1)
        for(int i = 0; i < len; ++i)
            y[i].r /= len;

    return ;
}

Complex x1[MAXN], x2[MAXN];
char str1[MAXN], str2[MAXN];
int len1, len2, len;
int ans[MAXN];

int main()
{
    //freopen("aa.in", "r", stdin);
    //freopen("bb.out", "w", stdout);

    while(scanf("%s %s", str1, str2) != EOF)
    {
        memset(ans, 0, sizeof(ans));
        len1 = strlen(str1), len2 = strlen(str2);
        len = 1;
        while(len < 2 * len1 || len < 2 * len2)
            len <<= 1;
        for(int i = 0; i < len1; ++i)
            x1[i] = Complex((str1[len1-i-1]-'0'), 0.0);
        for(int i = len1; i < len; ++i)
            x1[i] = Complex(0.0, 0.0);
        for(int i = 0; i < len2; ++i)
            x2[i] = Complex((str2[len2-i-1]-'0'), 0.0);
        for(int i = len2; i < len; ++i)
            x2[i] = Complex(0.0, 0.0);
        FFT(x1, len, 1); FFT(x2, len, 1);
        for(int i = 0; i < len; ++i)
            x1[i] = x1[i] * x2[i];
        FFT(x1, len, -1);
        for(int i = 0; i < len; ++i)
        {
            //cout << "Test: " << x1[i].r << endl;
            ans[i] = (int)(x1[i].r + 0.5);
        }
        for(int i = 0; i < len; ++i)
        {
            ans[i+1] += ans[i] / 10;
            ans[i] %= 10;
        }
        len = len1 + len2 - 1;
        while(ans[len] <= 0 && len > 0)
            len--;
        for(int i = len; i >= 0; --i)
            putchar((ans[i]+'0'));
        putchar('\n');
    }
    return 0;
}