快速 Fourier 变换（Fast Fourier Transform, FFT），是一种快速地将多项式从系数表示法转化为点值表示法的算法。借助 FFT 可以高效地实现卷积运算。

离散 Fourier 变换

可参考：（链接未公开）

离散 Fourier 变换（Discrete Fourier Transform, DFT）是对有限长数组的一种可逆变换。

定义

设 (x_0, \dots, x_{N - 1}) 为一个有限长数组，定义它的 DFT 为序列 (X_0, \dots, X_{N - 1})，满足：

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

\omega_N = \mathrm{e}^{\mathrm{i} \frac{2 \pi}{N}} 即为 N 次单位根，对应方程 z^N = 1 的解，因此 DFT 也可以写成：

X_k = \sum_{n = 0}^{N - 1} x_n \omega_N^{-kn}.

其逆变换称为逆离散 Fourier 变换（IDFT），公式为：

x_n = \frac{1}{N} \sum_{k = 0}^{N - 1} X_k \omega_N^{kn}.

对于一个 N - 1 次多项式 f(z) = x_0 + x_1 z + \cdots + x_{N - 1} z^{N - 1}，对其系数做 DFT 得到的值 X_k，对应多项式在单位根上的取值：

X_k = f \left(\omega_N^{-k} \right).

提示

DFT 也可以定义为

X_k = \sum_{n = 0}^{N - 1} x_n \mathrm{e}^{\mathrm{i} \frac{2 \pi}{N} kn} = \sum_{n = 0}^{N - 1} x_n \omega_N^{kn}.

定义的区别在于，这里单位根的幂次上没有负号。此时，相应地 IDFT 应该带上负号：

x_n = \frac{1}{N} \sum_{k = 0}^{N - 1} X_k \omega_N^{-kn}.

不考虑物理意义，只进行卷积运算时，两个定义没有很大的区别。

矩阵形式

DFT 可视为对向量 (x_0, \dots, x_{N - 1})^{\mathrm T} 的线性变换，因此可以写成矩阵形式：

\begin{bmatrix} X_{0} \\ X_{1} \\ X_{2} \\ \vdots \\ X_{N-1} \end{bmatrix} = \begin{bmatrix} 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{bmatrix} \begin{bmatrix} x_{0} \\ x_{1} \\ x_{2} \\ \vdots \\ x_{N-1} \end{bmatrix}

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

快速 Fourier 变换

FFT 是 DFT 的高效算法。对长度为 2 的幂（N = 2^l）的数组，FFT 可在 O(N \log N) 的时间复杂度内实现 DFT。

分治

对于一个长度为 2 的幂的数组 (a_0, a_1, \dots, a_{N-1})，计算其 FFT，就是计算单位根 \omega_N^k 上多项式 f(x) = a_0 + a_1 x + \cdots + a_{N - 1} x^{N - 1} 的取值。

以长度为 8 的多项式为例，考虑将 f(x) 的系数分奇偶组合：

\begin{aligned} f(x) &= (a_0+a_2x^2+a_4x^4+a_6x^6) + (a_1x+a_3x^3+a_5x^5+a_7x^7)\\ &= (a_0+a_2x^2+a_4x^4+a_6x^6) + x(a_1+a_3x^2+a_5x^4+a_7x^6) \end{aligned}

当数组的长度不是 2 的幂时，需要对数组补 0。

实现

using ld = long double;
using comp = complex<ld>;
using ll = long long;
constexpr size_t MAXN = 4e6 + 5;
constexpr ld pi = 3.14159265358979297537819636906696188820;

comp tmp[MAXN], f[MAXN], g[MAXN], h[MAXN];

void fft(comp *f, size_t n, int rev)
{
    if (n == 1)
    {
        return;
    }
    for (size_t i = 0; i < n; i++)
    {
        tmp[i] = f[i];
    }
    for (size_t i = 0; i < n / 2; i++)
    {
        f[i] = tmp[i * 2]; // g
        f[i + n / 2] = tmp[i * 2 + 1]; // h
    }
    fft(f, n / 2, rev);
    fft(f + n / 2, n / 2, rev);
    comp o(1, 0), d(cos(2 * pi / n), rev * sin(2 * pi / n));
    for (size_t i = 0; i < n / 2; i++)
    {
        tmp[i] = f[i] + o * f[i + n / 2];
        tmp[i + n / 2] = f[i] - o * f[i + n / 2];
        o *= d;
    }
    for (size_t i = 0; i < n; i++)
    {
        f[i] = tmp[i];
    }
}

应用

卷积

借助离散 Fourier 变换，可以计算环形数组的循环卷积 / 周期数列的卷积：

(x * y)_n = \sum_{k = 0}^{N - 1} x_k y _{n - k}.

该卷积的特点是只在一个周期内进行（无限个周期相加将得到无穷大），卷积后得到的依然是一个周期为 N 的数列。

此时，周期卷积的 DFT 等于 DFT 的乘积，即：

\mathscr{F}[x * y] = \mathscr{F}[x] \mathscr{F}[y].

但是，由于 FFT 只能处理长度为 2 的幂的情况，直接应用此方法计算循环卷积，也要求长度是 2 的幂。

多项式乘法

多项式乘法的过程，是对系数进行卷积。因此，可以使用 FFT 加速多项式的乘法。

算法

两个次数分别为 m, n 的多项式 f(x), g(x) 做卷积时，所得多项式的次数为 m + n，长度为 m + n + 1。为了从点值表示法还原出这个多项式的系数，要求至少有 m + n + 1 个点值。因此，FFT 的长度应选择 L = 2^{\lceil \log_2(m + n + 1) \rceil}。

初始时应对多项式 f(x), g(x) 进行高次补 0 至长度为 L，再对补 0 后的两个多项式的系数做 DFT，得到两个频域的序列 F_k, G_k。令 H_k = F_k G_k，再对序列 H_k 做一次 IDFT，就得到了两个多项式的乘积 h(x)。

代码

using comp = complex<double>;
comp f[MAXN], g[MAXN], h[MAXN];

void poly_mul(size_t m, size_t n)
{
    size_t len = 1;
    while (len < max(n + 1, m + 1))
    {
        len *= 2;
    }
    
}

模 x^N - 1 多项式环

当补 0 的长度不够时，两个多项式相乘后的系数在 FFT 会“存不下”，然而此时，FFT 和 iFFT 仍然可以算出一个系数，这个系数是什么？

联想前面提到的循环卷积，不难猜出这里同样是产生了某种循环。当序列的长度只有 N 时，x^N 和 1 将发生混叠，变得无法区分。

从数学上看，这相当于取了 x^N = 1 的一个等价类，也相当于多项式对 x^N - 1 取模。

例题

2021 澳门站 E. Pass the Ball!