FFT

转载自： http://cubercsl.cn/note/FFT/

简介

快速傅里叶变换（FFT）

什么是FFT？FFT（fast Fourier transform）是
用来计算离散傅里叶变换（discrete Fourier transform，FFT）及其逆变换（IDFT）的快速算法。如果没有数学分析基础，这里很难用严密的语言把它解释清楚；但如果你只需要一个直观感受的话，不妨记住这句话：DFT把时域信号变换为频域信号。

注意，时域和频域只是两种信号分析方法，并不是两种不同类别的信号。在谈论语言信号的时候，“音量小”指的是时域，而“声音尖”指的是频域。在录音的时候，音频通常按照时域形式保存，即各个时刻采样到的振幅；如果需要各个频率的数据，则需要DFT。

DFT有一个很有意思的性质：时域卷积，频域卷积；频域卷积，时域卷积。如果你不知道什么是卷积（convolution），也没关系，请再记住一句话：多项式乘法实际上是多项式系数向量的卷积。

多项式乘法

给定两个单变量多项式A(x)A(x)和B(x)B(x)，次数均不超过nn，如果快速计算两者的乘积呢？最简单的方法就是把系数两两相乘，再相加。这样计算的时间复杂度为O(n2)O(n2)，当nn很大时速度将会很慢。注意，高精度整数的乘法是多项式乘法的特殊情况（想一想，为什么），所以多项式乘法的快速乘法也可以用来做高精度乘法。

上述算法之所以慢，是因为表示方法不好。虽然“系数序列”是最自然的表示方法，但却不适合用来计算乘法。在多项式快速乘法中，需要用点值法来表示一个多项式。点值表示是一个“点-值”对的序列{(x0,y0),(x1,y1)…(xn−1,yn−1)}{(x0,y0),(x1,y1)…(xn−1,yn−1)}，它代表一个满足yk=A(xk)yk=A(xk)的多项式AA。可以证明：恰好有一个次数小于nn的多项式满足这个条件。

点值表示法非常适合做乘法：只要两个多项式的点集{xi}{xi}相同，则只需要把对应的值乘起来就可以了，只需O(n)O(n)时间。但问题在于输入输出的时候仍需要采用传统的系数表示法，因此需要快速地在系数表示和点值表示之间转换。还记得刚才那句话吗？时域卷积，频域卷积。也就是说，多项式的系数表示法对应于时域，而点值法对应于频域，因此只需要用FFT计算出一个DFT就可以完成转换。

单位根

这样算出来的点值表示法，对应的求值点是哪些呢？答案是2n2n次单位根。所谓“nn次单位根”是满足xn=1xn=1 的复数。xn=1xn=1 的话，xx岂不是就等于1？很遗憾，那只是实数域中的情况。在复数域中，1恰好有nn个单位根，e2kπi/ne2kπi/n，其中ii是虚数单位(i2=−1)(i2=−1), 而eiu=cos u+i sin ueiu=cos u+i sin u。单位根非常特殊，因此FFT才有办法在更短的时间内求出多项式在这些点的取值。

利用FFT进行快速多项式乘法的具体步骤

1. （补0）在两个多项式的最前面补0，得到两个2n2n次多项式，设系数向量分别为v1v1和v2v2。
2. （求值）用FFT计算f1=DFT(v1)f1=DFT(v1)和f2=DFT(v2)f2=DFT(v2)。这里得到的f1f1和f2f2分别是两个输入多项式在2n2n次单位根出的各个取值（即点值表示）。
3. （乘法）把两个向量f1f1和f2f2对应的每一维对应相乘，得带向量ff。它对应输入多项式乘积的点值表示。
4. （插值）用FFT计算v=IDFT(f)v=IDFT(f)，其中vv就是乘积的系数向量。

例题

题目来源 SHUOJ 1966

Description

给出A0,A1,…An−1A0,A1,…An−1和B0,B1,…Bn−1B0,B1,…Bn−1，求C0,C1,…Cn−1C0,C1,…Cn−1
其中

Ci+j=∑AiBjCi+j=∑AiBj

Input

多组数据，第一行为一个数字nn。
之后有两行数列，分别表示A0,A1,…An−1A0,A1,…An−1和B0,B1,…Bn−1B0,B1,…Bn−1。
1≤n≤1051≤n≤105
1≤Ai,Bi≤1001≤Ai,Bi≤100

Output

对于每组数据，输出一行共nn个数字。

Sample Input

1 2 3

3 1 2 3 1 2 3

Sample Output

1

1 4 10

代码：

#include <bits/stdc++.h>
using namespace std;

const double PI = acos(-1.0);
//复数结构体
struct Complex
{
    double x, y; //实部和虚部 x+yi
    Complex(double tx = 0.0, double ty = 0.0)
    {
        x = tx;
        y = ty;
    }
    Complex operator-(const Complex &b) const
    {
        return Complex(x - b.x, y - b.y);
    }
    Complex operator+(const Complex &b) const
    {
        return Complex(x + b.x, y + b.y);
    }
    Complex operator*(const Complex &b) const
    {
        return Complex(x * b.x - y * b.y, x * b.y + y * b.x);
    }
};
/*
* 进行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]);
        //交换互为小标反转的元素，i<j保证交换一次
        //i做正常的+1，j左反转类型的+1,始终保持i和j是反转的
        k = len / 2;
        while (j >= k)
        {
            j -= k;
            k /= 2;
        }
        if (j < k)
            j += k;
    }
}
/*
* 做FFT
* len必须为2^k形式，
* on==1时是DFT，on==-1时是IDFT
*/
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].x /= len;
}

const int MAXN = 400010;
Complex a[MAXN], b[MAXN], c[MAXN];
int x[MAXN], y[MAXN];
int sum[MAXN];

int main()
{
    int n;
    while (cin >> n)
    {
        int t = 1;
        for (int i = 0; i < n; i++)
            cin >> x[i];
        for (int i = 0; i < n; i++)
            cin >> y[i];
        for (int i = 0; i < n; i++)
        {
            a[i] = Complex(x[i], 0);
            b[i] = Complex(y[i], 0);
        }
        while (t < n * 2)
            t <<= 1;
        for (int i = n; i < t; i++)
        {
            a[i] = Complex(0, 0);
            b[i] = Complex(0, 0);
        }
        fft(a, t, 1);
        fft(b, t, 1);
        for (int i = 0; i < t; i++)
            c[i] = a[i] * b[i];
        fft(c, t, -1);
        for (int i = 0; i < t; i++)
            sum[i] = (int)(c[i].x + 0.5);
        for (int i = 0; i < n; i++)
        {
            if (i)
                cout << " " << sum[i];
            else
                cout << sum[i];
        }
        cout << endl;
    }
    return 0;
}