【学习】快速傅里叶变化(FFT)

知识点

声明：$\text{以下内容来自于goodqt的课件}$

复数

形如 $a + bi$，其中 $a∈R$, $b ∈ R$, $i^2 = −1$ 的数称为复数，记作$z = a + bi ∈ C$。

$$(a + bi) + (c + di) = (a + c) + (b + d)i$$

$$(a + bi) − (c + di) = (a − c) + (b − d)i$$

$$(a + bi) × (c + di) = (ac − bd) + (bc + ad)i$$

$$\frac{a + bi}{c + di} =\frac{(a + bi)(c − di)}{(c + di)(c − di)}=\frac{ac + bd}{c^2 + d^2}+\frac{bc − ad}{c^2 + d^2}$$
在复数域 C 中定义如下运算：
$$e^z =\sum^{∞}_{n=0}\frac{z^n}{n!}$$

$$\cos z =\sum^{∞}_{n=0}(−1)^n\frac{z^{2n}}{(2n)!}$$

$$sin z = \sum^{∞}_{n=0}(−1)^n \frac{z^{2n+1}}{(2n + 1)!}$$
可以证明，当 z ∈ R 时，这里的定义和之前我们在实数域中给出的定义完全等价。

可以证明，在复数域中，下式依然成立：

$$e^{z_1+z_2} = e^{z_1} e^{z_2}$$

令 $z = ix$，立即得到：

$$e^{ix} = \cos x + i\sin x$$

对任意复数 $z = a + bi \neq 0$，必定存在唯一的 $0 ≤ θ < 2π$ 和唯一的 $r > 0$ 满足：

$$z = a + bi = r(\cos θ + i\sin θ) = re^{iθ}$$

其中 $|z| = r = a^2 + b^2$称作复数 $z$ 的模，$θ$ 称作复数 $z$ 的辐角主值。
把复数 $z = a + bi$ 看作平面直角坐标系中的点 $(a, b)$，$r$, $θ$ 便有了明确的几何意义。此平面叫做复平面。

注意到

$$z_1z_2 = r_1e^{iθ_1}r_2e^{iθ_2} = r_1r_2e^{ i(θ_1+θ_2)}$$

$$z_1z_2 = \frac{r_1e^{iθ_1}}{r_2e^{iθ_2}}=\frac{r_1}{r_2}e^{i(θ_1−θ_2)}$$

所以在复平面中复数加减按平行四边形法则，复数乘除为模乘除、辐角加减。

单位根

考虑方程

$$x^n = 1\ \ \ \ \ (n ∈ N^∗)$$

有 $n$ 个两两不同的解分别为

$$x_k = e^{\frac{2kπi}{n}}\ \ \ \ \ (k ∈ [0, n − 1] ∩ N)$$
几何意义非常明显，它们被称为 n 次单位根。

多项式

给出 n − 1 次多项式有多种方法，最常见的是给出每一项前面的系数，即：

$$f(x) = \sum^{n−1}_{i=0}a_ix^i$$

还有其他的方法，即给出 n 个两两不同的位置处 f(x) 的函数值。
显然如果要做乘法，第一个的复杂度需要 $O(n^2)$，而第二个的复杂度只需 $O(n)$。
快速傅立叶变换是用 $O(n \log n)$ 的时间把前者转换到后者、把后者转化到前者，这样我们就可以 $O(n \log n)$ 完成一次前者的乘法了。

快速傅立叶变换

我们首先要把 $n$ 上调至 $2$ 的幂，后面本来没有的系数都设置为 $0$ 就好了。
然后我们取“$n$ 个两两不同的位置”为 $n$ 个 $n$ 次单位根，记录$m = \frac{n}{2}$，则

$$f(e^{\frac{2k\pi i}{n}}) = \sum_{j=0}^{n-1}a_j(e^{\frac{2k\pi i}{n}})^j = \sum_{j=0}^{m-1}a_{2j}(e^{\frac{2k\pi i}{m}})^j + e^{\frac{k\pi i}{m}}\sum_{j=0}^{n-1}a_{2j+1}(e^{\frac{2k\pi i}{m}})^j$$

分别拆出奇数项和偶数项，直接递归即可。时间复杂度$O(n \log n)$。

快速傅立叶逆变换

常数

虽然复杂度没啥问题，但是递归版的快速傅立叶变换常数还是非常大的，所以可能需要优化常数。
考虑分治的过程，不停的拆分奇数项和偶数项，所以我们直接把标号全部倒置，从低位开始不停的合并 $0/1$ 就好了，非递归版的常数还是可以接受的。

代码详解

递归版

#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;

const double pi = acos(-1);
//定义复数类型
struct Com{
    double x, y;
    Com(double X = 0, double Y = 0){x = X, y = Y;}
    Com operator + (const Com &t) const{return Com(x+t.x, y+t.y);}
    Com operator - (const Com &t) const{return Com(x-t.x, y-t.y);}
    Com operator * (const Com &t) const{return Com(x*t.x - y*t.y, y*t.x + x*t.y);}
};
//A一开始是系数数组，在调用完函数后，会变成记录答案的数组，答案记录在实部
void fft(Com *A, int n, int f){
    //仅有一个数的时候直接返回，这时的系数恰好是答案
    if(n == 1) return; 
    int m = n/2;
    Com C1[m], C2[m];
    //分奇偶划分，递归
    for(int i = 0; i < m; i ++) C1[i] = A[i*2], C2[i] = A[i*2+1];
    fft(C1, m, f), fft(C2, m, f);
    //k等与0的时候的单位元，k等于1的时候的数值
    Com t1(1, 0), t2(cos(pi/m), f*sin(pi/m));
    //当i大于m时，$e^{i\pi}=-1$，所以后面的项顺便算出来
    for(int i = 0; i < m; i ++){
        A[i] = C1[i] + t1 * C2[i];
        A[i+m] = C1[i] - t1 * C2[i];
        //让k增大
        t1 = t1 * t2;
    }
}

Com a[300010], b[300010];

int main(){
    int n, m;
    scanf("%d%d", &n, &m);
    for(int i = 0; i <= n; i ++) scanf("%lf", &a[i].x);
    for(int i = 0; i <= m; i ++) scanf("%lf", &b[i].x);
    //为了容纳答案的长度
    m += n;
    //最多有n+m+1项，"n <= m" '=' 不能省略
    for(n = 1; n <= m; n <<= 1);
    fft(a, n, 1);
    fft(b, n, 1);
    for(int i = 0; i < n; i ++) a[i] = a[i] * b[i];
    fft(a, n, -1); 
    //有浮点误差，+0.5消除
    for(int i = 0; i <= m; i ++) printf("%d ", int((a[i].x/n+0.5)));
    return 0;
}

非递归版

//与上边一样的就不打了
#include <cstdio>
#include <cstring>
#include <iostream>
#include <algorithm>
#include <cmath>
using namespace std;

const int maxn = 300010;
const double pi = acos(-1);
int Bit, rx[maxn];

struct Com{
    double x, y;
    Com(double X = 0, double Y = 0){x = X, y = Y;}
    Com(const Com &t){x = t.x, y = t.y;}
    Com operator + (const Com &t) const{return Com(x+t.x, y+t.y);}
    Com operator - (const Com &t) const{return Com(x-t.x, y-t.y);}
    Com operator * (const Com &t) const{return Com(x*t.x - y*t.y, y*t.x + x*t.y);}
};

void fft(Com *A, int n, int f){
    //预处理二进制翻转之后的值是多少，只会执行一次
    //原理：100111 进行反转，我们得到 10011 的翻转 110010，右移一位，空出最高位，把之前的那个补上即可
    if(rx[n-1] != n-1) for(int i = 0; i < n; i ++) rx[i] = (rx[i>>1]>>1) | ((i&1) << (Bit-1));
    //进行位数交换，只在有特定条件的时候进行
    for(int i = 0; i < n; i ++) if(i < rx[i]) swap(A[i], A[rx[i]]);
    //长度 -> 开头 -> 每一项
    for(int i = 1; i < n; i <<= 1){
        const Com Base(cos(pi/i), f*sin(pi/i));
        for(int j = 0; j < n; j += i << 1){
            Com mi(1, 0);
            for(int k = 0; k < i; k ++, mi = mi * Base){
                Com x = A[j+k], y = A[i+j+k] * mi;
                //这一句等价于：当i大于m时，$e^{i\pi}=-1$，所以后面的项顺便算出来
                A[j+k] = x+y, A[i+j+k] = x-y;
            }
        }
    }
}

Com a[maxn], b[maxn];

int main(){
    int n, m;
    scanf("%d%d", &n, &m);
    for(int i = 0; i <= n; i ++) scanf("%lf", &a[i].x);
    for(int i = 0; i <= m; i ++) scanf("%lf", &b[i].x);
    m += n;
    for(n = 1; n <= m; n <<= 1, Bit ++);
    fft(a, n, 1);
    fft(b, n, 1);
    for(int i = 0; i <= n; i ++) a[i] = a[i] * b[i];
    fft(a, n, -1); 
    for(int i = 0; i <= m; i ++) printf("%d ", int((a[i].x/n+0.5)));
    puts("");
    return 0;
}

练习

1. BZOJ 2179 FFT快速傅立叶

2. BZOJ 2194 快速傅立叶之二

3. BZOJ 3527 [Zjoi2014]力