「最喜欢的一集」- 卷积初步与 FFT

卷积#

引入#

让我们从一个非常简单的问题入手：假设有两个质量均匀的骰子 A 和 B，都分别写了 $1\sim 6$ 这 $6$ 个数字，现在同时扔出他们，投到的点数之间相互独立，则投到的点数和及对应的概率会是多少？

显然，点数一定在 $2\sim 12$ 之间，对应的概率不难计算。记骰子 A 投出 $i$ 的概率为 $A_i$，投出的数为 $A$；骰子 B 投出 $i$ 的概率为 $B_i$，投出的数为 $B$。那么，以点数和为 $8$ 为例，则 $P(A+B = 8) = A_2\times B_6 + A_3\times B_5+A_4\times B_4+A_5\times B_3+A_6\times B_6 = \dfrac{5}{36}$。同理，我们可以其他的点数和对应的式子及结果也写出来，可以统一写成：$P(A+B = X) = \sum\limits_{i+j = X} A_i \times B_j$。不难发现，如果骰子上的数很多，或质量不均匀，这个式子是恒成立的。

定义#

同理，假设 $a$ 和 $b$ 是两个数组，我们记 $a*b$ 为 $a$ 和 $b$ 的卷积，其中，$(a*b)_n = \sum\limits_{i+j = n} a_i\times b_j$（有时也写作 $(a*b)_n = \sum\limits_{i=0}^{n} a_i\times b_{n-i}$）。

例如，假设 $a = \{1,2,3\},b=\{3,2,1\},c=a*b$，可以得到 $c_0 = a_0\times b_0 = 3$，$c_1 = a_0\times b_1 + a_1\times b_0 = 8$，$c_2 = a_0\times b_2 + a_1\times b_1 + a_2\times b_0 = 14$，$\ldots$，最终得到 $c = \{3,8,14,8,3\}$。注意到我们的下标从 $0$ 开始，那是因为，如果此时我们作出 $a,b$ 和 $c$ 的生成函数（OGF），分别记作 $A(x),B(x)$ 与 $C(x)$，则我们会发现 $C(x)$ 恰好就是 $A(x)\times B(x)$。用上面的例子来解释即，$A(x) = 1+2x+3x^2,B(x) = 3+2x+x^2,A(x)\times B(x) = 3+8x+14x^2+8x^3+3x^4 = C(x)$。于是我们可以用接下来将要介绍的 FFT 进行数组上的卷积运算。

我们将数组的范围扩展到 $\Z$，就能得到卷积的离散形式：$(a*b)(x) = \sum\limits_{i} a_i\times b_{x-i}$。同理扩展到 $\R$，我们可以得到广义意义上的卷积形式：

$$\displaystyle{(f*g)(x) = \int_{-\infty}^{\infty}f(t)\times g(x-t)\,\mathrm{d}t}$$

FFT（Fast Fourier Transform，快速傅里叶变换）#

前置知识#

多项式#

众所周知，若干数字及字母的积称为单项式，而若干单项式的和称为多项式。本文中多项式特指形如 $\sum\limits_{i=0}^n a_i\times x^i$ 等的单变量的多项式。

同时，根据最高次项次数，称 $f(x) = \sum\limits_{i=0}^n a_i\times x^i$ 为一个 $n$ 次多项式，显然 $y=f(x)$ 就是一个 $n$ 次函数。

多项式的表示法#

多项式具有两种表示法：系数表示法与点值表示法。

对于多项式 $\sum\limits_{i=0}^n a_i\times x^i$，系数表示法就是用一个把每一项的系数按次数升序写在一起形成的一个数列 $\{a_0,a_1,\ldots,a_n\}$ 来表示这个多项式。例如 $\{1,2,3\}$ 就能够表示多项式 $1+2x+3x^2$，显然 $\{1,2,3\}$ 的生成函数即为 $f(x)=1+2x+3x^2$。有时，也把形如 $\sum\limits_{i=0}^n a_i\times x^i$ 的多项式直接称为系数表示法。

点值表示法，顾名思义，就是用 $(n+1)$ 个点的坐标来表示一个 $n$ 次函数代表的 $n$ 次多项式。例如：$(0,1)$ 和 $(2,5)$ 能够确定函数 $f(x) = 2x + 1$，即多项式 $1 + 2x$；$(-2,1),(1,4)$ 和 $(2,9)$ 能够确定函数 $f(x) = x^2 + 2x + 1$，即多项式 $1 + 2x + x^2$。可以证明这种用 $(n+1)$ 个点的坐标表示 $n$ 次多项式的方法是一一对应的。

多项式的运算#

多项式的加减法即将对应次数的项的系数项加减：$C(x)=A(x)\pm B(x)=\sum\limits_{i=0}^n (a_i\pm b_i) x^i$。若使用点值表示法，需要 $(n+1)$ 对横坐标相同的点，直接将纵坐标相加减即可。

多项式的乘法即将两个系数数列做卷积：$C(x)=A(x)\times B(x) = \sum\limits_{i=0}^n \sum\limits_{j=0}^n a_i\times b_j\times x^{i+j}$。注意此时需要至少 $(2n+1)$ 对横坐标相同的点。

复数#

众所周知，我们定义虚数单位 $i$，满足 $i^2 = -1$。形如 $z = a + bi(a,b\in \R)$ 的数称为复数，属于复数集 $C$。其中 $a,b$ 称为 $z$ 的实部和虚部。

数轴能够表示所有实数，对其进行扩充，增加一条垂直实轴的虚轴，形成类似平面直角坐标系的复平面，复平面中点 $(a,b)$ 对应了复数 $a+bi$。点 $(a,b)$ 到原点的距离 $\sqrt{a^2 + b^2}$ 称为复数 $z=a+bi$ 的模长，记作 $|z|$。以实轴为始边，$z$ 为终边的角 $\theta$ 称为 $z$ 的幅角，记作 $\operatorname{Arg} z$。一个复数的多个幅角之间相差 $2\pi$ 的整数倍，其中属于 $[0,2\pi)$ 的那个称为 $z$ 的辐角主值，记作 $\arg z$（一般讲到幅角通常就指辐角主值）。与 $z$ 关于实轴相互对称的复数称为 $z$ 的共轭复数，记作 $\bar z$。

同时，复数存在三角表示，形如 $z = r(\cos\theta + i\sin\theta)$，存在 $\begin{cases}r = |z|\\\theta = \operatorname{Arg} z\end{cases}$。

---

复数的四则运算与多项式无异。$(a+bi)\pm(c+di) = (a\pm c) + (b\pm d)i$，在复平面上表现为点的横纵坐标相加减，遵循向量加减的平行四边形法则与三角形法则。$(a+bi)\times(c+di)=(ac-bd)+(ad+bc)i$，在复平面上表现为模长相乘，幅角相加。

复数作指数时，欧拉公式指出 $e^{i\theta} = \cos\theta + i\sin\theta$，也就是说，$e^{i\theta}$ 在复平面上对应的复数模长为 $1$，幅角为 $\theta$。由此可以推出 $e^{i\pi} + 1 = 0$。

单位根#

在复平面内，和平面直角坐标系一样，可以以 $(0,0)$ 为圆心，$1$ 为半径，作出一个单位圆。

在复数范围内，$x^n=1$ 有 $n$ 个根，称为 $n$ 次单位根，记为 $\omega_n^1,\omega_n^2,\ldots,\omega_n^n$。一般称幅角大于 $0$ 的 $n$ 次单位根中幅角最小的那一个为 $\omega_n^1$，逆时针依次为 $\omega_n^2,\omega_n^3,\ldots,\omega_n^n$，可以发现 $\omega_n^n = 1$。根据复数乘法的定义，这些根的模都等于 $1$。进一步计算可以发现，$\omega_n^k = e^{i\cdot k\frac{2\pi}{n}},k\in \Z$，即，$n$ 次单位根对应的 $n$ 个点平分单位圆。证明：$(e^{i\cdot k\frac{2\pi}{n}})^n = e^{i\cdot 2k\pi} = \cos(2k\pi) + i\sin(2k\pi) = 1$，所以 $e^{i\cdot k\frac{2\pi}{n}}$ 是 $n$ 次单位根。又由于代数基本定理指出 $n$ 次方程有且仅有 $n$ 个根，于是这些就是该方程的所有根。

---

单位根具有很多很好的性质，例如：$n$ 为偶数时，记 $m = \dfrac{n}{2}$，有 $(\omega_n^k)^2 = \omega_m^k$ 和 $\omega_n^k = -\omega_{n}^{k+m}$。同时我们也可以循环定义单位根，即 $\omega_n^k = \omega_n^{k+Kn},K\in\Z$，这三条性质在 FFT 的过程中有重要作用。

FFT 过程#

省流：系数表示法 $\xrightarrow{DFT/FFT}$ 点值表示法 $\xrightarrow{IDFT/IFFT}$ 系数表示法。

简述#

在多项式乘法过程中，若使用系数表示法，复杂度是 $\mathcal{O}(n^2)$ 的；而若使用点值表示法，则复杂度是 $\mathcal{O}(n\log n)$ 的。那么我们当然倾向于使用点值表示法表示。但是，系数表示法 $\rightarrow$ 点值表示法（Discrete Fourier Transform，DFT，离散傅里叶变换）与点值表示法 $\rightarrow$ 系数表示法（Inverse Discrete Fourier Transform，IDFT，离散傅里叶逆变换）的过程都是 $\mathcal{O}(n^2)$ 的，于是我们分别使用 FFT 与 IFFT 来优化。

DFT & FFT#

显然暴力的 DFT 的复杂度是 $\mathcal{O}(n^2)$ 的，而 FFT 使用了分治的思想来优化 DFT 的过程。所谓分治，就是把大问题分成小问题求解，再递归求解小问题。在 DFT 的过程中，对于一个 $n-1$ 次的多项式（注意此后 $n$ 指多项式的次数加一，$m = \dfrac{n}{2}$），要取 $n$ 个点的点值，于是我们想到先把它分解为前 $m$ 个点和后 $m$ 个点。于是，若能根据前 $m$ 个点的点值 $\mathcal{O}(1)$ 求出后 $m$ 个点的点值，那么复杂度就有 $\mathcal{O}(n\log n)$ 的保证。又因为取点除不能重复外没有任何限制，所以我们考虑构造一种能够让我们事半功倍的取点方法。

既然想要找到这样一种 $\mathcal{O}(1)$ 的方式，必定离不开函数的对称性。但是大部分函数没有对称性，于是我们要将其配凑成具有对称性的形式。我们总可以某个多项式对应的函数 $f(x)$ 进行以下操作：

$$\begin{aligned} f(x)&=\sum\limits_{i=0}^{n-1}a_ix^i\\ &= \sum\limits_{i=0}^{m-1} a_{2i+1}x^{2i+1} + \sum\limits_{i=0}^{m-1} a_{2i}x^{2i}\\ &= x\times \sum\limits_{i=0}^{m-1} a_{2i+1}x^{2i} + \sum\limits_{i=0}^{m-1} a_{2i}x^{2i}\\ &= x\times \sum\limits_{i=0}^{m-1} a_{2i+1}(x^2)^i + \sum\limits_{i=0}^{m-1} a_{2i}(x^2)^i \end{aligned}$$

于是我们记 $g(x) = \sum\limits_{i=0}^{m-1} a_{2i+1}x^i,h(x) = \sum\limits_{i=0}^{m-1} a_{2i}x^i$，就有 $f(x) = h(x^2)+x\cdot g(x^2)$。此时我们发现，$f(-x)=h(x^2)-x\cdot g(x^2)$，只需要分别求出 $g(x^2)$ 与 $h(x^2)$ 即可一次性求出 $f(x)$ 与 $f(-x)$ 两个点。

此时总结一下，我们将多项式对应的 $f(x)$ 写成了 $h(x^2) + x\cdot g(x^2)$，故我们只需要选一组点 $\{x_1,x_2,\ldots,x_m\}(x_i > 0,i\in[1,m])$，求解的过程中因为得到了 $h(x_i^2)$ 与 $g(x_i^2)$，就能在算 $f(x_i)$ 时顺便算出 $f(-x_i)$，于是就可以优化时间复杂度。而且求解 $g(x_1^2),g(x_2^2),\ldots,g(x_m^2)$ 的过程还是求点值的过程，似乎可以分治解决了。

---

真的是这样吗？在 $f$ 的求点值过程中，我们可以任意取点，但是此时在 $g$ 的求解过程中，我们要求的是特定的几个值，所以不能递归求解。例如，我们选了求 $f(1),f(-1),f(2),f(-2)$，那么就要求 $g(1),g(4)$，而此时 $g(1),g(4)$ 是固定的而非我们自己选的，不满足 $1$ 和 $4$ 是相反数，所以不能再递归求解。

那么，既然实数不能这样递归，我们尝试使用复数。注意到 $n$ 次单位根 $\omega_n^k$ 与 $\omega_{n}^{k+m}$ 互为相反数的重要性质，我们考虑将 $n$ 次单位根作为需要的 $n$ 个点及点值，即要求 $f(\omega_n^1),f(\omega_n^2),\ldots,f(\omega_n^n)$，于是根据上面的分治规则，我们只需要求 $f(\omega_n^1),f(\omega_n^2),\ldots,f(\omega_n^m)$ 即可，因为 $\begin{cases}f(\omega_n^k) = h(\omega_m^k) + \omega_n^k\cdot g(\omega_m^k)\\f(\omega_n^{k+m}) = h(\omega_m^k) - \omega_n^k\cdot g(\omega_m^k)\end{cases}$，可以分治操作。而此时，递归到下一层要求的 $g(\omega_m^k)$ 刚好是所有 $m$ 次单位根，仍具有对称性，可以递归求解。这一步操作称为蝶形运算（蝴蝶操作）优化，即可以通过 $k$ 与 $k+m$ 处值得不断复写来得到答案。

因为分治递归的规模不断除以 $2$，而我们的公式只有在 $n$ 为偶数时才生效，于是我们需要将 $n$ 变为 $2$ 的正整数次幂，多余的项次数按照 $0$ 处理。

关于时间复杂度的证明：显然复杂度递推式为 $T(n) = 2T(\frac{n}{2}) + \mathcal{O}(n)$，根据主定理，其为 $\mathcal{O}(n\log n)$ 级别。

IDFT & IFFT#

现在只需要解决在 $\mathcal{O}(n\log n)$ 的时间内将点值表示法转换为系数表示法的问题，即，已知 $b_k = \sum\limits_{i=0}^{n-1} a_i\times (\omega_n^k)^i$，如何返回去解出 $a$ 数组。

先给出结论：$a_k = \dfrac{1}{n}\sum\limits_{i=0}^{n-1}b_i\times (\omega_n^{-k})^i$，接下来证明：

将 $b$ 数组代入：

$$\begin{aligned} \dfrac{1}{n}\sum\limits_{i=0}^{n-1}b_i\times (\omega_n^{-k})^i &= \dfrac{1}{n}\sum\limits_{i=0}^{n-1} (\sum\limits_{j=0}^{n-1} a_j\times \omega_n^{ij})\times\omega^{-ki}\\ &= \dfrac{1}{n}\sum\limits_{i=0}^{n-1} \sum\limits_{j=0}^{n-1} a_j\times \omega_n^{i(j-k)}\\ &= \dfrac{1}{n}\sum\limits_{j=0}^{n-1} a_j\times \sum\limits_{i=0}^{n-1} \omega_n^{i(j-k)} \end{aligned}$$

考虑到在 $j\neq k$ 时，$|j-k|\neq n$，于是 $\omega_n^{i(j-k)}\neq 1$，所以有 $\sum\limits_{i=0}^{n-1} \omega_n^{i(j-k)} = \dfrac{1-\omega_n^{n(j-k)}}{1-\omega_n^{j-k}} = \dfrac{1-\omega_n^n}{a-\omega_n^{j-k}} = \dfrac{1-1}{a-\omega_n^{j-k}} = 0$。所以 $\dfrac{1}{n}\sum\limits_{j=0}^{n-1} a_j\times \sum\limits_{i=0}^{n-1} \omega_n^{i(j-k)} = \dfrac{1}{n}a_k\times \sum\limits_{i=0}^{n-1}\omega_n^0 = a_k$，得证。

于是我们得到了 $a_k = \dfrac{1}{n}\sum\limits_{i=0}^{n-1}b_i\times (\omega_n^{-k})^i$，是为 DFT。

---

我们发现，在 DFT 的式子中，可以看作是另一个多项式 $\sum\limits_{i=0}^{n-1}b_i\times x^{-i}$ 在分别求 $n$ 次单位根处的点值，于是，我们可以再调用 FFT 来求解。注意到此处 $x$ 上方的指数由 $i$ 变为 $-i$，也就是说，FFT 求出的依次是 $\omega_n^0,\omega_n^1,\ldots,\omega_n^{n-1}$ 处的点值，而我们想要 $\omega_n^0,\omega_n^{-1},\ldots,\omega_n^{-(n-1)}$ 处的点值，即 $\omega_n^n,\omega_n^{n-1},\ldots,\omega_n^1$ 处的点值，只需要把枚举时乘 $\omega_n^1$ 改成乘 $\overline{\omega_n^1}$ 即可。

显然，这个过程的复杂度和 FFT 一样，是 $\mathcal{O}(n\log n)$ 的。

递归实现#

#include<iostream>
#include<cstdio>
#include<cmath>
#define maxn 2100005
using namespace std;
const double pi=3.1415926;
struct complex{
	double x,y; complex (double xx=0,double yy=0){x=xx; y=yy;} // 负数运算
	complex operator +(complex rhs){return complex(x+rhs.x,y+rhs.y);}
	complex operator -(complex rhs){return complex(x-rhs.x,y-rhs.y);}
	complex operator *(complex rhs){return complex(x*rhs.x-y*rhs.y,x*rhs.y+y*rhs.x);}
}a[maxn],b[maxn];
void fft(int N,complex *a,int sgn){
	if(N==1) return; complex a1[N],a2[N]; // 1 次单位根就是 1，于是 N=1 时直接返回
	for(int i=0;i<=(N>>1);i++) a1[i]=a[i<<1],a2[i]=a[i<<1|1]; // 按照奇偶分开
	fft(N>>1,a1,sgn); fft(N>>1,a2,sgn); // 分治求解
	complex e=complex(cos(pi*2/N),sgn*sin(pi*2/N)),z=complex(1,0); // e 即为 wn1，z 即为 wn0
	for(int i=0;i<(N>>1);i++,z=z*e) a[i]=a1[i]+z*a2[i],a[i+(N>>1)]=a1[i]-z*a2[i]; // 蝶形运算
}
int main(){
	int n,m; scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&a[i].x); for(int i=0;i<=m;i++) scanf("%lf",&b[i].x);
	int N=1; while(N<=n+m) N<<=1; fft(N,a,1); fft(N,b,1); // FFT
	for(int i=0;i<=N;i++) a[i]=a[i]*b[i]; fft(N,a,-1); // IFFT
	for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/N+0.5));
	return 0;
}

蝶形运算优化与位逆序置换#

递归的代码会存在爆栈等大问题，于是，我们需要进一步考虑优化。前面提到，蝶形运算优化基于迭代，于是，我们考虑用迭代来实现 FFT 的过程。

我们发现，以 $n=8$ 为例，我们可以写出各步中 $g(x)$ 和 $h(x)$ 的系数：

$n=8:\qquad a_0,a_1,a_2,a_3,a_4,a_5,a_6,a_7$
$n=4:\qquad a_0,a_2,a_4,a_6\qquad a_1,a_3,a_5,a_7$
$n=2:\qquad a_0,a_4\quad a_2,a_6\qquad a_1,a_5\quad a_3,a_7$
$n=1:\qquad a_0\quad a_4\quad a_2\quad a_6\quad\ a_1\quad a_5\quad a_3\quad a_7$

此时，我们观察下标：$0,4,2,6,1,5,3,7$，将他们写成二进制：$000,100,010,110,001,101,011,111$，再把每个二进制倒过来：$000,001,010,011,100,101,110,111$，刚好就是 $0,1,2,3,4,5,6,7$。可以证明，对于所有 $n$ 为 $2$ 的正整数次幂这个规律都成立。

于是我们可以在一开始就把 $a$ 按照这个顺序排列，然后直接对相邻的 $2$ 的正整数次幂个数进行蝶形运算即可。这一优化极大优化了空间复杂度与时间复杂度。

---

下面考虑如何求出这个下标顺序。找规律得到，$a_0 = 0$，$a_1 = a_0 | (n>>1)$，$\begin{cases}a_2 = a_0 | (n>>2)\\a_3 = a_1 | (n>>2)\end{cases}$，$\begin{cases}a_4 = a_0 | (n>>3)\\a_5 = a_1 | (n>>3)\\a_6 = a_2 | (n>>2)\\a_7 = a_3 | (n>>3)\end{cases}$。在 $n$ 为其他值时也是同样。

接下来考虑如何把数组按下标在 $a$ 中的顺序重组。记 $x$ 的二进制逆序对应的数为 $\operatorname{rev}(x)$，由于 $\operatorname{rev}(\operatorname{rev}(x))$，所以只需要把 $a_i$ 和 $a_{\operatorname{rev}(i)}$ 互换即可。

处理数组的时间复杂度显然是等于甚至远小于 $\mathcal{O}(n\log n)$ 的。这些操作称为位逆序置换。在位逆序置换后，就能迭代得进行蝶形运算优化。

迭代实现#

#include<iostream>
#include<cstdio>
#include<cmath>
#define maxn 4000005
using namespace std;
const double pi=3.14159265358;
int n,m,N=1,ord[maxn],cnt=0; struct complex{
	double x,y; complex (double xx=0,double yy=0){x=xx; y=yy;}
	complex operator +(complex rhs){return complex(x+rhs.x,y+rhs.y);}
	complex operator -(complex rhs){return complex(x-rhs.x,y-rhs.y);}
	complex operator *(complex rhs){return complex(x*rhs.x-y*rhs.y,x*rhs.y+y*rhs.x);}
}a[maxn],b[maxn];
void sorted(){
	ord[cnt++]=0; ord[cnt++]=N>>1; int thr=N>>2;
	for(int i=2;i<N;i<<=1,thr>>=1) for(int j=0;j<i;j++) ord[cnt++]=ord[j]|thr; // 求解下标数组
}
void fft(complex *a,int sgn){
	for(int i=1;i<N;i++) if(ord[i]>i) swap(a[i],a[ord[i]]); // 位逆序置换
	for(int len=2,gap=1;len<=N;gap<<=1,len<<=1){
		complex e=complex(cos(pi/gap),sgn*sin(pi/gap)),z=complex(1,0);
		for(int le=0,ri=len-1;ri<=N;le+=len,ri+=len){
			complex res=z; for(int i=le;i<le+gap;i++,res=res*e){
				complex z1=a[i]+res*a[i+gap],z2=a[i]-res*a[i+gap];
				a[i]=z1; a[i+gap]=z2;
			}
		}
	}
}
int main(){
	scanf("%d%d",&n,&m);
	for(int i=0;i<=n;i++) scanf("%lf",&a[i].x); for(int i=0;i<=m;i++) scanf("%lf",&b[i].x);
	while(N<=n+m) N<<=1; sorted(); fft(a,1); fft(b,1);
	for(int i=0;i<=N;i++) a[i]=a[i]*b[i]; fft(a,-1);
	for(int i=0;i<=n+m;i++) printf("%d ",(int)(a[i].x/N+0.5));
	return 0;
}

参考文章#

FFT·快速傅立叶变换
快速傅里叶变换
题解 P3803 【【模板】多项式乘法（FFT）】By attack
题解 P3803 【【模板】多项式乘法（FFT）】By 一扶苏一
OI-WIKI-FFT