发布日期: 2024-01-27

更新日期: 2024-02-04

文章字数: 3.8k

快速傅里叶变换 FFT

题意：

给定一个 $n$ 次多项式 $F(x)$，和一个 $m$ 次多项式 $G(x)$。

请求出 $F(x)$ 和 $G(x)$ 的加法卷积，即
$[F \cdot G](x) = \sum_{k=0}^{n+m}\sum_{i+j = k} f_ig_j\,x^k$
输入格式：

第一行两个整数 $n,m$。

接下来一行 $n+1$ 个数字，从低到高表示 $F(x)$ 的系数。

接下来一行 $m+1$ 个数字，从低到高表示 $G(x)$ 的系数。

输出格式：

一行 $n+m+1$ 个数字，从低到高表示 $F(x) \cdot G(x)$ 的系数。

数据范围：

保证输入中的系数大于等于 $0$ 且小于等于 $9$。

对于 $100\%$ 的数据：$1 \le n, m \leq {10}^6$。

快速傅里叶变换 (FFT) 是一种用于快速计算 DFT 的算法。

在此之前，我们需要先了解什么是 DFT 。不过理论上来讲，OI 只要记住 FFT 板子就行了。

傅里叶变换

本部分仅作了解即可，OI 完全不会涉及。

傅里叶变换（Fourier Transform）是一种分析信号的方法，它可分析信号的成分，也可用这些成分合成信号。许多波形可作为信号的成分，傅里叶变换用正弦波作为信号的成分。

设 $f(t)$ 是关于时间 $t$ 的函数，则傅里叶变换可以检测频率 $\omega$ 的周期在 $f(t)$ 出现的程度：

$F(\omega)=\mathbb{F}[f(t)]=\int_{-\infty}^{\infty} f(t) \mathrm{e}^{-\mathrm{i} \omega t} \mathrm{d} t$

它的逆变换为

$f(t)=\mathbb{F}^{-1}[F(\omega)]=\frac{1}{2 \pi} \int_{-\infty}^{\infty} F(\omega) \mathrm{e}^{\mathrm{i} \omega t} \mathrm{d} \omega$

逆变换的形式与正变换非常类似，分母 $2\pi$ 恰好是指数函数的周期。

傅里叶变换相当于将时域的函数与周期为 $2\pi$ 的复指数函数进行连续的内积，逆变换仍旧为一个内积。

傅里叶变换有相应的卷积定理，可以将时域的卷积转化为频域的乘积，同样也能逆转化。

离散傅里叶变换 DFT

离散傅里叶变换（Discrete Fourier transform, DFT）是傅里叶变换在时域和频域上都呈离散的形式，将信号的时域采样变换为其 DTFT（discrete-time Fourier transform）的频域采样。

也就是说，傅里叶变换是积分形式的连续的函数内积，离散傅里叶变换是求和形式的内积

设 $\left\{x_n\right\}_{n=0}^{N-1}$ 是某一满足有限性条件的序列，则它的离散傅里叶变换 (DFT) 为

$X_k=\sum_{n=0}^{N-1} x_n \mathrm{e}^{-\mathrm{i} \frac{2 \pi}{N} k n}$

通常以 $\mathcal{F}$ 表示这一变换，即 $\hat x = \mathcal{F}x$ 。

类似与积分形式，它的逆离散傅里叶变换 (IDFT) 为

$x_n=\frac{1}{N} \sum_{k=0}^{N-1} X_k \mathrm{e}^{\mathrm{i} \frac{2 \pi}{N} k n}$

记为 $x = \mathcal{F}^{-1}\hat{x}$ 。

实际上，DFT 和 IDFT 变换式中和式前面的归一化系数并不重要。

在上面的定义中，DFT 和 IDFT 前的系数分别为 $1$ 和 $\frac{1}{N}$ 。有时我们会将这两个系数都改为 $\frac{1}{\sqrt{N}}$ 。

离散傅里叶变换仍旧是时域到频域的变换。而其求和形式的特殊性，使其可以有其他的解释方法。

如果把序列 $x_n$ 看作多项式 $F(x)$ 的 $x^n$ 项系数。

则计算得到的 $X_k$ 恰好是多项式 $F(x)$ 代入单位根 $\mathrm{e}^{-\frac{2k \pi\mathrm{i}}{n}}$ 的点值 $F(\mathrm{e}^{-\frac{2k \pi\mathrm{i}}{n}})$ 。

这便构成了卷积定理的另一种解释方法，即对多项式进行特殊的插值操作。

补充知识：

插值 (Interpolation)，是一种通过已知的、离散的数据点，在范围内推求新数据点的过程或方法。

$n$ 次单位根是 $n$ 次幂为 $1$ 的复数，且 $n$ 次的单位根恰好有 $n$ 个，即 $z_k = \mathrm{e}^{\frac{2k\pi\mathrm{i}}{n}}\,(k=0,1,\cdots,n-1)$ 。

而 DFT 恰好是多项式在单位根处进行插值。

例如计算

$\binom{n}{3} + \binom{n}{7} + \binom{n}{11}+\binom{n}{15} + \cdots$

可以发现，这道题计算的组合数的 $k$ 都是模 $4$ 余 $3$ 的。

考虑定义函数 $F(x)$ 为

$F(x) = (1 + x)^n = \binom{n}{0}x^0 + \binom{n}{1}x^1 + \binom{n}{2}x^2+\cdots$

代入四次单位根 $\mathrm{i}$ 可得

$f(\mathrm{i}) = (1 + \mathrm{i})^{n} = \binom{n}{0} + \binom{n}{1}\mathrm{i}-\binom{n}{2}-\binom{n}{3}\mathrm{i}+\cdots$

于是下面的求和恰好可以把各项消掉

$f(1)+\mathrm{i} f(\mathrm{i})-f(-1)-\mathrm{i} f(-\mathrm{i}) =4\binom{n}{3}+4\binom{n}{7}+4\binom{n}{11}+4\binom{n}{15}+\ldots$

因此答案为

$\frac{2^n+\mathrm{i}(1+\mathrm{i})^n-\mathrm{i}(1-\mathrm{i})^n}{4}$

这道题在单位根处插值，恰好构成 DFT 。

矩阵公式

由于离散傅立叶变换是一个线性算子，所以它可以用矩阵乘法来描述。

在矩阵表示法中，离散傅立叶变换表示如下：

$\left[\begin{array}{c} X_0 \\ X_1 \\ X_2 \\ \vdots \\ X_{n-1} \end{array}\right]=\left[\begin{array}{ccccc} 1 & 1 & 1 & \cdots & 1 \\ 1 & \alpha & \alpha^2 & \cdots & \alpha^{n-1} \\ 1 & \alpha^2 & \alpha^4 & \cdots & \alpha^{2(n-1)} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \alpha^{n-1} & \alpha^{2(n-1)} & \cdots & \alpha^{(n-1)(n-1)} \end{array}\right]\left[\begin{array}{c} x_0 \\ x_1 \\ x_2 \\ \vdots \\ x_{n-1} \end{array}\right]$

其中 $\alpha = \mathrm{e}^{-\frac{2k \pi\mathrm{i}}{n}}$ 。

快速傅里叶变换 FFT

快速傅立叶变换 (Fast Fourier Transform, FFT) 是一种高效实现 DFT 的算法。

它对傅里叶变换的理论并没有新的发现，

但是对于在计算机系统或者说数字系统中应用离散傅立叶变换，可以说是进了一大步。

快速数论变换（NTT）是快速傅里叶变换（FFT）在数论基础上的实现，本文不会讲述。

在 1965 年，Cooley 和 Tukey 发表了快速傅里叶变换算法。事实上 FFT 早在这之前就被发现过了，但是在当时现代计算机并未问世，人们没有意识到 FFT 的重要性。一些调查者认为 FFT 是由 Runge 和 König 在 1924 年发现的。但事实上高斯早在 1805 年就发明了这个算法，但一直没有发表。

分治法实现（递归）

FFT 算法的基本思想是分治。就 DFT 来说，它分治地来求当 $x = \omega_n^{k}$ 的时候 $f(x)$ 的值。

这里 $\omega_n^{k}=\mathrm{e}^{\frac{2 k\pi\mathrm{i}}{n}}=\cos \frac{2 k\pi}{n}+\mathrm{i} \sin \frac{2 k\pi}{n}$ 即为 $n$ 次单位根的第 $k$ 个。

基-2 FFT 的分治思想体现在将多项式分为奇次项和偶次项处理。

举个例子，对于 $8$ 次整系数多项式：（FFT 强制钦定 $n=2^m$ 为多项式的次数，为了方便处理）

$F(x) = a_0+a_1 x+a_2 x^2+a_3 x^3+a_4 x^4+a_5 x^5+a_6 x^6+a_7 x^7$

把它按奇偶分成两组，并提一个 $x$ ，则

$\begin{aligned} F(x) & =\left(a_0+a_2 x^2+a_4 x^4+a_6 x^6\right)+\left(a_1 x+a_3 x^3+a_5 x^5+a_7 x^7\right) \\ & =\left(a_0+a_2 x^2+a_4 x^4+a_6 x^6\right)+x\left(a_1+a_3 x^2+a_5 x^4+a_7 x^6\right) \end{aligned}$

不妨记

$\begin{aligned} & F_e(x)=a_0+a_2 x+a_4 x^2+a_6 x^3 \\ & F_o(x)=a_1+a_3 x+a_5 x^2+a_7 x^3 \end{aligned}$

则

$F(x) = F_e(x^2) + x F_o{x^2}$

利用偶数次单位根的性质 $\omega_n^i=-\omega_n^{i+n / 2}$ ，以及 $F_e(x^2),F_o(x^2)$ 为偶函数，

可知复平面上 $\omega_n^{i}$ 和 $\omega_n^{i+n/2}$ 的 $F_e(x^2)$ 和 $F_o(x^2)$ 对应的值相等，得到

$\begin{aligned} F\left(\omega_n^k\right) & =F_e((\omega_n^k)^2)+\omega_n^k \times F_o((\omega_n^k)^2) \\[6pt]& =F_e(\omega_n^{2 k})+\omega_n^k \times F_o(\omega_n^{2 k}) \\[6pt]& =F_e(\omega_{n / 2}^k)+\omega_n^k \times F_o(\omega_{n / 2}^k) \end{aligned}$

和

$\begin{aligned} F(\omega_n^{k+n / 2}) &= F_e(\omega_n^{2 k+n})+\omega_n^{k+n / 2} \times F_o(\omega_n^{2 k+n}) \\[6pt]&= F_e(\omega_n^{2 k})-\omega_n^k \times F_o(\omega_n^{2 k}) \\[6pt]&= F_e(\omega_{n / 2}^k)-\omega_n^k \times F_o(\omega_{n / 2}^k) \end{aligned}$

因此我们求出了 $F_e(\omega_{n/2}^{k})$ 和 $F_o(\omega_{n/2}^{k})$ 后，就可以同时求出 $F\left(\omega_n^k\right)$ 和 $F(\omega_n^{k+n / 2})$ 。

于是对 $F_e$ 和 $F_o$ 分别递归即可。时间复杂度 $\mathcal{O}(n \log n)$ 。

但是不放代码了，因为这个实现方式是过不了 $10^6$ 的，常数太大了。

倍增法实现（迭代）

所谓迭代，指的是将原来递归的“拆分”改为了合并答案。

依然以 $8$ 次多项式为例，模拟拆分的过程：

$\left\{x_0, x_1, x_2, x_3, x_4, x_5, x_6, x_7\right\} \\\left\{x_0, x_2, x_4, x_6\right\},\left\{x_1, x_3, x_5, x_7\right\} \\\left\{x_0, x_4\right\}\left\{x_2, x_6\right\},\left\{x_1, x_5\right\},\left\{x_3, x_7\right\} \\\left\{x_0\right\},\left\{x_4\right\},\left\{x_2\right\},\left\{x_6\right\},\left\{x_1\right\},\left\{x_5\right\},\left\{x_3\right\},\left\{x_7\right\}$

结论：可以发现，每个位置就是下标的二进制左右翻转一下，比如 $x_1$ 是 $\mathtt{001}$ ，翻转一下是 $\mathtt{100}$ ，即 $4$ 。

这个变换被称作位逆序置换，证明略。位逆序置换可以做到 $\mathcal{O}(n)$ 从小到大递推实现。

记 $L = 2^k$ ，其中 $k$ 表示二进制数的长度，并设 $R(x)$ 表示翻转后的数（高位补 $0$ ）。显然 $R(0) = (0)_2$ 。

考虑递推。在求 $R(x)$ 时，显然 $R\left(\left\lfloor\frac{x}{2}\right\rfloor\right)$ 是已知的。

因此我们把 $x$ 右移一位后，翻转，再右移一位，就得到了 $x$ 除了二进制下个位以外其他位的翻转结果

考虑个位的翻转结果。如果个位是 $0$ ，翻转后最高位还是 $0$ ，反之需要加上 $\frac{L}{2} = 2^{k-1}$ ，则

$R(x)=\left\lfloor\frac{R\left(\left\lfloor\frac{x}{2}\right\rfloor\right)}{2}\right\rfloor+(x \bmod 2) \times \frac{L}{2}$

这意味着我们把原来 $\mathcal{O}(n\log n)$ 的拆分转变为了 $O(n)$ 的翻转。

实际上总时间复杂度还是 $\mathcal{O}(n\log n)$ ，但是我们却省去了大量常数。

另外还有个蝶形运算优化，和位逆序置换类似，不过我没写这种优化的代码，想看可以看参考文献[1]。

等会儿，这还没结束呢，我们通过 FFT 计算出了原序列的点值形式，我们还要把它转回系数呢

快速傅里叶逆变换 IFFT

傅里叶逆变换可以用傅里叶变换表示。对此我们有两种理解方式。

1. 线性代数角度

IDFT（傅里叶反变换）的作用，是把目标多项式的点值形式转换成系数形式。

而 DFT 本身是个线性变换，可以理解为将目标多项式当作向量，左乘一个矩阵得到变换后的向量，以模拟把单位复根代入多项式的过程：

$\left[\begin{array}{c} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{array}\right]=\left[\begin{array}{cccccc} 1 & 1 & 1 & 1 & \cdots & 1 \\ 1 & \omega_n^1 & \omega_n^2 & \omega_n^3 & \cdots & \omega_n^{n-1} \\ 1 & \omega_n^2 & \omega_n^4 & \omega_n^6 & \cdots & \omega_n^{2(n-1)} \\ 1 & \omega_n^3 & \omega_n^6 & \omega_n^9 & \cdots & \omega_n^{3(n-1)} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \omega_n^{3(n-1)} & \cdots & \omega_n^{(n-1)^2} \end{array}\right]\left[\begin{array}{c} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{array}\right]$

而我们已经得到左侧的结果了，现在要求最右边的。显然这是个矩阵求逆的问题。

由于这个矩阵的元素非常特殊，它的逆矩阵也有特殊的性质：

该矩阵每一项取倒数，再除以变换的长度 $n$ ，就可以得到它的逆矩阵。这里不作证明了。那么

$\frac{1}{\omega_n}=\omega_n^{-1}=\mathrm{e}^{-\frac{2 \pi \mathrm{i}}{n}}=\cos \left(\frac{2 \pi}{n}\right)+\mathrm{i} \sin \left(-\frac{2 \pi}{n}\right)$

好吧，我知道你们可能也不是很会线性代数，所以我们还是讲第二种理解方式吧。

2. 单位根的周期性

利用单位根的周期性同样可以理解 IDFT 与 DFT 之间的转换关系。

考虑构造法。我们已知 $y_i=f\left(\omega_n^i\right),~i = 0,1,\cdots,n-1$ ，求 $\{a_0,a_1,\cdots,a_{n-1}\}$ ，构造多项式如下：

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

相当于把 $\left\{y_0, y_1, y_2, \cdots, y_{n-1}\right\}$ 当作 $A$ 的系数表示法。

设 $b_i = \omega_n^{-i}$ （注意不是 $\mathrm{i}$ ），则 $A$ 在 $x = b_0,b_1,\cdots,b_{n-1}$ 处的点值表示法为 $\{A(b_0),A(b_1),\cdots,A(b_{n-1})\}$ 。

对 $A(x)$ 的定义式做一下变换，可以将 $A(b_k)$ 表示为

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

记 $S\left(\omega_n^a\right)=\sum_{i=0}^{n-1}\left(\omega_n^a\right)^i$

当 $a=0\pmod{n}$ 时，$S(\omega_n^{a}) = n$ 。

当 $a \ne 0 \pmod{n}$ 时，考虑错位相减

$\begin{aligned} S\left(\omega_n^a\right) & =\sum_{i=0}^{n-1}\left(\omega_n^a\right)^i \\ \omega_n^a S\left(\omega_n^a\right) & =\sum_{i=1}^n\left(\omega_n^a\right)^i \\ S\left(\omega_n^a\right) & =\frac{\left(\omega_n^a\right)^n-\left(\omega_n^a\right)^0}{\omega_n^a-1}=0 \end{aligned}$

即

$S\left(\omega_n^a\right)= \begin{cases}n, & a=0 \\ 0, & a \neq 0\end{cases}$

代回原式可得

$A\left(b_k\right)=\sum_{j=0}^{n-1} a_j S\left(\omega_n^{j-k}\right)=a_k \cdot n$

也就是说，给定点 $b_i = \omega_n^{-i}$ ，则 $A$ 的点值表示法为

$\begin{aligned} & \left\{\left(b_0, A\left(b_0\right)\right),\left(b_1, A\left(b_1\right)\right), \cdots,\left(b_{n-1}, A\left(b_{n-1}\right)\right)\right\} \\ = & \left\{\left(b_0, a_0 \cdot n\right),\left(b_1, a_1 \cdot n\right), \cdots,\left(b_{n-1}, a_{n-1} \cdot n\right)\right\} \end{aligned}$

因此我们取其单位根为其倒数，对 $\left\{y_0, y_1, y_2, \cdots, y_{n-1}\right\}$ 跑一遍 FFT ，然后除以 $n$ 即可得到 $f(x)$ 的系数表示。

代码实现

至此，我们终于讲完了 FFT 算法的全部，放张思维导图

时间复杂度 $\mathcal{O}(n \log n)$

代码：

#include <bits/stdc++.h>
using namespace std;
// #define int long long
// #define INF 0x3f3f3f3f3f3f3f3f
typedef complex<double> Complex;
void up(int &x,int y) { x < y ? x = y : 0; }
void down(int &x,int y) { x > y ? x = y : 0; }
#define N ((int)(1e7 + 15))

const double pi = acos(-1.0);
Complex a[N], b[N];
int l, r[N], limit = 1;
void FFT(Complex *A, int type)
{
    for(int i = 0; i < limit; i++) if(i < r[i]) swap(A[i], A[r[i]]);
    for(int mid = 1; mid < limit; mid *= 2)
    {
        Complex Wn(cos(pi / mid), type * sin(pi / mid));
        for(int j = 0; j < limit; j += (mid * 2))
        {
            Complex w(1, 0);
            for(int k = 0; k < mid; k++, w = w * Wn)
            {
                Complex x = A[j + k], y = w * A[j + k + mid];
                A[j + k] = x + y; A[j + k + mid] = x - y;
            }
        }
    }
}
signed main()
{
    ios::sync_with_stdio(0);
    cin.tie(0);cout.tie(0);
    // freopen("check.in","r",stdin);
    // freopen("check.out","w",stdout);
    int n, m; cin >> n >> m;
    for(int i = 0, x; i <= n; i++) cin >> x, a[i] = x;
    for(int i = 0, x; i <= m; i++) cin >> x, b[i] = x;
    for(; limit <= n + m; l++, limit *= 2);
    for(int i = 0; i < limit; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1));
    FFT(a, 1); FFT(b, 1);
    for(int i = 0; i < limit; i++) a[i] = a[i] * b[i];
    FFT(a, -1);
    for(int i = 0; i <= n + m; i++)
        cout << (int)(a[i].real() / limit + 0.5) << " \n"[i == n + m];
    return 0;
}