多项式与 FFT 备忘录

2022-07-24

很多两年前学过的东西已经忘了, 虽然很大一部分原因在于当时就没有学明白, 但还是有一种原有的认知落空的感觉.

综述

在 OI/XCPC 中, FFT 常用于加速卷积的计算. 此外, 配合牛顿迭代得到的一系列多项式算法提供了对生成函数进行快速运算的方法, 因此在许多计数问题中也 FFT 也很有用.

定义与前置知识

多项式和形式幂级数

对于环 $R$, $a_1, a_2, \ldots, a_n \in R$ 且 $a_n \neq 0$, 称

为环 $R$ 上次数为 $n$ 的多项式 $f(x)$. 其中次数也称作度, 记作 $\deg f$.

如果允许无限项存在, 则可以得到 $R$ 上的形式幂级数

环 $R$ 上的多项式构成了一个环, 记为 $R[x]$. 类似地, 环 $R$ 上的形式幂级数也构成了一个环, 记作 $R[[x]]$.

先来讨论在复数域 $\mathbb{C}$ 上的情况. 当然也可以是其他环.

多项式的点值表示

对于度数为 $n$ 的多项式 $f(x) = \sum_{k=0}^n a_kx^k$, 带入 $n+1$ 个点 $x_0, x_1, \ldots, x_n$ 可得到 $f(x)$ 对应的点值表示 $f(x_0), f(x_1), \ldots, f(x_n)$. 表示成矩阵形式为

如果这 $n+1$ 个点互不相同, 那么中间的矩阵可逆, 根据 $n+1$ 个点对 $(x_k, f(x_k))$ 就可以唯一确定 $f(x)$.

卷积和循环卷积

对于一个长为 $n$ 的序列 $a$ 和一个长为 $m$ 的序列 $b$, 有 $a$ 和 $b$ 的卷积为

卷积可以用来表示多项式相乘的结果. 如果把 $a$ 和 $b$ 分别看作多项式 $f(x)$ 和 $g(x)$ 的各项系数, 那么两多项式之积 $f(x)g(x)$ 可表示为

对于两个长度为 $n$ 的序列 $a$, $b$, 有 $a$ 和 $b$ 的循环卷积为

离散傅里叶变换及其逆变换

离散傅里叶变换 (Discrete Fourier Transform, DFT) 将一个长为 $N$ 的序列 $x_{0}, x_{1}, \ldots, x_{N-1}$ 变换为一个长为 $N$ 的复数序列 $X_{0}, X_{1}, \ldots, X_{N-1}$. 具体而言

同时有逆变换 (Inverse Discrete Fourier Transform, IDFT)

逆变换的正确性可以通过代入验证.

记 $\omega_n = e^{i\frac{2\pi}{n}}$, 根据复数的性质 $1, \omega_n, \ldots, \omega_n^{n-1}$ 组成了方程 $x^n = 1$ 的 $n$ 个不同的根, 称这 $n$ 个根为单位根.

如果将 $x_n$ 看作多项式 $f(x)$ 的系数, 那么 $X_k$ 就是 $f(x)$ 在 $x = e^{-i\frac{2\pi}{N}n}$ 处的点值. 因此 DFT 可以理解为多项式在 $n$ 个单位根处的点值表示. 同时, 实现上在计算 DFT 时往往使用 $1, \omega_n, \ldots, \omega_n^{n-1}$, 而非它们的倒数.

卷积定理

对于两个长度为 $n$ 的序列 $a$ 和 $b$, 记 $a$ 和 $b$ 的循环卷积为 $c$, 那么有

也就是说, 对 $a$ 和 $b$ 分别做 DFT, 两个结果逐项相乘就得到 $\operatorname{DFT}(c)$, 再做一次 IDFT 就可以得到 $c$.

证明可以从 DFT 与多项式在 $n$ 个单位根处的点值表示的关系入手. 设以 $a$, $b$ 和 $c$ 作为各项系数的多项式分别为 $f(x)$, $g(x)$ 和 $h(x)$. 对于 $s\in \mathbb{Z}$, 有

也就是说 $\operatorname{DFT}(a) \cdot \operatorname{DFT}(b) = \operatorname{DFT}(c)$. 等式两端同时做 IDFT 就得到卷积定理.

利用卷积定理得到的计算结果实际上是循环卷积, 而非卷积. 但实际上可以通过补 0 使 $a$ 和 $b$ 的项数变为 $2n$, 得到项数为 $2n$ 的 $c$. 此时的循环卷积相当于普通的卷积. 在 $a$ 和 $b$ 长度不同且需要计算卷积时同样可以这样做.

快速傅里叶变换

快速傅里叶变化 (Fast Fourier Transform, FFT) 能够在 $O(n\log n)$ 的时间复杂度内实现 DFT/IDFT 的计算, 再利用卷积定理即可在 $O(n\log n)$ 的时间复杂度内计算卷积.

对于 $n\in\mathbb{N}_+$, $k\in\mathbb{Z}$, 单位根 $\omega_n$ 有性质 $\omega_n^k=\omega_{2n}^{2k}$, $\omega_{2n}^{k+n}=-\omega_{2n}^k$. 接下来的推导将用到这些性质.

出于简便考虑, 只讨论对长度为 2 的幂次的序列做 FFT 的情况. 其余情况可以通过补 0 让长度变为 2 的幂次.

先来讨论 DFT. 考虑分治, 将当前序列分作奇数项和偶数项, 假定已经得到了奇数项对应 DFT 的各项值 $g(\omega_n^k)$, 以及偶数项对应的各项值 $h(\omega_n^k)$. 现在需要合并两者得到 $f(\omega_{2n}^k)$. 具体有

对于 IDFT, 一个直观的做法是将代入的单位根替换为原来的倒数, 再将所有值除以序列的长度 $n$. 但还有另一种做法. 考虑单位根的周期性, 将后 $n-1$ 个元素翻转就等同于代入原有单位根的倒数. 翻转后再将所有值除以 $n$ 就好了.

此时就得到了一个可以用递归实现且时间复杂度为 $O(n\log n)$ 的算法. 但递归实现常数较大, 接下来利用一些技巧将其更改为非递归实现.

蝴蝶变换 / 位逆序置换

每一个位置 $i$, 在递归的过程中位置 $i$ 最终到达的位置恰好是 $i$ 二进制下各位翻转后得到的数. 这个翻转后的值可以在 $O(n)$ 的时间内递推出来.

1
2
3
rev[0] = 0;
for (int i = 1; i < n; ++i) // n = (1 << l)
  rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));

借助于 rev[i] 就可以实现非递归的 FFT.

代码实现

> 点击显示 / 隐藏 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
// UOJ #34
// DeP
#include <cstdio>
#include <complex>
#include <algorithm>
using namespace std;

using Complex = complex<double>;
constexpr int LOG = 18, MAXN = 1 << LOG | 1;
constexpr double PI = acos(-1.0);

namespace Poly {
  int rev[MAXN];
  void init(int lim, int l) {
    for (int i = 1; i < lim; ++i)
      rev[i] = (rev[i >> 1] >> 1) | ((i & 1) << (l - 1));
  }

  void FFT(Complex* f, int lim, int type) {
    for (int i = 1; i < lim; ++i)
      if (i < rev[i]) swap(f[i], f[rev[i]]);
    for (int mid = 1; mid < lim; mid <<= 1) {
      Complex wn = exp(Complex(0.0, PI / mid));
      for (int i = 0; i < lim; i += mid << 1) {
        Complex w = 1.0;
        for (int j = 0; j < mid; ++j, w *= wn) {
          Complex f0 = f[i+j], f1 = w * f[i+j+mid];
          f[i+j] = f0 + f1, f[i+j+mid] = f0 - f1;
        }
      }
    }
    if (type < 0) {
      for (int i = 0; i < lim; ++i) f[i] /= lim;
      reverse(f + 1, f + lim);
    }
  }
}

int n, m;
Complex f[MAXN], g[MAXN];

int main() {
  scanf("%d%d", &n, &m), ++n, ++m;
  for (int x, i = 0; i < n; ++i)
    scanf("%d", &x), f[i] = Complex(x);
  for (int x, i = 0; i < m; ++i)
    scanf("%d", &x), g[i] = Complex(x);

  int lim = 1, l = 0;
  while (lim < n + m - 1) lim <<= 1, ++l;
  Poly::init(lim, l), Poly::FFT(f, lim, 1), Poly::FFT(g, lim, 1);
  for (int i = 0; i < lim; ++i) f[i] = f[i] * g[i];
  Poly::FFT(f, lim, -1);

  for (int i = 0; i < n + m - 1; ++i)
    printf("%.0f%c", f[i].real(), " \n"[i == n + m - 2]);
  return 0;
}

预处理单位根

在代码实现中, 单位根的幂次被反复计算, 可以将这部分预处理出来. 在多次进行 DFT/IDFT 时可以带来明显的常数优化.

快速数论变换

数论变换 (Number Theoretic Transform, NTT) 是 DFT 在模质数 $p$ 意义下的域 $\mathbb{F}_p$ 上的情况, 对模数 $p$ 有一定要求. 不过很多人也把此种情况下的 FFT 直接称作 NTT.

常用的 NTT 模数满足 $p=a\cdot 2^{b} + 1$, 要求序列长度不超过 $2^b$. 记 $p$ 的原根为 $g$, 那么 $g^{\frac{p-1}{n}}$ 就相当于 $\omega_n$, 且具有相似的性质.

是原根为 $3$ 的 NTT 模数. 写三模数 NTT 的时候可能都会用到.

应用

加速卷积的计算

后文的其他应用本质上也是将问题转化为卷积的计算, 从而利用 FFT 优化时间复杂度. 普通的卷积当然可以直接算, 下面是两个不那么显然的卷积.

二项式反演中常出现的式子, 也有人把这个叫做减法卷积. 记 $j = n + k - i$, 则有

反转 $a$ 得到 $a’$, $a’$ 与 $b$ 卷积的第 $n+k$ 项即为 $c_k$.

同 (1) 类似, 记 $j = n - k - i$, 则有

反转 $a$ 得到 $a’$, $a’$ 与 $b$ 卷积的第 $n-k$ 项即为 $c_k$.

字符串匹配

对于长为 $n$ 字符串 $s$, $s$ 在字符串 $t$ 的位置 $k$ 处完全匹配当且仅当

展开后是卷积的形式, 利用 FFT 计算就好了.

分治与 FFT

在 $g_{1\ldots n}$ 已知的情况下, 对于形同

的式子, 可以利用分治的技巧在 $O(n\log^2 n)$ 的时间复杂度计算 $f_{1\ldots n}$. 具体而言, 对于一段区间, 先计算左半区间的值, 然后将左半区间上值对右半区间的贡献累计到右半区间上, 再计算右半区间的值. 其实就是 CDQ 分治.

1
2
3
4
5
6
7
8
9
void solve(int l, int r) { // [l, r]
  if (l == r) return;
  int mid = (l + r) / 2;
  solve(l, mid);
  Poly::Mul(f + l, mid - l + 1, g, r - l + 1, h);
  for (int i = mid + 1; i <= r; ++i)
    f[i] = (f[i] + h[i - l]) % P;
  solve(mid + 1, r);
}

另外, 这个式子比较特殊, 经过一定的转化可以得到利用多项式求逆计算的方法, 时间复杂度 $O(n\log n)$.

多项式的运算

记 $[x^k]f(x)$ 表示多项式 $f(x)$ 中 $x^k$ 的系数.

模多项式

首先介绍多项式除法. 对于多项式 $f(x)$ 和 $g(x)$, 如果存在多项式 $Q(x)$ 和 $R(x)$ 满足

且 $\deg R < \deg g$. 那么这样的 $Q(x)$ 和 $R(x)$ 唯一, 并称 $Q(x)$ 是 $f(x)$ 除以 $g(x)$ 的商, $R(x)$ 为余数.

同时, 称 $f(x)$ 和 $R(x)$ 在模 $g(x)$ 意义下同余, 记作

一个显然的性质是, 对于 $g(x)$ 的根 $x_0$, 考虑到 $g(x_0) = 0$, 则一定会有

特别地, 当 $g(x) = x^n - 1$ 时, 上式中的 $x_0$ 可以取单位根 $\omega_n$. 所以 DFT 和 IDFT 的可以理解为是在模 $x^n - 1$ 意义下的变换.

在后文的运算中往往只关注模 $x^n$ 时的结果, 形式上就是多项式系数的前 $n$ 项.

多项式牛顿迭代

已知多项式 $g(x)$, 多项式牛顿迭代提供了计算满足以下条件 $f(x)$ 的方法.

考虑倍增. 当 $n=1$ 时, 单独计算满足 $[x^0]g(f(x)) = 0$ 的 $[x^0]f(x)$.

当 $n>1$ 时. 记 $f_0(x)$ 是在模 $x^n$ 意义下的结果. 假定 $f_0(x)$ 已知, 现在需要计算模 $x^{2n}$ 下的结果 $f(x)$. 对 $g(f(x))$ 在 $f_0(x)$ 处进行泰勒展开, 有

其中 $g^{(k)}(f_0(x))$ 指 $g(f(x))$ 对 $f(x)$ 求 $k$ 阶导.

又因为 $f(x) - f_0(x)$ 最低的非零项次数大于等于 $x^n$, 所以当 $k\ge2$ 时有

皆大欢喜. 保留 $k=0$ 和 $k=1$ 两项就好了. 整理后得到

其中导数依旧是对 $f(x)$ 求导. 通过构造不同的 $g(f(x))$, 就可以实现不同的运算.

多项式求逆

已知多项式 $h(x)$, 求模 $x^m$ 下多项式的 $f(x)$, 满足 $h(x)f(x) \equiv 1 \pmod{x^m}$. 有时也把 $f(x)$ 记作 $h^{-1}(x)$.

首先 $[x^0]f(x) = \left([x^0]h(x)\right)^{-1}$. 令 $g(x) = f(x) - h^{-1}(x)$, 则有

倍增计算就好了, 时间复杂度 $O(n\log n)$.

多项式开方

已知多项式 $h(x)$, 求模 $x^m$ 下多项式的 $f(x)$, 满足 $f^2(x) \equiv h(x) \pmod{x^m}$. 有时也把 $f(x)$ 记作 $\sqrt{h(x)}$.

如果是在 $\mathbb{F}_p$ 中考虑, $[x^0]f(x)$ 可以通过计算模 $p$ 下的二次剩余得到. 令 $g(x) = f^2(x) - h(x)$, 则有

倍增计算就好了, 时间复杂度 $O(n\log n)$.

多项式对数函数

对于多项式 $A(x) = \sum_{k\ge 1} a_kx^k$, 定义多项式对数函数

如果多项式 $h(x)$ 满足 $[x^0]h(x) = 1$, 那么有

如果在积分时线性预处理逆元, 那么时间复杂度为 $O(n)$.

多项式指数函数

对于多项式 $A(x) = \sum_{k\ge 1} a_kx^k$, 定义多项式指数函数

已知多项式 $h(x)$, 满足 $[x^0]h(x) = 0$, 求模 $x^m$ 下的 $\exp h(x)$. 也就是求 $f(x)$, 满足 $\ln f(x) \equiv h(x) \pmod{x^m}$.

首先 $[x^0]f(x) = 1$. 令 $g(x) = \ln f(x) - h(x)$, 则有

倍增计算就好了, 时间复杂度 $O(n\log n)$.

多项式幂函数

当 $[x^0]f(x) = 1$ 时, 有

如果 $f(x)$ 在 $\mathbb{F}_p$ 上, 有 $f^p(x) = 1$. 在 $k$ 非常大的时候可能会用到.

当 $[x^0]f(x) \neq 0$ 时, 找到 $f(x)$ 次数最低的不为零项 $f_ix^i$, 将 $f_ix^i$ 提取出来就好了.

时间复杂度 $O(n\log n)$.

参考资料

© 2019-2023 DepletedPrism's Blog - Powered by Hexo - Theme Chan