快速傅里叶变换
点值表示法
在做多项式乘法时, 类似于竖式乘法, 一个朴素的想法是将一个多项式的每一项与另一个多项式的每一项相乘, 得到新的多项式, 时间复杂度 $\Theta(n^2)$。
我们考虑对多项式看作向量, 对其做线性变换, 其和系数表达是一一对应的, 从而可以在其变换上考虑。因此, 我们引入了点值表示。
初中我们就学过, 3 个点可以确定一个二次函数。将其扩展可以得到
定理
$n+1$ 个点可以确定一个 $n$ 阶多项式。
反证法。
假设存在两个不相同的 $n$ 阶多项式 $f(x),g(x)$, 使得它们有 $n+1$ 个交点, 则对于 $n$ 次函数 $h(x)=f(x)-g(x)$ 在 $x$ 轴上有 $n+1$ 个交点, 而由 $n$ 次方程恰有 $n$ 个根, 这是不可能的。故假设不成立, 即原命题得证。
因此, 我们可以用 $n+1$ 个点唯一表示一个多项式, 即
$$
f(x)=\{(x_0,y_0),(x_1,y_1),\dots,(x_n,y_n)\}
$$
我们将这种表示方法称作多项式的点值表示。
显然, 如果得出了两个多项式的点值表示法, 那么多项式的乘法将变的非常简单, 即将纵坐标相乘。至此, 问题转化成高效实现这个线性变换。我们将这种变换称为多项式的离散傅里叶变换。
但是, 我们确定了 $n+1$ 个点的横坐标时, 仍需要 $\Theta(n^2)$ 时间得到整个点集, 同时逆变换更是一个棘手的问题。因此, 我们引入了快速傅里叶变换。接下来介绍其中一种最常见的算法—库利—图基算法 (Cooley-Tukey)。
想要优化这个变换过程, 我们可以从减少选点的数量入手, 即期望选取更少的点, 同时这些点可以表达整个多项式。我们可以考虑, 选取了一个点 $(x,y)$ 后, 存在 $(-x, y)$ 或 $(-x, -y)$ 在这个多项式上, 那么选取的点的数量就可以减少了。
不难发现, 以上的性质对应于函数的奇偶性, 因此, 我们可以将这个多项式分解成奇偶两部分, 即
$$
f(x)=x\times g(x^2)+h(x^2)
$$
通过这种方式, 我们可以选取至多 $\dfrac{n}{2}$ 个点 $\{\pm x_0,\pm x_1,\dots,\pm x_{n/2-1}\}$, 并递归分治求解两个子问题, 向上合并时不难得到
$$
\forall i\in\left[0,\frac{n}{2}-1\right],
\begin{cases}
f(x_i)=x\times g(x_i^2)+h(x_i^2)\\
f(-x_i)=-x\times g(x_i^2)+h(x_i^2)
\end{cases}
$$
我们分析一下复杂度, 这个算法将问题分为了两个子问题, 并递归求解, 则有 $T(n)=2T\left(\dfrac n2\right)+\Theta(n)=\Theta(n\log n)$。至此, 这个算法可以在比较快的时间里解决多项式的离散傅里叶变换。
但是, 在分治过程中会存在一个问题, 即分治的子函数 $g,h$ 的横坐标均为正数, 不满足 $f$ 的横坐标互为相反数。而这就是 FFT 最精妙的思想—将这些数取成复数, 使得平方后仍然是互为相反数。
单位复根
在研究复数时, 我们通常与复平面结合起来研究问题。我们先来观察一个特例。
考虑 $n=4$, 先取横坐标为 1 的点, 则显然 $(-1)^2=1,1^2=1$, 再次向下递归时, 有 $i^2=-1,(-i)^2=-1$。不难发现, 这个过程相当于求解了一个方程 $x^n=1$。我们知道, 复数的乘法运算是“模长相乘, 辐角相加”, 则不难将其转移到单位圆上, 如图所示表示 $n=8$ 的情况。
There is a graph(TODO).
观察上图, 我们可以发现 $x^n=1$ 的解相当于将单位圆均分成 $n$ 份得到 $n$ 个点。
定义
$\omega_n$ 表示单位圆上 $n$ 个点中最小的正角, 则 $\omega_n=\cos \left(\dfrac{2\pi}{n}\right)+i\sin\left(\dfrac{2\pi}{n}\right)$, 将 $\omega_n$ 称为单位复根。
我们不难发现, $\forall k\in[0,n-1],\omega_n^k$ 相当于单位复根在单位圆上旋转, 其模长均为 1。
定理
$\omega_n^k(k\in[0,n-1])$ 是 $x^n=1$的 $n$ 个根。
证明
根据欧拉公式 (Euler Format) , $\omega_n=\cos \left(\dfrac{2\pi}{n}\right)+i\sin\left(\dfrac{2\pi}{n}\right)=e^{\frac{2\pi i}{n}}$。
则 $\left(\omega_n^k\right)^n=\left(\left(e^{\frac{2\pi i}{n}}\right)^k\right)^n=\left(e^{2\pi i}\right)^k=\left((-1)^2\right)^k=1^k=1$
再次观察上图我们发现, 两个根互为相反数, 则两个根在一条直线上, 即 $\omega_n^k=-\omega_n^{k+n/2}$, 分治过程中我们对一个根进行平方, 即 $\left(\omega_n^k\right)^2=\left(e^{\frac{4\pi i}{n}}\right)^k=\left(e^{\frac{2\pi i}{n/2}}\right)^k=\omega_{n/2}^k$。据此, 我们带入到上述递归式中, 得到
$$
\begin{aligned}
f\left(\omega_n^k\right)&=\omega_n^k\times g\left(\left(\omega_n^k\right)^2\right)+h\left(\left(\omega_n^k\right)^2\right)\\&=\omega_n^k\times g\left(\omega_{n/2}^k\right)+h\left(\omega_{n/2}^k\right)
\end{aligned}
$$
$$
\begin{aligned}
f\left(\omega_n^{k+n/2}\right)&=f\left(-\omega_n^k\right)\\&=-\omega_n^k\times g\left(\left(\omega_n^k\right)^2\right)+h\left(\left(\omega_n^k\right)^2\right)\\&=-\omega_n^k\times g\left(\omega_{n/2}^k\right)+h\left(\omega_{n/2}^k\right)
\end{aligned}
$$
至此, 我们通过引入复数, 成功解决了多项式的离散傅里叶变换。
快速傅里叶逆变换
在多项式乘法中, 我们还需要进行逆变换。而事实上, 快速傅里叶逆变换与快速傅里叶变换在形式上非常相像, 这里用线性代数及拉格朗日插值 (Lagrange polynomial) 来证明。
考虑分别用列向量 $\mathbf{A}=[a_0\ a_1\ a_2\ \cdots\ a_{n-1}]^{\intercal}, \mathbf{Y}=[y_0\ y_1\ y_2\ \cdots\ y_{n-1}]^{\intercal}$ 表示多项式及其离散傅里叶变换, 那么这个变换过程相当于 $\mathbf{A}$ 左乘一个 $n$ 阶方阵 $\mathbf{T}$, 有
$$
\begin{bmatrix}
y_0\\
y_1\\
y_2\\
\vdots\\
y_{n-1}
\end{bmatrix}=
\begin{bmatrix}
1 & 1 & 1 & \cdots & 1\\
1 & \omega_n & \omega_n^2 & \cdots & \omega_n^{n - 1}\\
1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n - 1)}\\
\vdots & \vdots & \vdots & \ddots & \vdots\\
1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n - 1)^2}
\end{bmatrix}
\begin{bmatrix}
a_0\\
a_1\\
a_2\\
\vdots\\
a_{n-1}
\end{bmatrix}
$$
上式可表达成 $\mathbf{Y}=\mathbf{T}\times \mathbf{A}$, 则有 $\mathbf{A}=\mathbf{T}^{-1}\times \mathbf{Y}$, 即逆变换相当于 $\mathbf{Y}$ 左乘 $n$ 阶方阵 $\mathbf{T}$ 的逆。其中 $\mathbf{T}$ 称为范德蒙矩阵 (Vandermonde matrix), 其有一个绝妙的性质。
定理
对于 $n$ 阶方阵 $\mathbf{T}=$
$$
\begin{bmatrix}1 & 1 & 1 & \cdots & 1\\1 & \omega_n & \omega_n^2 & \cdots & \omega_n^{n - 1}\\1 & \omega_n^2 & \omega_n^4 & \cdots & \omega_n^{2(n - 1)}\\\vdots & \vdots & \vdots & \ddots & \vdots\\1 & \omega_n^{n-1} & \omega_n^{2(n-1)} & \cdots & \omega_n^{(n - 1)^2}\end{bmatrix}
$$
其的逆 $\mathbf{T}^{-1}=$
$$
\frac{1}{n}\begin{bmatrix}1 & 1 & 1 & \cdots & 1\\1 & \omega_n^{-1} & \omega_n^{-2} & \cdots & \omega_n^{-(n - 1)}\\1 & \omega_n^{-2} & \omega_n^{-4} & \cdots & \omega_n^{-2(n - 1)}\\\vdots & \vdots & \vdots & \ddots & \vdots\\1 & \omega_n^{-(n-1)} & \omega_n^{-2(n-1)} & \cdots & \omega_n^{-(n - 1)^2}\end{bmatrix}
$$
即原方阵的每一项取倒数, 再除以 $n$。
证明
考虑这个逆变换在做什么。它将 $n$ 个点依次插值得到系数表达, 考虑拉格朗日插值基函数 (Lagrange basis polynomials)
$$
f(x)=\sum_{i=0}^{n-1}y_i\prod_{j\neq i}\frac{x-x_j}{x_i-x_j}
$$
分析可知,
$$
\mathbf{T}^{-1}[i,j]=[x^j]\prod_{k\neq i}\frac{x-x_k}{x_i-x_k}=[x^j]\frac{\prod\limits_{k\neq i}x-x_k}{\prod\limits_{k\neq i}x_i-x_k}
$$
设 $\lambda(x)=\prod\limits_{i=0}^{n-1}x-\omega_n^i=x^n-1$, 则根据洛必达法则 (L’Hôpital’s rule), 分母
$$
\prod_{k\neq i}x_i-x_k=\lim_{x\to \omega_n^i}\frac{\lambda(x)}{x-\omega_n^i}=\lim_{x\to \omega_n^i}\lambda’(x)=n\times \omega_n^{i(n-1)}=n\times \omega_n^{-i}
$$
根据泰勒级数 (Taylor series), 分子
$$
\begin{aligned}\prod_{k\neq i}x-x_j&=\frac{\lambda(x)}{x-\omega_n^i}=\frac{x^n-1}{x-\omega_n^i}=\omega_n^{-i}\left(\frac{x^n}{1-x\omega_n^{-i}}-\frac{1}{1-x\omega_n^{-i}}\right)\\&=\omega_n^{-i}\left(\sum_{j=0}^{\infty}\omega_n^{-i(j-n)}x^j-\sum_{j=0}^{\infty}\omega_n^{-ij}x^j\right)=\omega_n^{-i}\sum_{j=0}^{n-1}\omega_n^{-ij}x^j\end{aligned}
$$
则原式
$$
\mathbf{T}^{-1}[i,j]=[x^j]\frac{\omega_n^{-i}\sum_{j=0}^{n-1}\omega_n^{-ij}x^j}{n\times \omega_n^{-i}}=\frac{\omega_n^{-ij}}{n}
$$
至此, 这和快速傅里叶变换的差别仅仅只有一处, 即 $\omega_n$ 的取值。
快速数论变换
在快速傅里叶变换的实现过程中, 由于引入了三角函数, 使得存在大量的浮点数运算, 在数据规模较大时误差是难以避免的。快速数论变化 (Number-theoretic transform, NTT) 的范德蒙矩阵将全部取到整数, 即不牵涉浮点数运算, 确保了精度。
阶与原根
在快速傅里叶变换时, 我们在求解方程 $x^n=1$ 时引入了复数, 这里同理可以想到求解同余方程 $x^n\equiv 1\pmod p$, 据此我们可以自然地引入原根。
根据原根的性质, 对于原根 $g$, $g^t$ 每 $\varphi(p)$ 一周期, 故我们要满足 $n=2^m|\varphi(p)$, 考虑 $p$ 取 $q\times 2^m+1$ 且 $p$ 是质数, 常见的大质数有
$$
\begin{aligned}
p=998244353&=7\times 17\times 2^{23}+1\to g=3\\
p=1004535809&=479\times 2^{21}+1\to g=3
\end{aligned}
$$
同样, 类似于单位复根的概念, 定义 $\left(g_n^k\right)^n=1(k\in[0,n-1])$, 为了满足单位复根周期性的性质, 可以得到 $g_n=g^{\frac{p-1}{n}}$, 那么, 这里的 $g_n^k$ 可以对应于快速傅里叶变换中的 $\omega_n^k$。接下来我们证明 $g_n$ 与 $\omega_n$ 具有同样的性质。
定理
$g_n^{n/2}=-1,g_n^k=-g_n^{k+n/2},\left(g_n^k\right)^2=g_{n/2}^k$
证明
$\left(g_n^{n/2}\right)^2=g_n^n=1\implies g_n^{n/2}=\pm 1$。
根据阶的性质, $a^n=1\implies p-1|n$。
故 $g_n^{n/2}\neq 1\implies g_n^{n/2}=-1$。
$g_n^{k+n/2}=g_n^{n/2}\times g_n^k=-g_n^k$。$\left(g_n^k\right)^2=\left(g^{\frac{p-1}{n/2}}\right)^k=g_{n/2}^k$。
至此, 我们也只要将快速傅里叶变换的单位根的取值做出改动, 在 $p$ 的模意义运算即可得到快速数论变换。
应用
在算法竞赛中, 涉及线性变换的题目层出不穷, 下面通过一道经典例题来感受这一经典算法。
求 $n$ 个点的有标号无向连通图个数。$n\le 10^5$。
首先设 $f(n)$ 表示 $n$ 个点无向图的个数, 由于每条边可以连或不连, 则 $f(n)=2^{\frac{n(n-1)}{2}}$。设 $g(n)$ 表示 $n$ 个点的有标号无向连通图个数, 即答案。考虑枚举与 1 号节点相连的联通块大小, 剩下点可以任意连边, 同时还需乘上组合数, 则有
$$
f(n)=\sum_{i=1}^n \binom{n-1}{i-1}g(i)f(n-i)
$$
展开得到
$$
2^{\frac{n(n-1)}{2}}=\sum_{i=1}^n\frac{(n-1)!}{(i-1)!(n-i)!}g(i)\times 2^{\frac{(n-i)(n-i-1)}{2}}
$$
$$
\frac{g(n)}{(n-1)!}=\frac{2^{\frac{n(n-1)}{2}}}{(n-1)!}-\sum_{i=1}^{n-1}\frac{g(i)}{(i-1)!}\times \frac{2^{\frac{(n-i)(n-i-1)}{2}}}{(n-i)!}
$$
初始时 $g(1)=1$, 这样直接暴力转移可以做到 $\Theta(n^2)$, 由于这显然是一个卷积的性质, 可以用 FFT 优化, 但我们并不知道 $g(i)$ 的值, 考虑 CDQ 分治, 计算出前 $\dfrac n2$ 项对后 $\dfrac n2$ 项的贡献, 时间复杂度 $T(n)=2T\left(\dfrac n2\right)+\Theta(n\log n)=\Theta\left(n\log^2 n\right)$。
通过这道例题, 我们可以发现 FFT 相关算法可以作为优化计数题时间复杂度的重要手段, 同时也会与其他算法相结合。但事实上, 就笔者的经验来说, 推出递推式和化简才是计数题的难点和重中之重, FFT/NTT 仅能作为锦上添花的辅助手段。
声明：本文摘自本人学校由本人所作的探究性学习，作者均为本人，详情可 E-mail。