本节点所需要的前置研发项目:多项式基础

在上一篇文章中,我们知道多项式有两种表示方法,这两种表示方法分别对应了求多项式卷积的两种方式,然而效率都比较低。现在我们考虑对第二种方法——点值相乘法进行优化。

如果我们可以找到一组特殊的点值,它满足一些奇妙的性质使得问题可以被分治求解,那么问题就解决了。

世界上恰好就存在这么一种东西,它叫单位根。首先我们需要介绍一些概念。

前置知识

复数

考虑下面这个方程的解。
$$
x^2+1=0
$$
显而易见它没有实数解。我们强行令它有解,称这个解为 $\mathrm{i}$,即虚数单位。(注意这里使用的是罗马体,而一般的循环变量 $i$ 是斜体)

有一种数,形如 $a+b\mathrm{i}$,其中 $a,b\in\mathbb{R}$,我们称这种数为复数,a 是实部,b 是虚部。一个复数在复平面上对应一个向量,它与 x 轴正方向的夹角称为幅角,长度称为。以下是复数的一些运算法则。
$$
(a_1+b_1\mathrm{i})\pm(a_2+b_2\mathrm{i})=(a_1\pm a_2)+(b_1\pm b_2)\mathrm{i}
$$

$$
(a_1+b_1\mathrm{i})\times(a_2+b_2\mathrm{i})=(a_1a_2-b_1b_2)+(a_1b_2+a_2b_1)\mathrm{i}
$$

这两个式子暴力合并同类项即可得到。(记住 $\mathrm{i}^2=-1$)

复数相乘有一个重要的性质,即幅角相加,模相乘

单位根

满足 $x^n-1=0$ 的 x 称为 n 次单位根。若 $\omega$ 为单位根,且 $\omega^0, \omega^1, \omega^2, \cdots, \omega^{n-1}$对应所有的 n 次单位根,那么称 $\omega$ 为本原单位根

由上面复数相乘的性质可以知道,$\omega_n=\cos\frac{2\pi}{n}+\mathrm{i}\sin\frac{2\pi}{n}$ 是一个本原单位根,并且本原单位根的所有次幂 n 等分单位圆。

DFT 和 IDFT

先放上公式。

设 $b$ 是 $a$ 的离散傅里叶变换,记为 $b=\mathcal{F(a)}$

DFT:$b_k=\sum\limits_{i=0}^{n-1}a_i\omega_n^{ki}(0\leq k<n)$

IDFT:$a_k=\frac{1}{n}\sum\limits_{i=0}^{n-1}b_i\omega_n^{-ki}(0\leq k<n)$

不难发现 DFT 就是带入点值,IDFT 就是解出原来的多项式系数表示。

证明

将 DFT 的式子代入 IDFT 的式子可以得到 $a_k=\frac{1}{n}\sum\limits_{i=0}^{n-1}\omega_n^{-ki}\sum\limits_{j=0}^{n-1}a_i\omega_n^{ij}$

改变求和顺序得到 $a_k=\frac{1}{n}\sum\limits_{j=0}^{n-1}a_i\sum\limits_{i=0}^{n-1}\omega_n^{ij}\omega_n^{-ki}=\frac{1}{n}\sum\limits_{j=0}^{n-1}a_i\sum\limits_{i=0}^{n-1}\omega_n^{i(j-k)}$

考虑 $\sum\limits_{i=0}^{n-1}\omega_n^{i(j-k)}$ 的值。

  • 若 $j=k$,原式等于 n。
  • 若 $j\ne k$,原式等于 $\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i$ ,显然这是一个等比数列,根据等比数列求和公式可以得到 $\sum\limits_{i=0}^{n-1}(\omega_n^{j-k})^i=\dfrac{1-(\omega_n^{j-k})^n}{1-\omega_n^{j-k}}=\dfrac{1-(\omega_n^{n})^{j-k}}{1-\omega_n^{j-k}}$。由于 $\omega_n^n=1$,原式等于 0。

也就是说后面的式子只在 $j=k$ 时为 n,其余时候都为 0。等式右边的式子也就等于 $a_k$。

至此我们证明了 IDFT 和 DFT 是互逆的。

循环卷积

对于两个长度为 $n$ 的数列 $a$ 和 $b$,定义一个数列 $c$ 满足
$$
c_k=\sum\limits_{i+j\mod n=k}a_ib_j
$$
则称 $c$ 是 $a$ 和 $b$ 的循环卷积,记作 $c=a*b$

循环卷积满足这个性质:$\mathcal{F}(ab)=\mathcal{F}(a)\mathcal{F(b)}$

因此,要求两个多项式的卷积,我们可以先将他们的长度拓展到原来长度之和(缺位补零),然后分别做一次 DFT,乘起来再 IDFT 就可以了。如果直接按照定义计算,时间复杂度为 $O(n^2)$。

至于为什么要拓展长度,因为 DFT 所求的是循环卷积。需要拓展长度才能保证正确性。

FFT(快速傅里叶变换)

FFT 是快速计算 DFT 和 IDFT 的算法。

关于单位根,我们还有以下两条性质。

  • $(\omega_{2n}^k)^2=\omega_{n}^k$
  • $\omega_{2n}^{n+k}=-\omega_{2n}^k$

这两个性质到单位圆上去看看就可以发现了。

设 $n = 2m$,我们将要做 DFT / IDFT 的式子次数奇偶分类。

$A(x)=\sum\limits_{i=0}^{n-1}a_ix^i=\sum\limits_{i=0}^{m-1}a_{2i}x^{2i}+\sum\limits_{i=0}^{m-1}a_{2i+1}x^{2i+1}=\sum\limits_{i=0}^{m-1}a_{2i}x^{2i}+x\sum\limits_{i=0}^{m-1}a_{2i+1}x^{2i}$

设 $A_0(x)=\sum\limits_{i=0}^{m-1}a_{2i}x^{i}$ ,$A_1(x)=\sum\limits_{i=0}^{m-1}a_{2i+1}x^{i}$ ,那么 $A(x)=A_0(x^2)+xA_1(x^2)$。

下面考虑带入单位根,对于 $0\leq k <m$

  • $A(\omega_n^k)=A_0((\omega_n^k)^2)+\omega_n^kA_0((\omega_n^k)^2)=A_0(\omega_m^k)+\omega_n^kA_0(\omega_m^k)$
  • $A(\omega_n^{m+k})=A_0((\omega_n^{m+k})^2)+\omega_n^{m+k}A_0((\omega_n^k)^2)=A_0(\omega_m^{k})-\omega_n^{k}A_0(\omega_m^{k})$

以上两个式子称为蝴蝶操作

因此,如果我们得到了 $A_0(x)$ 和 $A_1(x)$,我们可以 $O(n)$ 计算出 $A(x)$。这样递归计算就可以了。

由于每次分成两个子任务,所以$T(n)=T(n/2)+O(n)=O(l\log n)$。

注意 FFT 要求每次的 n 都为偶数,所以我们需要找到一个 $lim=2^l \geq deg(a)+deg(b)$

IDFT 操作同理。然而在实际编程中,对比 DFT 和 IDFT 的式子,可以发现只有两个位置不一样:

  1. 多了一个 n 分之一。
  2. 单位根的指数变成了负的。

那么我们在函数中带入一个系数,如果系数为 1 就是 DFT,系数为 -1 就是 IDFT。另外如果带入 -1 还要给结果除掉 n。

Q:代码呢?

A:由于递归很慢,所以说我实际上没有写过一个递归 FFT 的板子,也就没有了QwQ。然而接下来会讲到另一种更好的实现方法。

位逆序置换

来看下 FFT 的运算过程,下标用二进制表示(图片来自网络)

不难发现每一行左边到右边的二进制倒过来了。所以我们预处理出数组 $r[i]$,表示 i 二进制反转后得到的数,这样就得到了最终状态。然后我们第一次将距离为 1 的位置合并,第二次将距离为 2 的位置合并,第三次将距离为 4 的位置合并······这样就可以将递归变成递推了。

二进制反转可以 $O(n)$ 动态规划得到。设 $lim=2^l \geq deg(a)+deg(b)$ 状态转移方程为:$r_i=(r_{i>>1}>>1) \or ((i \and 1)<<(l-1))$,其中 $i\in[0, lim)$。(感性理解)

代码

下面给出 FFT 的 C++ 实现。

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
#include<bits/stdc++.h>
using namespace std;

inline int get() {
int x = 0, f = 1; char c = getchar();
while(!isdigit(c)) { if(c == '-') f = -1; c = getchar(); }
while(isdigit(c)) { x = x * 10 + c - '0'; c = getchar(); }
return x * f;
}

const double Pi = acos(-1.0);
const int N = 5e6 + 5;
struct Complex { // 复数结构体
double x, y;
Complex(double a = 0, double b = 0) { x = a, y = b; }
friend Complex operator +(Complex a, Complex b) { return Complex(a.x + b.x, a.y + b.y); }
friend Complex operator -(Complex a, Complex b) { return Complex(a.x - b.x, a.y - b.y); }
friend Complex operator *(Complex a, Complex b) { return Complex(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }
} a[N], b[N];
int n, m, lim = 1, l;
int r[N];

inline void FFT(Complex *A, int type) { // type=1 表示 DFT,type=2 表示 IDFT
for(register int i = 0; i < lim; i++) if(i < r[i]) swap(A[i],A[r[i]]); // 进行位逆序置换得到最终状态
for(register int mid = 1; mid < lim; mid <<= 1) { // 递推,枚举 m
Complex Wn = Complex(cos(Pi / mid), type * sin(Pi / mid)); // 本原单位根
for(register int r = (mid << 1), j = 0; j < lim; j += r) { // 枚举合并位置
Complex w = Complex(1, 0);
for(register int k = 0; k < mid; k++, w = w * Wn) { // 将 j + k 和 j + mid + k 合并
Complex x = A[j + k], y = w * A[j + mid + k]; // 见上方式子
A[j + k] = x + y;
A[j + mid + k] = x - y;
}
}
}
}

int main() {
n = get(), m = get();
for(register int i = 0; i <= n; i++) a[i].x = get();
for(register int i = 0; i <= m; i++) b[i].x = get();
while(lim <= n + m) lim <<= 1, l++; // 找到 lim
for(register int i = 0; i < lim; i++) r[i] = (r[i >> 1] >> 1) | ((i & 1) << (l - 1)); // DP 出反转数组
FFT(a, 1);
FFT(b, 1); // DFT
for(register int i = 0; i < lim; i++) a[i] = a[i] * b[i];
FFT(a, -1); // IDFT
for(register int i = 0; i <= n + m; i++) printf("%d ", (int)(a[i].x / lim + 0.5)); // 结果除掉 n
return 0;
}