快速傅里叶变换模板

快速傅里叶变换模板

#include<iostream>
#include<cstdio>
#include<cmath>
#define endl "\n"
#define gc getchar()
#define rep(a,b,c) for (int i = a; i <= b; i += c)

using namespace std;

const int N = 1e7 + 10;
const double pi = acos(-1.0);

template <class T>
inline void read(T& x)
{
    x = (T)0;
    T tag = 1;
    char c = getchar();
    while (c < '0' || c > '9')
    {
        if (c == '-')
            tag = -1;
        c = getchar();
    }
    while (c >= '0' && c <= '9')
    {
        x = (x << 1) + (x << 3) + (c ^ 48);
        c = getchar();
    }
    x *= tag;
}

struct cp
{
    double x, y;
    cp(double xx = 0, double yy = 0) { x = xx, y = yy; }
    cp operator + (cp b) { return cp(x + b.x, y + b.y); }
    cp operator - (cp b) { return cp(x - b.x, y - b.y); }
    cp operator * (cp b) { return cp(x * b.x - y * b.y, x * b.y + y * b.x); } 
}a[N], b[N];

int n, m;
int len, r[N];
int limit = 1;

void fast_fast_tle(cp* A, int type)
{
    rep(0, limit - 1, 1)
        if (i < r[i]) swap(A[i], A[r[i]]);
    //迭代实现操作，从最小长度的序列自底向上求出要迭代的序列 
    for (int mid = 1; mid < limit; mid <<= 1)//待合并区间的中点
    {
        cp Wn(cos(pi / mid), type * sin(pi / mid)); //单位根 
        for (int R = mid << 1, j = 0; j < limit; j += R)//R是区间的右端点，j表示前已经到哪个位置了 
        {
            cp w(1, 0);//幂 
            for (int k = 0; k < mid; k++, w = w * Wn)//枚举左半部分 
            {
                cp x = A[j + k], y = w * A[j + mid + k];//蝴蝶效应 
                A[j + k] = x + y;
                A[j + mid + k] = x - y;
            }
        }
    }
}
int main()
{
    read(n), read(m);
    rep(0, n ,1) read(a[i].x);
    rep(0, m, 1) read(b[i].x);
    while (limit <= n + m) limit <<= 1, len++;
    //len表示最大的二进制长度
    cout << limit << " " << len << endl;
    //对index进行二进制翻转
    rep(0, limit - 1, 1)
        r[i] = (r[i >> 1] >> 1) | ((i & 1) << (len - 1));
    // 在原序列中 i 与 i/2 的关系是 ： i可以看做是i/2的二进制上的每一位左移一位得来
    // 那么在反转后的数组中就需要右移一位，同时特殊处理一下复数 
    rep(0, limit - 1, 1)
        cout << r[i] << " ";
    cout << endl;

    fast_fast_tle(a, 1);
    rep(0, limit - 1, 1)
        cout << r[i] << " ";

    fast_fast_tle(b, 1);
    rep(0, limit - 1, 1)
        cout << r[i] << " ";
    cout << endl;
    rep(0, limit, 1) a[i] = a[i] * b[i];
    fast_fast_tle(a, -1);
    rep(0, limit - 1, 1)
        cout << r[i] << " ";
    cout << endl;
    rep(0, n + m, 1)
        printf("%d ", (int)(a[i].x / limit + 0.5));
    return 0;
}