Radix-2 迭代FFT

转载自http://blog.csdn.net/sxy_cnyali/article/details/52143276

#include <bits/stdc++.h>

using namespace std;

#define REP(i, a, b) for (int i = (a), _end_ = (b); i <= _end_; ++i)
#define debug(...) fprintf(stderr, __VA_ARGS__)
#define mp make_pair
#define x first
#define y second
#define pb push_back
#define SZ(x) (int((x).size()))
#define ALL(x) (x).begin(), (x).end()

template<typename T> inline bool chkmin(T &a, const T &b){ return a > b ? a = b, 1 : 0; } 
template<typename T> inline bool chkmax(T &a, const T &b){ return a < b ? a = b, 1 : 0; }

typedef long long LL;

const int dmax = 300100 << 1, oo = 0x3f3f3f3f;

const double PI = acos(-1);

int N, M, n;

complex<double> a[dmax], b[dmax], p[dmax];

int c[dmax], rev[dmax];

inline int read()
{
    int sum = 0, fg = 1, c = getchar();
    while (c < '0' || c > '9') {if (c == '-')  fg = -1; c = getchar(); }
    while (c >= '0' && c <= '9') sum = sum * 10 + c - '0', c = getchar();
    return sum * fg;
}

inline void init()
{
    int dis = 0;
    for (n = 1; n <= N + M; n <<= 1) ++dis;
    REP(i, 0, n)
        rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (dis - 1));
}

void FFT(complex<double> x[], const int &N, const int &type = 1)
{
    REP(i, 0, N)
        if (i < rev[i])
            swap(x[i], x[rev[i]]);
    for (int k = 1; k < N; k <<= 1)
    {
        complex<double> wn(cos(PI / k), sin(type * PI / k));
        for (int i = 0; i <= N; i += (k << 1))
        {
            complex<double> w(1, 0);
            for (int j = 0; j < k; ++j)
            {
                complex<double> X = x[i + j],
                                Y = x[i + j + k] * w;
                x[i + j] = X + Y;
                x[i + j + k] = X - Y;
                w *= wn;
            }
        }
    }
}

int main()
{
#ifndef ONLINE_JUDGE
    freopen("input.txt", "r", stdin);
    freopen("output.txt", "w", stdout);
#endif
    N = read(), M = read();
    REP(i, 0, N) a[i] = read();
    REP(i, 0, M) b[i] = read();
    init();
    FFT(a, n), FFT(b, n);
    REP(i, 0, n) p[i] = a[i] * b[i];
    FFT(p, n, -1);
    REP(i, 0, N + M) c[i] = p[i].real() / n + 0.1;
    REP(i, 0, N + M) printf("%d ", c[i]);
    puts("");
    return 0;
}