快速傅里叶变换是离散傅氏变换的快速算法，它是根据离散傅氏变换的奇、偶、虚、实等特性，对离散傅立叶变换的算法进行改进获得的。

多项式

我们将一个次数为 $n-1$ 的多项式 $A(x)$ 表示为：

$$ A(x)=\sum_{i=0}^{n-1} a_ix^i $$

系数表示

多项式 $A(x)$ 的系数组成的 $n$ 维列向量（由于篇幅限制，本处写成行向量形式）$\vec a=(a_0,a_1,\dots,a_{n-1})$ 就是它的系数表示。

点值表示

我们将 $n$ 个互不相同的 $x$ 代入多项式得到对应的 $n$ 个值，这 $n$ 个二元组 $(x_i,y_i)$ 可以唯一确定这个多项式。我们将点值表示记为：

$$ \{(x_0,y_0),(x_1,y_1),\dots,(x_{n-1},y_{n-1})\} $$

求值与插值

我们有定理：一个 $n-1$ 次多项式可以由 $n$ 个不同点的取值唯一确定。

可以在 $\mathcal O(n^2)$ 时间内求出一个多项式的点值表达；已知点值表达也可以在 $\mathcal O(n^2)$ 时间内通过插值求出其系数表示。

多项式乘法

对于多项式的加减法这里不再赘述，考虑将两个 $n-1$次的多项式 $A(x)$ 和 $B(x)$ 相乘。其中 $A(x)=\sum_{i=0}^{n-1}a_ix^i,B(x)=\sum_{i=0}^{n-1}b_ix^i$。

系数乘法

多项式的加减法就不说了，考虑将两个 $n-1$ 次的多项式 $A(x)$ 和 $B(x)$ 相乘，那么得到的多项式 $C(x)$ 为：

$$ C(x)=A(x)\times B(x)=\sum_{k=0}^{2n-2}\left(\sum_{k=i+j}a_ib_j\right)x^k $$

显然这个复杂度为 $\mathcal O(n^2)$。

点值乘法

设多项式 $A(x),B(x)$ 的在第 $i$ 个点的点值表示为 $(x_{A_i},y_{A_i}),(x_{B_i},y_{B_i})$，那么有：

$$ y_{C_i}=y_{A_i}\times y_{B_i} $$

此时复杂度就优秀很多了，只需要 $\mathcal O(n)$。

复数

定义

设实数 $a,b$，定义虚数单位 $i$，有 $i^2=-1$，将形如 $a+bi$ 的数叫做复数。复数域是目前已知最大的域。

复平面

在复平面内，$x$ 轴表示实数，$y$ 轴表示虚数。从 $(0,0)$ 到 $(a,b)$ 的向量表示复数 $a+bi​$。

模长

复数 $z=a+bi$ 的模定义为其向量的长度 $\sqrt{a^2+b^2}$，记为 $\vert z\vert​$。

幅角

该表示从 $x$ 轴正半轴到该向量的转角的有向角（逆时针为正方向）叫做幅角。

因此复数也可以用 $(a,\theta)$ 表示，其中 $a,\theta$ 分别为模长、幅角。

共轭复数

复数 $z=a+bi$ 的共轭复数 $\overline z$ 和 $z$ 的乘积为实数，显然有 $\overline z=a-bi$（虚部符号取反）。

运算法则

加法

复数相加遵循平行四边形定则，设复数 $z_1=a+bi,z_2=c+di$，则有：

$$ z_1+z_2=(a+c)+(b+d)i $$

乘法

复数相乘有以下两种定义：

几何定义：设 $z_1=(a_1,\theta_1),z_=(a_2,\theta_2)$，则有 $z_1\times z_2=(a_1\times a_2,\theta_1+\theta_2)$。模长相乘，幅角相加。

代数定义：设 $z_1=a+bi,z_2=c+di$，则有 $z_1\times z_2=(ac-bd)+(bc+ad)i$。

单位根

下文若没有特别声明，默认 $n$ 为 $2$ 的正整数次幂。

在复平面上，以原点为圆心作单位圆（半径为 $1$ 的圆）。以原点为起点，将单位圆的 $n$ 等分点为终点，作出 $n$ 个向量。设其中幅角为正且最小的向量对应的复数为 $\omega_n$，将其称为 $n$ 次单位根。

根据复数乘法的几何定义，我们可以得到其余的 $n-1$ 个向量对应的复数分别为 $\omega_n^2,\omega_n^3,\dots,\omega_n^n$，其中 $\omega_n^n=\omega_n^0=1$。

由于单位根的幅角为圆周的 $\frac{1}{n}$，那么我们得到一个计算单位根及其幂的公式（这里融合的欧拉公式 $e^{ix}=\cos x+i\sin x$）：

$$ e^{k\frac{2\pi i}{n}}=\omega_n^k=\cos k\frac{2\pi}{n}+i\sin k\frac{2\pi}{n} $$

性质

消去引理

$$ \omega_{dn}^{dk}=\omega_n^k\quad (n,k,d\ge 0) $$

几何证明

在复平面上两者表示的向量终点相同。

代数证明

直接由单位根的定义式即可推导：

$$ \omega_{dn}^{dk}=e^{dk\frac{2\pi i}{dn}}=e^{k\frac{2\pi i}{n}}=\omega_n^k $$

折半引理

$$ \omega_n^{k+\frac{n}{2}}=-\omega_n^k\quad (n>0,n\ \text{为偶数}) $$

几何证明

这两个向量关于原点对称，因此为相反数关系。

代数证明

$$ \omega_n^{\frac{n}{2}}=e^{\pi i}=-1 $$

求和引理

这个引理是后文使用 $\text{FFT}$ 进行快速插值的基础：

$$ \sum_{j=0}^{n-1}\left(\omega_n^k\right)^j=0\quad (n\ge 1,n\not\mid k) $$

根据等比数列求和的知识可以得到：

$$ \sum_{j=0}^{n-1}\left(\omega_n^k\right)^j=\frac{\left(\omega_n^k\right)^n-1}{\omega_n^k-1}=\frac{1^k-1}{\omega_n^k-1}=0 $$

快速傅里叶变换 FFT

讲完了这么多前置知识，终于到重点啦！

考虑之前的多项式系数表示和点值表示，发现点值表示的乘法复杂度只有 $\mathcal O(n)$，是不是可以从这里入手呢？因此快速傅里叶变换的本质就是把多项式变成点值表示。

我们需要一些好的 $x$ 使得其若干次幂为 $\pm 1$ 或者 $\pm i$。这样会使得计算非常方便，这刚好运用了单位根良好的性质。我们把 $n$ 次单位根的 $0$ 到 $n-1$ 次幂代入多项式的系数表示，得到的点值向量 $(A(\omega_n^0),A(\omega_n^1),\dots,A(\omega_n^{n-1}))$ 就称为其系数向量 $(a_0,a_1,\dots,a_{n-1})$ 的离散傅里叶变换（$\text{DFT}$）。

但是这个过程按照朴素的变换还是 $\mathcal O(n^2)$ 的。

我们考虑将多项式按照系数下标的奇偶性分成 $2$ 部分：

$$ A(x)=(a_0+a_2x^2+\cdots+a_{n-2}x^{n-2})+(a_1x+a_3x^3+\cdots+a_{n-1}x^{n-1}) $$

右边提出一个 $x$，令：

$$ A_1(x)=a_0+a_2x+a_4x^2+\cdots+a_{n-2}x^{\frac{n}{2}-1} \\ A_2(x)=a_1+a_3x+a_5x^2+\cdots+a_{n-1}x^{\frac{n}{2}-1} $$

则有：

$$ A(x)=A_1(x^2)+xA_2(x^2) $$

假设 $k<\frac{n}{2}​$，现在我们将 $\omega_n^k​$ 代入 $A(x)​$ 得到：

$$ \begin{aligned} A(\omega_n^k)&=A_1(\omega_n^{2k})+\omega_n^k A_2(\omega_n^{2k}) \\ &=A_1(\omega_{\frac{n}{2}}^k)+\omega_n^k A_2(\omega_{\frac{n}{2}}^k) \end{aligned} $$

接下来，我们再将 $\omega_{n}^{k+\frac{n}{2}}$ 代入 $A(x)$ 得到：

$$ \begin{aligned} A(\omega_n^{k+\frac{n}{2}})&=A_1(\omega_n^{2k+n})+\omega_n^{k+\frac{n}{2}} A_2(\omega_n^{2k+n}) \\ &=A_1(\omega_{n}^{2k}\times\omega_n^n)-\omega_n^k A_2(\omega_{n}^{2k}\times\omega_n^n) \\ &=A_1(\omega_{n}^{2k})-\omega_n^k A_2(\omega_{n}^{2k}) \\ &=A_1(\omega_{\frac{n}{2}}^{k})-\omega_n^k A_2(\omega_{\frac{n}{2}}^{k}) \end{aligned} $$

注意到这这两个式子右侧的第二项只有符号不同；又注意到：当 $k$ 取遍 $[0,\frac{n}{2})$ 时，$k+\frac{n}{2}$ 取遍 $[\frac{n}{2},n)$，这意味着如果已知 $A_1(x),A_2(x)$ 在 $\omega_{\frac{n}{2}}^0,\omega_{\frac{n}{2}}^1,\cdots,\omega_{\frac{n}{2}}^{\frac{n}{2}-1}$ 处的点值，就可以在 $\mathcal O(n)$ 的值见内求出 $A(x)$ 在 $\omega_n^0,\omega_n^1,\cdots,\omega_n^{n-1}$ 处的点值。而关于 $A_1(x)$ 和 $A_2(x)$ 的问题都是相对于原问题规模缩小了一半的子问题，分治的边界为一个常数项 $a_0$。

根据主定理可以计算出这个分治算法的复杂度：

$$ T(n)=2T\left(\frac{n}{2}\right)+\mathcal O(n)=\mathcal O(n\log n) $$

快速傅里叶逆变换 IFFT

以上的 $\text{FFT}$ 只不过是求出了点值表示，而我们的最终目的是求这个多项式的系数表示。因此快速傅里叶逆变换的本质就是将点值表示重新转化为系数表示。这个过程就是傅里叶逆变换。

我们设 $(y_0,y_1,\cdots,y_{n-1})$ 为 $(a_0,a_1,\cdots,a_{n-1})$ 的傅里叶变换。考虑另一个向量 $(c_0,c_1,\cdots,c_{n-1})$ 满足：

$$ c_k=\sum_{i=0}^{n-1}y_i\left(w_n^{-k}\right)^i $$

即多项式 $P(x)=y_0+y_1x+y_2x^2+\cdots+y_{n-1}x^{n-1}$ 在 $\omega_{n}^0,\omega_n^{-1},\omega_n^{-2},\cdots,\omega_n^{-(n-1)}$ 处的点值表示。

将上式展开，得到：

$$ \begin{aligned} c_k&=\sum_{i=0}^{n-1} y_i\left(\omega_n^{-k}\right)^i \\ &=\sum_{i=0}^{n-1}\left(\sum_{j=0}^{n-1}a_j\left(\omega_n^i\right)^j\right)\left(\omega_n^{-k}\right)^i \\ &=\sum_{i=0}^{n-1}\sum_{j=0}^{n-1}a_j\left(\omega_{n}^{j-k}\right)^i \\ &=\sum_{j=0}^{n-1}a_j\sum_{i=0}^{n-1}\left(\omega_{n}^{j-k}\right)^i \end{aligned} $$

我们考虑一个式子：

$$ S(\omega_n^k)=1+\omega_n^k+\left(\omega_n^k\right)^2+\cdots+\left(\omega_n^k\right)^{n-1} $$

根据单位根的求和引理，我们可以得到：

$$ S(\omega_n^k)=0 $$

那么我们把上式转化为：

$$ c_k=\sum_{j=0}^{n-1}a_jS(\omega_n^{j-k}) $$

根据求和引理可以得到 $S(\omega_n^{j-k})$ 不等于 $0$ 当且仅当 $j=k$，此时 $S(\omega_n^{j-k})=n$。因此我们得到：

$$ \begin{aligned} c_i&=na_i \\ a_i&=\frac{1}{n}c_i \end{aligned} $$

所以我们用单位根的倒数代替单位根做一次类似 $\text{FFT}$ 的过程，再将结果的每个数都除以 $n$，就可以得到系数表示了！

实现

我们可以使用 $\text{STL}$ 中的 complex，但是为了减小常数建议手写复数类！

但是为什么没有 O2 的 $\text{STL}$ 比开 O2 的手写类快啊 QAQ

递归版

我们直接按照上面的分析来实现就行了。

代码

#include <cstdio>
#include <cmath>
#include <algorithm>

const int N = 3e6 + 5;
const double pi = acos(-1);

int n1, n2, a1[N], a2[N], ret[N];

struct Complex {
    double real, imag;
    Complex(double _real = 0,double _imag = 0) {
        real = _real, imag = _imag;
    }
    double len() const {
        return real * real + imag * imag;
    }
    Complex operator!() const {
        return Complex(real, -imag);
    }
    Complex operator+(const Complex &b) const {
        return Complex(real + b.real, imag + b.imag);
    }
    Complex operator-(const Complex &b) const {
        return Complex(real - b.real, imag - b.imag);
    }
    Complex operator*(const double &b) const {
        return Complex(real * b, imag * b);
    }
    Complex operator*(const Complex &b) const {
        return Complex(real * b.real - imag * b.imag, real * b.imag + imag * b.real);
    }
    Complex operator/(const double &b) const {
        return Complex(real / b, imag / b);
    }
    Complex operator/(const Complex &b) const {
        return *this * !b / b.len();
    }
} c1[N], c2[N];

struct FFT {
    int extend(int len) {
        int n = 1;
        for (; n < len; n <<= 1);
        return n;
    }
    Complex omega(int n, int k) {
        return Complex(cos(2 * pi * k / n), sin(2 * pi * k / n));
    }
    void trans(Complex *a, int n, bool inv) {
        if (n == 1) return;
        int m = n >> 1;
        Complex a1[m], a2[m];
        for (int i = 0; i < m; i++) {
            a1[i] = a[i << 1], a2[i] = a[i << 1 | 1];
        }
        trans(a1, m, inv);
        trans(a2, m, inv);
        for (int i = 0; i < m; i++) {
            Complex x = omega(n, i);
            if (inv) x = !x;
            a[i] = a1[i] + x * a2[i];
            a[i + m] = a1[i] - x * a2[i];
        }
    }
    void dft(Complex *a, int n) {
        trans(a, n, 0);
    }
    void idft(Complex *a, int n) {
        trans(a, n, 1);
        for (int i = 0; i < n; i++) {
            a[i] = a[i] / n;
        }
    }
} fft;

void mul(int *a1, int n1, int *a2, int n2, int *ret) {
    int n = fft.extend(n1 + n2 - 1);
    std::fill(c1, c1 + n, Complex(0, 0));
    std::fill(c2, c2 + n, Complex(0, 0));
    for (int i = 0; i < n1; i++) {
        c1[i].real = a1[i];
    }
    for (int i = 0; i < n2; i++) {
        c2[i].real = a2[i];
    }
    fft.dft(c1, n), fft.dft(c2, n);
    for (int i = 0; i < n; i++) {
        c1[i] = c1[i] * c2[i];
    }
    fft.idft(c1, n);
    for (int i = 0; i < n; i++) {
        ret[i] = (int)(c1[i].real + 0.5);
    }
}
int main() {
    scanf("%d%d", &n1, &n2);
    for (int i = 0; i < n1; i++) {
        scanf("%d", &a1[i]);
    }
    for (int i = 0; i < n2; i++) {
        scanf("%d", &a2[i]);
    }
    mul(a1, n1, a2, n2, ret);
    int n = n1 + n2 - 1;
    for (int i = 0; i < n; i++) {
        printf("%d%c", ret[i], " \n"[i == n - 1]);
    }
    return 0;
}

迭代版

分治边界顺序

然而递归的常数实在是太大了，需要优化这个过程。我们考虑递归 $\text{FFT}​$ 分治到边界时每个数的顺序及其二进制位。

000 001 010 011 100 101 110 111
 0   1   2   3   4   5   6   7
 0   2   4   6 | 1   3   5   7
 0   4 | 2   6 | 1   5 | 3   7
 0 | 4 | 2 | 6 | 1 | 5 | 3 | 7
000 100 010 110 001 101 011 111

发现规律：分治到边界后的下标等于原下标的二进制位翻转。

蝴蝶操作

考虑将两个子问题合并的过程，假如 $A_1(\omega_{\frac{n}{2}}^k)$ 和 $A_2(\omega_{\frac{n}{2}}^k)$ 分别存在 $a(k)$ 和 $a(k+\frac{n}{2})$ 中，$A(\omega_n^k)$ 和 $A(\omega_{n}^{k+\frac{n}{2}})$ 将要被存在 $b(k)$ 和 $b(k+\frac{n}{2})$ 中，那么合并的操作可以表示为：

$$ \begin{aligned} b(k)&\leftarrow a(k)+\omega_n^k\times a(k+\frac{n}{2}) \\ b(k+\frac{n}{2})&\leftarrow a(k)-\omega_n^k\times a(k+\frac{n}{2}) \end{aligned} $$

考虑加入临时变量使得该操作可以在原地完成，而不需要另外添加数组 $b$。换言之，我们把 $A(\omega_n^k)$ 和 $A(\omega_n^{k+\frac{n}{2}})$ 直接覆盖原来 $a(k)$ 和 $a(k+\frac{n}{2})$ 的值。

$$ \begin{aligned} x&\leftarrow a(k) \\ y&\leftarrow \omega_n^k\times a(k+\frac{n}{2}) \\ b(k)&\leftarrow x+y \\ b(k+\frac{n}{2})&\leftarrow x-y \end{aligned} $$

这里有一个小 $\text{trick}$ 可以使得不需要每次计算 $\omega$。如果要将长度为 $\frac{l}{2}$ 的序列答案合并成长度为 $l$ 的。其中二进制翻转的过程也很巧妙，细节详见代码。

代码

#include <cstdio>
#include <complex>

typedef std::complex<double> Complex;

const int N = 3e6 + 5;
const double pi = acos(-1);

int n1, n2, a1[N], a2[N], ret[N];
Complex c1[N], c2[N];

struct FFT {
    int rev[N];
    Complex omega(int n, int k) {
        return Complex(cos(2 * pi * k / n), sin(2 * pi * k / n));
    }
    int extend(int len) {
        int n = 1;
        for (; n < len; n <<= 1);
        return n;
    }
    void reverse(Complex *a, int n) {
        int k = 0;
        for (; (1 << k) < n; k++);
        for (int i = 0; i < n; i++) {
            rev[i] = rev[i >> 1] >> 1 | (i & 1) << (k - 1);
            if (i < rev[i]) std::swap(a[i], a[rev[i]]);
        }
    }
    void trans(Complex *a, int n, int inv) {
        reverse(a, n);
        for (int l = 2; l <= n; l <<= 1) {
            int m = l >> 1;
            Complex w = omega(l, 1);
            if (inv) w = conj(w);
            for (int j = 0; j < n; j += l) {
                Complex wk(1, 0);
                for (int i = 0; i < m; i++, wk *= w) {
                    Complex p = a[i + j], q = wk * a[i + j + m];
                    a[i + j] = p + q, a[i + j + m] = p - q;
                }
            }
        }
    }
    void DFT(Complex *a, int n) {
        trans(a, n, 0);
    }
    void IDFT(Complex *a, int n) {
        trans(a, n, 1);
        for (int i = 0; i < n; i++) a[i] /= n;
    }
} fft;

void mul(int *a1, int n1, int *a2, int n2, int *ret) {
    int n = fft.extend(n1 + n2 - 1);
    std::fill(c1, c1 + n, Complex(0, 0));
    std::fill(c2, c2 + n, Complex(0, 0));
    for (int i = 0; i < n1; i++) c1[i] = a1[i];
    for (int i = 0; i < n2; i++) c2[i] = a2[i];
    fft.DFT(c1, n), fft.DFT(c2, n);
    for (int i = 0; i < n; i++) c1[i] = c1[i] * c2[i];
    fft.IDFT(c1, n);
    for (int i = 0; i < n; i++) ret[i] = (int)(c1[i].real() + 0.5);
}
int main() {
    scanf("%d%d", &n1, &n2);
    for (int i = 0; i < n1; i++) scanf("%d", &a1[i]);
    for (int i = 0; i < n2; i++) scanf("%d", &a2[i]);
    mul(a1, n1, a2, n2, ret);
    int n = n1 + n2 - 1;
    for (int i = 0; i < n; i++) {
        printf("%d%c", ret[i], " \n"[i == n - 1]);
    }
    return 0;
}

总结

我们用下图来总结一下 $\text{FFT}$ 的本质：

习题

最后修改：2021 年 05 月 14 日
© 允许规范转载
如果觉得我的文章对你有用，请随意赞赏