蒙哥马利算法

本文章中，如果出现 $x^{-1}\pmod p$ 的形式，则我们认为 $x^{-1}$ 表示 $x$ 再模 $p$ 意义下的逆元。

共分为 3 Part，分别为蒙哥马利模乘，蒙哥马利约减，蒙哥马利模幂。

蒙哥马利模乘，简单地说，该算法能加快计算形如 $ab \pmod N$ 形式的计算。

• 蒙哥马利形式：为了计算 $ab \pmod N$，找一个 $R$，然后使得 $a’\equiv aR \pmod N,b’\equiv bR \pmod N$，其为蒙哥马利形式。

上述条件中 $R$ 并不是任意一个数，需要满足如下两个条件：

• $R=2^k>N$，$k$ 为满足条件的最小正数。

• $R\perp N$

为了计算 $ab\pmod N$，我们需要用到蒙哥马利形式。详细的，令 $X=a’b’$，设计一个函数计算 $XR^{-1}\pmod N$，不难得到 $X_1\equiv XR^{-1}\equiv abR\pmod N$ 的结果，那么我们再做一次 $X_1R^{-1}\pmod N$ 就能够得到 $ab\pmod N$ 的结果了。我们称这个算法为蒙哥马利约减。

不难看出蒙哥马利模乘的核心在于蒙哥马利约减，所以我们接着看如何解决蒙哥马利约减。

蒙哥马利约减是计算 $XR^{-1}\pmod N$，等价于 $\frac{X}{R}\mod N$，由前面 $R$ 的定义，我们可以直接运用位运算，将 $X$ 右移 $k$ 位。但是显然右移操作会自动下取整，我们不希望他出现下取整的操作，于是我们可以转而计算 $(X+mN)\cdot R^{-1}\pmod N$，使得 $R\mid (X+mN)$，那么我们就可以直接位运算了。

现在问题转向如何求 $m$，由于 $R \perp N$，由裴属定理可知我们必然能够找到 $N’$ 与 $R’$ 使得 $RR’-NN’=1$，于是便有 $N’=-N^{-1}\pmod R$。

列出式子，开始转移：

$$X+mN\equiv 0\pmod R$$ $$XN’+mNN’\equiv0\pmod R$$ $$XN’+m(RR’-1)\equiv 0\pmod R$$ $$XN’\equiv m\pmod R$$

转移过程是显然的，我们便这样求出了符合要求的 $m$。

上述蒙哥马利算法是用于加速模乘时的方案，实际上 $R$ 取任意大于 $N$ 的数都能够得到有意义的结果，下面我们给出一个例题。

[ARC148F] 998244353 → 1000000007

本题主要运用蒙哥马利约减，$XR^{-1}\equiv(X+mN)\cdot R^{-1}\pmod N$。

由于我们可以知道 $m\equiv -XN^{-1}\pmod R$，所以在知道 $N$ 和 $R$ 的前提下，我们只需要 $\bmod R$ 的操作就可以得到 $\bmod N$ 的结果。

于是在这个题中，我们令 $R=998244353$，$N=1000000007$，就能得到任意 $X \bmod N$ 的结果。

我们定义函数 $f(x)$ 为一次蒙哥马利约减后的结果，可以得到：

$$f(x)=(x+mN)\cdot R^{-1}$$

令 $x\leftarrow AB \bmod R$，我们有：

$$m<R$$ $$f(AB)<\lfloor\frac{AB}{R}\rfloor+N$$

故单次操作会把 $AB$ 映射到 $[1,\lfloor\frac{AB}{R}\rfloor+N]$ 的范围之内，那么我们进行 $\log_RAB$ 次操作后便能够把 $AB$ 映射到 $[1,N] $ 这个区间。

对于一个 $64$ 位整型变量 $x$，$\log_Rx$ 约等于 $3$，为了确保范围落在 $[1,N]$ 之间，我们进行 $4$ 次 $f$ 操作。

综上所述，我们便可以这样进行操作：

• $C\leftarrow AB\bmod R$
• $C\leftarrow f©$
• $C\leftarrow C \cdot (R^4 \bmod N)$
• $C\leftarrow f©$
• $C\leftarrow f©$
• $C\leftarrow f©$

最后我们便能够得到 $C=AB\bmod N$。