快速取模算法
Barret Reduction, Lemire Reduction 和 Montgomery Reduction.
问题
取模运算非常慢. 现代 C++ 编译器会自动优化对编译期常数的取模. 然而, 一些题目要求对运行期常量进行大量取模操作. 快速取模算法可以解决这类问题.
以下记这个运行期常量为
\[P \in [3, 2^{31}) \cap \mathbb{Z}.\]Barret Reduction 和 Lemire Reduction 解决 $x \bmod p$ 的问题. Montgomery Reduction 解决 $a_1a_2\cdots a_n \bmod P$ 的问题.
Barret Reduction
现代 CPU 中, 加减和位运算的时钟周期最短, 乘法次之, 除法和模运算显著慢于它们.
注意到
\[x \bmod P = x - \left\lfloor\frac{x}{P}\right\rfloor P.\]所以求解模运算只需要求解除法. 求解 $2^w$ 对应的除法非常快. Barret Reduction 的核心想法就是用
\[\left\lfloor\frac{xB}{2^w}\right\rfloor\]来近似
\[\left\lfloor\frac{x}{P}\right\rfloor.\]实际上, 这就是寻找
\[\frac{1}{P}\]在 $w$ 位二进制定点数格式下的表示.
令
\[B = \left\lfloor\frac{2^w}{P}\right\rfloor\]现在估计精度, 有
\[\begin{gather*} 0 \leq \frac{2^w}{P} - B < 1;\\ 0 \leq \frac{x}{P} - \frac{xB}{2^w} < \frac{x}{2^w}. \end{gather*}\]取 $w = 63$, 则大多数场景下有 $x < 2^w$, 所以
\[0 \leq \left\lfloor\frac{x}{P}\right\rfloor - \left\lfloor\frac{xB}{2^w}\right\rfloor \leq 1.\]使用时, 注意在关键的节点手动处理最后的一点精度.
这是工程上最典范的做法, Lemire Reduction 在大多数场景下弱于 Barret Reduction.
Lemire Reduction
不感兴趣可以跳过.
这是 Barret Reduction 的一种变体. Lemire Reduction 直接估计 $x \bmod P$ 而不是先求除法.
设
\[x = qP + r.\]其中 $r$ 就是所求的余数, 则
\[\begin{gather*} \frac{x}{P} = q + \frac{r}{P};\\ r = \left\{\frac{x}{P}\right\}P. \end{gather*}\]现在不妨模仿 Barret Reduction 的估计, 设
\[B = \left\lceil\frac{2^w}{P}\right\rceil.\](此处用上取整而非下取整可以简化下面的某处过程, 不需要特殊处理 0 和 P 的边界问题. 读者可以自行尝试.)
有
\[\left\{\frac{xB}{2^w}\right\} = \frac{xB \bmod 2^w}{2^w}.\]此时, 如果
\[\left\lfloor\frac{x}{P}\right\rfloor = \left\lfloor\frac{xB}{2^w}\right\rfloor,\]那么
\[0 \leq \frac{(xB \bmod 2^w) P}{2^w} - \left\{\frac{x}{P}\right\}P < \frac{xP}{2^w}.\]不同于 Barret Reduction, 此处
\[r = \left\{\frac{x}{P}\right\}P\]必定为整数, 因此只要
\[\frac{xP}{2^w} \leq 1\]就有
\[r = \left\lfloor\frac{(xB \bmod 2^w) P}{2^w}\right\rfloor.\]不存在精度问题.
如果
\[\left\lfloor\frac{x}{P}\right\rfloor = \left\lfloor\frac{xB}{2^w}\right\rfloor - 1,\]那么必然有
\[P - \frac{xP}{2^w} < r \leq P\]而由于 $r \in \mathbb{Z}$, 在刚才的条件下, 只能有 $r = P$. 这是不可能的.
为了解决 $x \in [0, 10^{18}], P \in [0, 10^9]$ 的问题, 一般地, Lemire Reduction 需要 $w = 90$, 这也是它的劣势.
Montgomery Reduction
先前的算法都利用了对 $2^w$ 作除法和取模运算非常快速的特点, Montgomery Reduction 也不例外. 特别地, 这个算法需要 $P$ 为奇数. 它可以快速地计算一系列数的乘积 ${}\bmod P$ 的值.
Montgomery 算法的核心是所谓的 REDC 操作, 下面首先介绍它.
REDC 操作
考虑取 $R = 2^w > P$, 由于 $P$ 是奇数, 因此 $P, R$ 互素.
任取
\[T \in [0, PR) \cap \mathbb{Z},\]首先用扩展欧几里德算法求解
\[P^{-1} \pmod R,\]随后令
\[Q = \left(-TP^{-1}\right) \bmod R.\]则有
\[T + QP \equiv 0 \pmod R.\]所以
\[S = \frac{T + QP}{R} \in [0, 2P) \cap \mathbb{Z}.\]且
\[S \equiv TR^{-1}\pmod P.\]要将 $S$ 放到 $[0, P)$ 中, 只需进行一次减法.
在模 $R$ 意义下 $P$ 的乘法逆可以预先求得, 以上过程只涉及关于 $R$ 的除法和取模, 可以用位运算快速完成.
以后记 $S = \mathrm{REDC}(T)$
求解连乘积
如果直接用 REDC 操作求 $ab \bmod P$, 那么每次会积累下一个 $R^{-1}$. 因此, Montgomery Reduction 首先对每个涉及的 $x$ 求解 Montgomery 表示
\[\overline{x} = (xR) \bmod P\]则有
\[\mathrm{REDC}\left(\overline{a} \cdot \overline{b}\right) = (aR \cdot bR \cdot R^{-1}) \bmod P = (abR) \bmod P = \overline{ab}.\]这样, 就可以用 REDC 操作求出
\[D = \prod_{i = 1}^n a_i\]的 Montgomery 表示 $\overline{D}$. 最后
\[D = \mathrm{REDC}\left(\overline{D}\right)\]即求解完毕.
现在还有一个问题: 如何求解
\[\overline{x} = xR \bmod P\]的值? 实际上, 通过预处理
\[G = R^2 \bmod P,\]则有
\[\overline{x} = \mathrm{REDC}(xG).\]于是解决.