如何高效进行模乘、模幂运算？——蒙哥马利算法（Montgomery Algorithm）从入门到精通

最新推荐文章于 2021-12-15 00:02:38 发布

程序猿小宅

最新推荐文章于 2021-12-15 00:02:38 发布

阅读量1.8w

点赞数 48

分类专栏：算法学习文章标签：算法抽象代数密码学

本文链接：https://blog.csdn.net/a675115471/article/details/107553091

版权

算法学习专栏收录该内容

1 篇文章 1 订阅

订阅专栏

蒙哥马利算法（Montgomery Algorithm）从入门到精通

加密算法中，模运算（包括模乘、模幂运算）是难以避免的，如何高效地进行模运算，是提高算法效率的一个关键。

直观的想法

在数学上，模运算相当于是取余数的过程。以 $\div n = c \cdots\cdots d$ 为例，其中 $\leqslant d < n$ ，则称 $x$ 模 $n$ 等于 $d$ ，表示为 $\mod n = d$ ，并称 $n$ 为模数（modulus）。此外，若存在一个 $y$ ，使得 $\mod n = d$ 同时成立，则称 $x$ 和 $y$ 在模 $n$ 的情况下是同余的，并表示为 $\equiv y \ (\mod n)$ . 注意到，在一般的加密算法中，所有的操作数都是非负数，所以本文提到的模数均为正整数。

根据模运算的含义，以计算 $\mod b = d$ 为例，我们可以很容易地提出两种模运算的计算方法。

计算方法1

显然，存在一个 $c$ ，使得 $\div b = c \cdots\cdots d$ 成立，那么将 $a$ 减去 $\cdot c$ 即可得到 $d$ . $c$ 可以通过 $\lfloor a / b \rfloor$ 计算直接取得，这里的 $\lfloor \cdot \rfloor$ 表示的是向下取整的操作，比如 $\lfloor 7 / 4 \rfloor = \lfloor 1.75 \rfloor = 1$ . 该计算方法可直接表示为
$\cdot \lfloor a / b \rfloor \enspace .$

计算方法2

可以通过一个循环，从 $a$ 中不断地减去 $b$ ，直到结果落到区间 $[0, b)$ 内为止，该计算方法表示为伪代码

Input: $a$ , $b$

Output: $\mod b$

Algorithm 1:

$\leftarrow a$
while $d > = b$ do
$\qquad d \leftarrow d - b$
end
return $d$ .

评估

在计算机以及其他的硬件设备中，比起加法、乘法运算，除法运算的效率相当慢。故，计算方法1虽然表达简洁，但是效率不高。计算方法2利用减法操作，取代了除法运算，但是当 $a$ 和 $b$ 相差较大时，while循环次数将明显增多，显然也不是个高效的实现。

蒙哥马利算法

蒙哥马利算法解决的是模乘运算的效率问题，即给定模数 $n$ 和两个自然数 $a, b < n$ ，得到 $\cdot b \mod n$ 的值。

引入

先考虑 $\cdot b$ 如何运算。

假设在机器中所有的操作数均以 $r$ 进制的形式进行存储。比如，在传统的计算机中，操作数均为二进制的形式，则此时 $r = 2$ . 若模数 $n$ 可用 $k$ 位 $r$ 进制数表示，则 $a$ 和 $b$ 的位数也均不超过 $k$ 位。参考两个十进制数的竖式计算的方法，对 $b$ 按进制数分解，假设 $b$ 的 $k$ 位 $r$ 进制数形式表示为 $b_{k-1} b_{k-2} \dots b_0$ ，显然， $\sum^{k-1}_{i = 0}{b_i \cdot r^i}$ . 从而 $\cdot b$ 可表示为
$\cdot b := \sum^{k-1}_{i = 0}{a \cdot b_i \cdot r^i} \enspace .$
等号的右侧可以看做是关于 $r$ 的一个 $k - 1$ 次多项式，根据霍纳法则（Horner’s Rule）， $\cdot b$ 可进一步转化为
$\cdot b := \sum^{k-1}_{i = 0}{a \cdot b_i \cdot r^i} = (a \cdot b_0 + r \cdot (\cdots (a \cdot b_{k-2} + r \cdot (a \cdot b_{k-1})) \cdots ))$
的形式。显然，如果直接先计算完 $\cdot b$ 后，最后对乘积模 $n$ ，需要占用很多存储空间。好在对于乘法和加法，计算后统一取模的结果，与在计算过程中边取模边计算的结果是一样的，因此可以把模运算放在每次循环之中，从而减少资源的占用。伪代码表示为

Input: $a, b < n$

Output: $\cdot b \mod n$

Algorithm 2:

$\leftarrow 0$
for $i = k - 1$ downto $0$ step $- 1$ do
$\qquad d \leftarrow (a \cdot b_i + r \cdot d) \mod n$
end
return $d$ .

从伪代码来看，涉及到的运算全部集中于第3行中。而其中， $\cdot b_i$ 中 $b_i <r$ ，不难计算； $\cdot d$ 相当于在以 $r$ 进制表示的情况下，直接将 $d$ 左移1位，也很简单；然后将前后两个乘积相加，也不复杂；最后剩下的就是模 $n$ 计算了，又回到了一开始所提到的问题上来了。那么蒙哥马利是如何解决这一步的呢？

蒙哥马利模乘算法的实现

若将伪代码中第3行中的 $\cdot b_i + r \cdot d$ 看做一个整体，并将其记为 $z$ 的话，则问题等价于如何计算 $\mod n$ . 但是，直接计算 $\mod n$ 较为复杂，因此蒙哥马利取而代之的是计算 $\cdot r^{-1} \mod n$ . 若每次循环均计算的是 $\cdot r^{-1} \mod n$ ，那么 $k$ 次循环后，最终输出的就不是 $\cdot b \mod n$ 了，而是 $\cdot b \cdot r^{-k} \mod n$ .

由于
$\begin{aligned} a \cdot b \cdot r^{-k} &:= \sum^{k-1}_{i = 0}{a \cdot b_i \cdot r^i} \cdot r^{-k} \\ & = \sum^{k-1}_{i = 0}{a \cdot b_i \cdot r^{i-k}} \\ & = \sum^{k}_{i = 1}{a \cdot b_{k-i} \cdot {(r^{-1})}^i} \\ & = (a \cdot b_{k-1} + (\cdots(a \cdot b_1 + (a \cdot b_0) \cdot r^{-1}) \cdot r^{-1} \cdots ) \ ) \cdot r^{-1} \enspace ， \end{aligned}$
从而修改后的伪代码表示为

Input: $a, b < n$

Output: $\cdot b \cdot r^{-k} \mod n$

Algorithm 3:

$\leftarrow 0$
for $i = 0$ upto $k - 1$ step $+ 1$ do
$\qquad d \leftarrow (a \cdot b_i + d) \cdot r^{-1} \mod n$
end
return $d$ .

可以看到第3行已经由 $\mod n$ 转化成 $\cdot r^{-1} \mod n$ 的形式，这样的好处在于，若 $z$ 恰好是 $r$ 的倍数时，只需将 $z$ 右移1位即可，以此来取代模 $n$ 运算。但是，显然， $z$ 并不一定恰好为 $r$ 的倍数，因此，需要将 $z$ 加上一个值，在不改变模 $n$ 运算结果的前提下，得到 $r$ 的倍数。不改变模 $n$ 的运算结果，那么加上的这个值必为模数 $n$ 的倍数，不妨设其为 $\cdot n$ . 相加结果为 $r$ 的倍数，等效为
$\cdot n \mod r = 0 \enspace .$
若记 $z_0$ 为 $z$ 在 $r$ 进制表示下的最低位（例如十进制下的个位数、而二进制下的最低位），则当 $\cdot n \equiv -z_0 \mod r$ 时，等式成立。从而 $-z_0 \cdot n^{-1} \mod r$ ，而模 $r$ 是很容易计算的，当 $q$ 以 $r$ 进制表示时，模 $r$ 相当于取最低位，即直接取 $q_0$ . 那么 $-z_0 \cdot n^{-1} \mod r$ 中的 $n^{-1} \mod r$ 该如何计算呢？

欧拉定理告诉我们，当 $n$ 和 $r$ 互质的时候， $n^{\phi(r)-1} \equiv n^{-1} \mod r$ 成立，式中 $\phi(r)$ 称为欧拉函数，指的是比 $r$ 小的所有正整数中，与 $r$ 互质的数的总个数。通常在加密过程中，操作数都是二进制串，将两两一组看，可以看做是4进制；三个三个一组，可以看做是8进制；四个二进制四个二进制一组，可以看做是16进制，因此，前面提到的进制数 $r$ 可以表示为 $2^\nu$ ，对于 $2^\nu$ 进制而言，相当于就是将二进制串 $\nu$ 个 $\nu$ 个一组地看。而模数 $n$ ，通常都是取质数，我们知道，除了2以外的所有质数均为奇数，因此在实际操作过程中， $r$ 和 $n$ 显然是互质的，满足欧拉定理的前提条件。当 $2^\nu$ 时，有 $\phi(r) = \phi(2^\nu) = 2^\nu - 2^{\nu-1} = 2^{\nu-1}$ 的结论成立，因此只需计算 $n^{2^{\nu-1}-1} \mod r$ 即可得到 $n^{-1} \mod r$ 。

那么，该如何计算 $n^{2^{\nu-1} - 1} \mod r$ ？一种直观的想法是，将 $2^{\nu-1} - 1$ 个 $n$ 直接相乘。这是非常耗时的，试想一下，当操作数按十个十个一组来看，对应的 $\nu$ 为10，那么则需要计算510次乘法。快速模幂算法告诉我们，若即底数为 $n$ ，指数 $e$ 的二进制表示形式为 $e_{k-1}e_{k-2}\dots e_0$ ，模数为 $r$ ，那么快速计算 $n^e \mod r$ 的伪代码表示为

Input: $n$ , $e_{k-1}e_{k-2}\dots e_0$ , $r$

Output: $t = n^e \mod r$

Algorithm 4:

$\leftarrow 1$
for $i = k - 1$ downto $0$ step $- 1$ do
$\qquad t \leftarrow t \cdot t \mod r$
$\qquad$ if $e_i == 1$ then
$\qquad \qquad t \leftarrow t \cdot n \mod r$
$\qquad$ end
end
return t .

注意到， $2^{\nu-1} - 1$ 为 $\nu-1$ 位全为1的二进制数，因此，在计算 $n^{2^{\nu-1} - 1} \mod r$ 过程中，伪代码中的第4行恒成立，不需判断直接计算第5行。从快速模幂算法的伪代码我们可以看到，算法仅需执行 $2\log_2{e}$ 次乘法，仍是对于十个十个一组来看的情况来看的话，那么仅需计算18次乘法即可。

那么，蒙哥马利模乘算法基本可以实现了，借助快速模幂算法得到 $n^{-1} \mod r$ ，进而计算出每次循环时所要加上的 $\cdot n$ ，最终根据 Algorithm 3 就能计算得出 $\cdot b \cdot r^{-k} \mod n$ 这个最终结果了。记模数 $n$ 为 $k$ 位 $r$ 进制数，进制数 $2^\nu$ ，则伪代码表示为

Input: $a, b < n$

Output: $\cdot b \cdot r^{-k} \mod n$

Algorithm 5: $\mathrm{MontMul}$

$\qquad$ \\ 第1-5行，计算 $n^{-1} \mod r$

$\leftarrow 1$

for $i = 1$ upto $\nu-1$ step $+ 1$ do $\quad$ \\ 循环 $\nu-1$ 次

$\qquad t \leftarrow t \cdot t \cdot n \mod r$

end $\quad$ \\ 此时的 $t = n^{-1} \mod r$

$\leftarrow r - t \quad$ \\ 此时的 $t = -n^{-1} \mod r$

$\qquad$ \\ 第6-14行，计算 $\cdot b \cdot r^{-k} \mod n$

$\leftarrow 0$
for $i = 0$ upto $k - 1$ step $+ 1$ do
$\qquad q \leftarrow (a_0 \cdot b_i + d'_0) \cdot t \mod r \quad$ \\ $q = -z_0n^{-1} \mod r$ ，这里 $(a_0 \cdot b_i + d'_0)$ 相当于 $z_0$ ，而 $t$ 即是 $n^{-1} \mod r$
$\qquad d' \leftarrow (a \cdot b_i + d' + q\cdot n) \cdot r^{-1} \quad$ \\ 相较于伪代码 Algorithm 3，此时省略了 $\mod n$ 的操作，下面会进行解释
end $\quad$ \\ 此时的 $b^{'} < 2 n$
if $d^{'} > = n$ then $\quad$ \\ 相较于伪代码 Algorithm 3，此处新增了第11-13行，下面会进行解释
$\qquad d' \leftarrow d' - n$
end
return d’ .

在伪代码中，第3、8行均涉及到了 $\mod r$ 的操作，相当于是取 $r$ 进制的最低位，对于 $2^\nu$ 而言，即直接取二进制情况下的末 $\nu$ 位；第9行涉及到了 $⋅r−1 \cdot r^{-1}$ 的操作，相当于是在 $r$ 进制的形式下右移1位，对于 $2^\nu$ 而言，即直接右移 $\nu$ 位；而其他仅剩下简单的乘法、加法操作了。

注意到相较于 Algorithm 3，第9行，计算 $⋅r−1 \cdot r^{-1}$ 后直接省略了 $\mod n$ 操作。我们注意到一个事实，在第7-10行的循环中， $d^{'}$ 始终小于 $2 n$ ，即 $d^{'} < 2 n$ 恒成立，那么第10行退出循环时的 $d^{'}$ 无非就两种情况，一种是 $\in [0, n)$ ，一种是 $\in [n, 2n)$ 。对于前者， $d^{'}$ 无需再次模 $n$ ；对于后者，从 $d^{'}$ 中减去一个 $n$ 即是模 $n$ 的结果。所以在第9行中去掉 $\mod n$ 并不影响算法的正确性。下面我们对这个事实进行证明。

事实1 在 Algorithm 5 中，第7-10行中的 $d^{'}$ 始终不大于 $2 n$ ，即 $d^{'} < 2 n$ 恒成立。

事实1的证明 下面结合数学归纳法对事实1进行证明。

① 进入 for 循环前，第6行 $d^{'}$ 赋值为0，显然， $d^{'} < 2 n$ .

② 程序进入第7-10行的循环后，对于 $\cdots, k-1$ . 考察第9行，假设 $\leftarrow$ 右边的 $d^{'} < 2 n$ 成立。而根据输入， $\leqslant n -1 < n$ ，且 $b_i, q \leqslant r-1$ ，因此有
$\cdot b_i + d' + q\cdot n< n \cdot (r-1) + 2n + (r-1)\cdot n = 2nr .$
从而 $\cdot b_i + d' + q\cdot n) \cdot r^{-1} < 2nr \cdot r^{-1} = 2n$ ，经过第9行的赋值语句后，第 $i$ 次循环中得到的 $d^{'}$ 满足 $d^{'} < 2 n$ .

综上所述，由①②可知，直至第10行结束 for 循环， $d^{'} < 2 n$ 恒成立。

至此，蒙哥马利模乘算法解析完毕，我们记 Algorithm 5 为 $\mathrm{MontMul}$ ，调用为 $\leftarrow \mathrm{MontMul}(a, b)$ ，返回的 $\cdot b \cdot r^{-k} \mod n$ . 最后回忆一下，其中， $a, b$ 为乘数， $n$ 为模数， $a, b < n$ ， $r$ 为计算过程中采用的进制数， $k$ 为模数 $n$ 在 $r$ 进制表示下的位数，即 $r^{k-1} \leqslant n < r^k$ . 另外，通常取 $n$ 为质数，取 $2^\nu$ ，其中， $\nu$ 为正整数。注意到 $\mathrm{MontMul}$ 中只使用了移位、乘法、加法这三种计算机等硬件所擅长的操作，而消除了模运算、除法运算这两种计算机等硬件所不擅长的操作。此外，可以观察到，第1-5行求 $n^{-1} \mod r$ 的过程仅和 $n$ 与 $r$ 有关，因此可以事先计算好，且在同一场景下（模数 $n$ 相同时）重复使用。

回到问题的一开始，若要计算 $\cdot b \mod n$ ，则仅需调用 $\mathrm{MontMul}$ 4次即可。令 $\rho = r^{2k} \mod n$ ：

$\mathrm{MontMul}(a, \rho) = a \cdot r^k \mod n$ ，
$\mathrm{MontMul}(b, \rho) = b \cdot r^k \mod n$ ，
$\mathrm{MontMul}(A, B) = a \cdot b \cdot r^k \mod n$ ，
$\mathrm{MontMul}(D, 1) = a \cdot b \mod n$ .

第4步得到的 $d$ 即为所需结果。其中，1、2步称之为进入蒙哥马利域，即将操作数 $x$ 变为 $\cdot r^k \mod n$ 的形式；第4步称之为退出蒙哥马利域，也叫做蒙哥马利约减，即由 $\cdot r^k \mod n$ 的形式恢复出 $x$ . 另外， $\rho = r^{2\cdot k} \mod n$ 是个仅与 $r, n$ 有关的变量，可事先计算好，在同一场景下（模数 $n$ 相同时）重复使用。

蒙哥马利模幂算法的实现

实际上，蒙哥马利模幂运算（ $\mathrm{MontExp}$ ）可直接由快速模幂算法 Algorithm 4 结合蒙哥马利模乘算法 Algorithm 5 ( $\mathrm{MontMul}$ ) 稍作修改得到，伪代码表示为

Input: $e_{k-1}e_{k-2}\dots e_0$

Output: $t = a^e \mod n$

Algorithm 6: $\mathrm{MontExp}$

$\leftarrow \mathrm{MontMul}(1, \rho) \quad$ \\ 进入蒙哥马利域： $\cdot r^k \mod n$
$\leftarrow \mathrm{MontMul}(a, \rho) \quad$ \\ 进入蒙哥马利域： $\cdot r^k \mod n$
for $i = k - 1$ downto $0$ step $- 1$ do
$\qquad T \leftarrow \mathrm{MontMul}(T, T)$
$\qquad$ if $e_i == 1$ then
$\qquad \qquad T \leftarrow \mathrm{MontMul}(T, A)$
$\qquad$ end
end
return $\mathrm{MontMul}(T, 1) \quad$ \\ 退出蒙哥马利域 .

以上，即是关于蒙哥马利算法的全部内容。

程序猿小宅

关注

48
点赞
踩
148

收藏

觉得还不错? 一键收藏
26
评论
如何高效进行模乘、模幂运算？——蒙哥马利算法（Montgomery Algorithm）从入门到精通

蒙哥马利算法（Montgomery Algorithm）从入门到精通加密算法中，模运算（包括模乘、模幂运算）是难以避免的，如何高效地进行模运算，是提高算法效率的一个关键。直观的想法在数学上，模运算相当于是取余数的过程。以 x÷n=c⋯⋯dx \div n = c \cdots\cdots dx÷n=c⋯⋯d 为例，其中 0⩽d<n0 \leqslant d < n0⩽d<n，则称 xxx 模 nnn 等于 ddd，表示为 xmod n=dx \mod n = dxmodn
复制链接

扫一扫