发布日期: 2024-01-27

更新日期: 2024-02-04

文章字数: 3.9k

快速傅里叶变换 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;
}