蒙哥马利算法(Montgomery Algorithm)从入门到精通
加密算法中,模运算(包括模乘、模幂运算)是难以避免的,如何高效地进行模运算,是提高算法效率的一个关键。
直观的想法
在数学上,模运算相当于是取余数的过程。以 x ÷ n = c ⋯ ⋯ d x \div n = c \cdots\cdots d x÷n=c⋯⋯d 为例,其中 0 ⩽ d < n 0 \leqslant d < n 0⩽d<n, 则称 x x x 模 n n n 等于 d d d,表示为 x m o d n = d x \mod n = d xmodn=d,并称 n n n 为模数(modulus)。此外,若存在一个 y y y,使得 y m o d n = d y \mod n = d ymodn=d 同时成立,则称 x x x 和 y y y 在模 n n n 的情况下是同余的,并表示为 x ≡ y ( m o d n ) x \equiv y \ (\mod n) x≡y (modn) . 注意到,在一般的加密算法中,所有的操作数都是非负数,所以本文提到的模数均为正整数。
根据模运算的含义,以计算 a m o d b = d a \mod b = d amodb=d 为例,我们可以很容易地提出两种模运算的计算方法。
计算方法1
显然,存在一个
c
c
c,使得
a
÷
b
=
c
⋯
⋯
d
a \div b = c \cdots\cdots d
a÷b=c⋯⋯d 成立,那么将
a
a
a 减去
b
⋅
c
b \cdot c
b⋅c 即可得到
d
d
d.
c
c
c 可以通过
⌊
a
/
b
⌋
\lfloor a / b \rfloor
⌊a/b⌋ 计算直接取得,这里的
⌊
⋅
⌋
\lfloor \cdot \rfloor
⌊⋅⌋ 表示的是向下取整的操作,比如
⌊
7
/
4
⌋
=
⌊
1.75
⌋
=
1
\lfloor 7 / 4 \rfloor = \lfloor 1.75 \rfloor = 1
⌊7/4⌋=⌊1.75⌋=1 . 该计算方法可直接表示为
d
:
=
a
−
b
⋅
⌊
a
/
b
⌋
.
d := a - b \cdot \lfloor a / b \rfloor \enspace .
d:=a−b⋅⌊a/b⌋.
计算方法2
可以通过一个循环,从 a a a 中不断地减去 b b b,直到结果落到区间 [ 0 , b ) [0, b) [0,b) 内为止,该计算方法表示为伪代码
Input: a a a, b b b
Output: d = a m o d b d = a \mod b d=amodb
Algorithm 1:
- d ← a d \leftarrow a d←a
- while d > = b d >= b d>=b do
- d ← d − b \qquad d \leftarrow d - b d←d−b
- end
- return d d d .
评估
在计算机以及其他的硬件设备中,比起加法、乘法运算,除法运算的效率相当慢。故,计算方法1虽然表达简洁,但是效率不高。计算方法2利用减法操作,取代了除法运算,但是当 a a a 和 b b b 相差较大时,while循环次数将明显增多,显然也不是个高效的实现。
蒙哥马利算法
蒙哥马利算法解决的是模乘运算的效率问题,即给定模数 n n n 和两个自然数 a , b < n a, b < n a,b<n,得到 a ⋅ b m o d n a \cdot b \mod n a⋅bmodn 的值。
引入
先考虑 a ⋅ b a \cdot b a⋅b 如何运算。
假设在机器中所有的操作数均以
r
r
r 进制的形式进行存储。比如,在传统的计算机中,操作数均为二进制的形式,则此时
r
=
2
r = 2
r=2 . 若模数
n
n
n 可用
k
k
k 位
r
r
r 进制数表示,则
a
a
a 和
b
b
b 的位数也均不超过
k
k
k 位。参考两个十进制数的竖式计算的方法,对
b
b
b 按进制数分解,假设
b
b
b 的
k
k
k 位
r
r
r 进制数形式表示为
b
k
−
1
b
k
−
2
…
b
0
b_{k-1} b_{k-2} \dots b_0
bk−1bk−2…b0,显然,
b
=
∑
i
=
0
k
−
1
b
i
⋅
r
i
b = \sum^{k-1}_{i = 0}{b_i \cdot r^i}
b=∑i=0k−1bi⋅ri . 从而
a
⋅
b
a \cdot b
a⋅b 可表示为
a
⋅
b
:
=
∑
i
=
0
k
−
1
a
⋅
b
i
⋅
r
i
.
a \cdot b := \sum^{k-1}_{i = 0}{a \cdot b_i \cdot r^i} \enspace .
a⋅b:=i=0∑k−1a⋅bi⋅ri.
等号的右侧可以看做是关于
r
r
r 的一个
k
−
1
k - 1
k−1 次多项式,根据霍纳法则(Horner’s Rule),
a
⋅
b
a \cdot b
a⋅b 可进一步转化为
a
⋅
b
:
=
∑
i
=
0
k
−
1
a
⋅
b
i
⋅
r
i
=
(
a
⋅
b
0
+
r
⋅
(
⋯
(
a
⋅
b
k
−
2
+
r
⋅
(
a
⋅
b
k
−
1
)
)
⋯
)
)
a \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 ))
a⋅b:=i=0∑k−1a⋅bi⋅ri=(a⋅b0+r⋅(⋯(a⋅bk−2+r⋅(a⋅bk−1))⋯))
的形式。显然,如果直接先计算完
a
⋅
b
a \cdot b
a⋅b 后,最后对乘积模
n
n
n,需要占用很多存储空间。好在对于乘法和加法,计算后统一取模的结果,与在计算过程中边取模边计算的结果是一样的,因此可以把模运算放在每次循环之中,从而减少资源的占用。伪代码表示为
Input: a , b < n a, b < n a,b<n
Output: d = a ⋅ b m o d n d = a \cdot b \mod n d=a⋅bmodn
Algorithm 2:
- d ← 0 d \leftarrow 0 d←0
- for i = k − 1 i = k-1 i=k−1 downto 0 0 0 step − 1 -1 −1 do
- d ← ( a ⋅ b i + r ⋅ d ) m o d n \qquad d \leftarrow (a \cdot b_i + r \cdot d) \mod n d←(a⋅bi+r⋅d)modn
- end
- return d d d .
从伪代码来看,涉及到的运算全部集中于第3行中。而其中, a ⋅ b i a \cdot b_i a⋅bi 中 b i < r b_i <r bi<r,不难计算; r ⋅ d r \cdot d r⋅d 相当于在以 r r r 进制表示的情况下,直接将 d d d 左移1位,也很简单;然后将前后两个乘积相加,也不复杂;最后剩下的就是模 n n n 计算了,又回到了一开始所提到的问题上来了。那么蒙哥马利是如何解决这一步的呢?
蒙哥马利模乘算法的实现
若将伪代码中第3行中的 a ⋅ b i + r ⋅ d a \cdot b_i + r \cdot d a⋅bi+r⋅d 看做一个整体,并将其记为 z z z 的话,则问题等价于如何计算 z m o d n z \mod n zmodn . 但是,直接计算 z m o d n z \mod n zmodn 较为复杂,因此蒙哥马利取而代之的是计算 z ⋅ r − 1 m o d n z \cdot r^{-1} \mod n z⋅r−1modn . 若每次循环均计算的是 z ⋅ r − 1 m o d n z \cdot r^{-1} \mod n z⋅r−1modn,那么 k k k 次循环后,最终输出的就不是 d = a ⋅ b m o d n d = a \cdot b \mod n d=a⋅bmodn 了,而是 d ′ = a ⋅ b ⋅ r − k m o d n d' = a \cdot b \cdot r^{-k} \mod n d′=a⋅b⋅r−kmodn .
由于
a
⋅
b
⋅
r
−
k
:
=
∑
i
=
0
k
−
1
a
⋅
b
i
⋅
r
i
⋅
r
−
k
=
∑
i
=
0
k
−
1
a
⋅
b
i
⋅
r
i
−
k
=
∑
i
=
1
k
a
⋅
b
k
−
i
⋅
(
r
−
1
)
i
=
(
a
⋅
b
k
−
1
+
(
⋯
(
a
⋅
b
1
+
(
a
⋅
b
0
)
⋅
r
−
1
)
⋅
r
−
1
⋯
)
)
⋅
r
−
1
,
\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}
a⋅b⋅r−k:=i=0∑k−1a⋅bi⋅ri⋅r−k=i=0∑k−1a⋅bi⋅ri−k=i=1∑ka⋅bk−i⋅(r−1)i=(a⋅bk−1+(⋯(a⋅b1+(a⋅b0)⋅r−1)⋅r−1⋯) )⋅r−1,
从而修改后的伪代码表示为
Input: a , b < n a, b < n a,b<n
Output: d ′ = a ⋅ b ⋅ r − k m o d n d' = a \cdot b \cdot r^{-k} \mod n d′=a⋅b⋅r−kmodn
Algorithm 3:
- d ← 0 d \leftarrow 0 d←0
- for i = 0 i = 0 i=0 upto k − 1 k-1 k−1 step + 1 +1 +1 do
- d ← ( a ⋅ b i + d ) ⋅ r − 1 m o d n \qquad d \leftarrow (a \cdot b_i + d) \cdot r^{-1} \mod n d←(a⋅bi+d)⋅r−1modn
- end
- return d d d .
可以看到第3行已经由
z
m
o
d
n
z \mod n
zmodn 转化成
z
⋅
r
−
1
m
o
d
n
z \cdot r^{-1} \mod n
z⋅r−1modn 的形式,这样的好处在于,若
z
z
z 恰好是
r
r
r 的倍数时,只需将
z
z
z 右移1位即可,以此来取代模
n
n
n 运算。但是,显然,
z
z
z 并不一定恰好为
r
r
r 的倍数,因此,需要将
z
z
z 加上一个值,在不改变模
n
n
n 运算结果的前提下,得到
r
r
r 的倍数。不改变模
n
n
n 的运算结果,那么加上的这个值必为模数
n
n
n 的倍数,不妨设其为
q
⋅
n
q \cdot n
q⋅n . 相加结果为
r
r
r 的倍数,等效为
z
+
q
⋅
n
m
o
d
r
=
0
.
z + q \cdot n \mod r = 0 \enspace .
z+q⋅nmodr=0.
若记
z
0
z_0
z0 为
z
z
z 在
r
r
r 进制表示下的最低位(例如十进制下的个位数、而二进制下的最低位),则当
q
⋅
n
≡
−
z
0
m
o
d
r
q \cdot n \equiv -z_0 \mod r
q⋅n≡−z0modr 时,等式成立。从而
q
=
−
z
0
⋅
n
−
1
m
o
d
r
q = -z_0 \cdot n^{-1} \mod r
q=−z0⋅n−1modr,而模
r
r
r 是很容易计算的,当
q
q
q 以
r
r
r 进制表示时,模
r
r
r 相当于取最低位,即直接取
q
0
q_0
q0 . 那么
q
=
−
z
0
⋅
n
−
1
m
o
d
r
q = -z_0 \cdot n^{-1} \mod r
q=−z0⋅n−1modr 中的
n
−
1
m
o
d
r
n^{-1} \mod r
n−1modr 该如何计算呢?
欧拉定理告诉我们,当 n n n 和 r r r 互质的时候, n ϕ ( r ) − 1 ≡ n − 1 m o d r n^{\phi(r)-1} \equiv n^{-1} \mod r nϕ(r)−1≡n−1modr 成立,式中 ϕ ( r ) \phi(r) ϕ(r) 称为欧拉函数,指的是比 r r r 小的所有正整数中,与 r r r 互质的数的总个数。通常在加密过程中,操作数都是二进制串,将两两一组看,可以看做是4进制;三个三个一组,可以看做是8进制;四个二进制四个二进制一组,可以看做是16进制,因此,前面提到的进制数 r r r 可以表示为 2 ν 2^\nu 2ν,对于 2 ν 2^\nu 2ν 进制而言,相当于就是将二进制串 ν \nu ν 个 ν \nu ν 个一组地看。而模数 n n n,通常都是取质数,我们知道,除了2以外的所有质数均为奇数,因此在实际操作过程中, r r r 和 n n n 显然是互质的,满足欧拉定理的前提条件。当 r = 2 ν r = 2^\nu r=2ν 时,有 ϕ ( r ) = ϕ ( 2 ν ) = 2 ν − 2 ν − 1 = 2 ν − 1 \phi(r) = \phi(2^\nu) = 2^\nu - 2^{\nu-1} = 2^{\nu-1} ϕ(r)=ϕ(2ν)=2ν−2ν−1=2ν−1 的结论成立,因此只需计算 n 2 ν − 1 − 1 m o d r n^{2^{\nu-1}-1} \mod r n2ν−1−1modr 即可得到 n − 1 m o d r n^{-1} \mod r n−1modr。
那么,该如何计算 n 2 ν − 1 − 1 m o d r n^{2^{\nu-1} - 1} \mod r n2ν−1−1modr ?一种直观的想法是,将 2 ν − 1 − 1 2^{\nu-1} - 1 2ν−1−1 个 n n n 直接相乘。这是非常耗时的,试想一下,当操作数按十个十个一组来看,对应的 ν \nu ν 为10,那么则需要计算510次乘法。快速模幂算法告诉我们,若即底数为 n n n,指数 e e e 的二进制表示形式为 e k − 1 e k − 2 … e 0 e_{k-1}e_{k-2}\dots e_0 ek−1ek−2…e0,模数为 r r r,那么快速计算 n e m o d r n^e \mod r nemodr 的伪代码表示为
Input: n n n, e = e k − 1 e k − 2 … e 0 e = e_{k-1}e_{k-2}\dots e_0 e=ek−1ek−2…e0, r r r
Output: t = n e m o d r t = n^e \mod r t=nemodr
Algorithm 4:
- t ← 1 t \leftarrow 1 t←1
- for i = k − 1 i = k - 1 i=k−1 downto 0 0 0 step − 1 -1 −1 do
- t ← t ⋅ t m o d r \qquad t \leftarrow t \cdot t \mod r t←t⋅tmodr
- \qquad if e i = = 1 e_i == 1 ei==1 then
- t ← t ⋅ n m o d r \qquad \qquad t \leftarrow t \cdot n \mod r t←t⋅nmodr
- \qquad end
- end
- return t .
注意到, 2 ν − 1 − 1 2^{\nu-1} - 1 2ν−1−1 为 ν − 1 \nu-1 ν−1 位全为1的二进制数,因此,在计算 n 2 ν − 1 − 1 m o d r n^{2^{\nu-1} - 1} \mod r n2ν−1−1modr 过程中,伪代码中的第4行恒成立,不需判断直接计算第5行。从快速模幂算法的伪代码我们可以看到,算法仅需执行 2 log 2 e 2\log_2{e} 2log2e 次乘法,仍是对于十个十个一组来看的情况来看的话,那么仅需计算18次乘法即可。
那么,蒙哥马利模乘算法基本可以实现了,借助快速模幂算法得到 n − 1 m o d r n^{-1} \mod r n−1modr,进而计算出每次循环时所要加上的 q ⋅ n q \cdot n q⋅n,最终根据 Algorithm 3 就能计算得出 a ⋅ b ⋅ r − k m o d n a \cdot b \cdot r^{-k} \mod n a⋅b⋅r−kmodn 这个最终结果了。记模数 n n n 为 k k k 位 r r r 进制数,进制数 r = 2 ν r = 2^\nu r=2ν,则伪代码表示为
Input: a , b < n a, b < n a,b<n
Output: d ′ = a ⋅ b ⋅ r − k m o d n d' = a \cdot b \cdot r^{-k} \mod n d′=a⋅b⋅r−kmodn
Algorithm 5: M o n t M u l \mathrm{MontMul} MontMul
\qquad \\ 第1-5行,计算 − n − 1 m o d r -n^{-1} \mod r −n−1modr
t ← 1 t \leftarrow 1 t←1
for i = 1 i = 1 i=1 upto ν − 1 \nu-1 ν−1 step + 1 +1 +1 do \quad \\ 循环 ν − 1 \nu-1 ν−1 次
t ← t ⋅ t ⋅ n m o d r \qquad t \leftarrow t \cdot t \cdot n \mod r t←t⋅t⋅nmodr
end \quad \\ 此时的 t = n − 1 m o d r t = n^{-1} \mod r t=n−1modr
t ← r − t t \leftarrow r - t \quad t←r−t\\ 此时的 t = − n − 1 m o d r t = -n^{-1} \mod r t=−n−1modr
\qquad \\ 第6-14行,计算 d ′ = a ⋅ b ⋅ r − k m o d n d' = a \cdot b \cdot r^{-k} \mod n d′=a⋅b⋅r−kmodn
- d ′ ← 0 d' \leftarrow 0 d′←0
- for i = 0 i = 0 i=0 upto k − 1 k-1 k−1 step + 1 +1 +1 do
- q ← ( a 0 ⋅ b i + d 0 ′ ) ⋅ t m o d r \qquad q \leftarrow (a_0 \cdot b_i + d'_0) \cdot t \mod r \quad q←(a0⋅bi+d0′)⋅tmodr\\ q = − z 0 n − 1 m o d r q = -z_0n^{-1} \mod r q=−z0n−1modr, 这里 ( a 0 ⋅ b i + d 0 ′ ) (a_0 \cdot b_i + d'_0) (a0⋅bi+d0′) 相当于 z 0 z_0 z0,而 t t t 即是 − n − 1 m o d r -n^{-1} \mod r −n−1modr
- d ′ ← ( a ⋅ b i + d ′ + q ⋅ n ) ⋅ r − 1 \qquad d' \leftarrow (a \cdot b_i + d' + q\cdot n) \cdot r^{-1} \quad d′←(a⋅bi+d′+q⋅n)⋅r−1\\ 相较于伪代码 Algorithm 3,此时省略了 m o d n \mod n modn 的操作,下面会进行解释
- end \quad \\ 此时的 b ′ < 2 n b' < 2n b′<2n
- if d ′ > = n d' >= n d′>=n then \quad \\ 相较于伪代码 Algorithm 3,此处新增了第11-13行,下面会进行解释
- d ′ ← d ′ − n \qquad d' \leftarrow d' - n d′←d′−n
- end
- return d’ .
在伪代码中,第3、8行均涉及到了 m o d r \mod r modr 的操作,相当于是取 r r r 进制的最低位,对于 r = 2 ν r = 2^\nu r=2ν 而言,即直接取二进制情况下的末 ν \nu ν 位;第9行涉及到了 ⋅ r − 1 \cdot r^{-1} ⋅r−1 的操作,相当于是在 r r r 进制的形式下右移1位,对于 r = 2 ν r = 2^\nu r=2ν 而言,即直接右移 ν \nu ν 位;而其他仅剩下简单的乘法、加法操作了。
注意到相较于 Algorithm 3,第9行,计算 ⋅ r − 1 \cdot r^{-1} ⋅r−1 后直接省略了 m o d n \mod n modn 操作。我们注意到一个事实,在第7-10行的循环中, d ′ d' d′ 始终小于 2 n 2n 2n,即 d ′ < 2 n d' < 2n d′<2n 恒成立,那么第10行退出循环时的 d ′ d' d′ 无非就两种情况,一种是 d ′ ∈ [ 0 , n ) d' \in [0, n) d′∈[0,n),一种是 d ′ ∈ [ n , 2 n ) d' \in [n, 2n) d′∈[n,2n)。对于前者, d ′ d' d′ 无需再次模 n n n;对于后者,从 d ′ d' d′ 中减去一个 n n n 即是模 n n n 的结果。所以在第9行中去掉 m o d n \mod n modn 并不影响算法的正确性。下面我们对这个事实进行证明。
事实1 在 Algorithm 5 中,第7-10行中的 d ′ d' d′ 始终不大于 2 n 2n 2n,即 d ′ < 2 n d' < 2n d′<2n 恒成立。
事实1的证明 下面结合数学归纳法对事实1进行证明。
① 进入 for 循环前,第6行 d ′ d' d′ 赋值为0,显然, d ′ < 2 n d' < 2n d′<2n .
② 程序进入第7-10行的循环后,对于
i
=
0
,
1
,
⋯
,
k
−
1
i = 0, 1, \cdots, k-1
i=0,1,⋯,k−1 . 考察第9行,假设
←
\leftarrow
← 右边的
d
′
<
2
n
d' < 2n
d′<2n 成立。而根据输入,
a
⩽
n
−
1
<
n
a \leqslant n -1 < n
a⩽n−1<n,且
b
i
,
q
⩽
r
−
1
b_i, q \leqslant r-1
bi,q⩽r−1,因此有
a
⋅
b
i
+
d
′
+
q
⋅
n
<
n
⋅
(
r
−
1
)
+
2
n
+
(
r
−
1
)
⋅
n
=
2
n
r
.
a \cdot b_i + d' + q\cdot n< n \cdot (r-1) + 2n + (r-1)\cdot n = 2nr .
a⋅bi+d′+q⋅n<n⋅(r−1)+2n+(r−1)⋅n=2nr.
从而
(
a
⋅
b
i
+
d
′
+
q
⋅
n
)
⋅
r
−
1
<
2
n
r
⋅
r
−
1
=
2
n
(a \cdot b_i + d' + q\cdot n) \cdot r^{-1} < 2nr \cdot r^{-1} = 2n
(a⋅bi+d′+q⋅n)⋅r−1<2nr⋅r−1=2n,经过第9行的赋值语句后,第
i
i
i 次循环中得到的
d
′
d'
d′ 满足
d
′
<
2
n
d' < 2n
d′<2n .
综上所述,由①②可知,直至第10行结束 for 循环, d ′ < 2 n d' < 2n d′<2n 恒成立。
至此,蒙哥马利模乘算法解析完毕,我们记 Algorithm 5 为 M o n t M u l \mathrm{MontMul} MontMul,调用为 d ′ ← M o n t M u l ( a , b ) d' \leftarrow \mathrm{MontMul}(a, b) d′←MontMul(a,b),返回的 d ′ = a ⋅ b ⋅ r − k m o d n d' = a \cdot b \cdot r^{-k} \mod n d′=a⋅b⋅r−kmodn . 最后回忆一下,其中, a , b a, b a,b 为乘数, n n n 为模数, a , b < n a, b < n a,b<n, r r r 为计算过程中采用的进制数, k k k 为模数 n n n 在 r r r 进制表示下的位数,即 r k − 1 ⩽ n < r k r^{k-1} \leqslant n < r^k rk−1⩽n<rk . 另外,通常取 n n n 为质数,取 r = 2 ν r = 2^\nu r=2ν,其中, ν \nu ν 为正整数。注意到 M o n t M u l \mathrm{MontMul} MontMul 中只使用了移位、乘法、加法这三种计算机等硬件所擅长的操作,而消除了模运算、除法运算这两种计算机等硬件所不擅长的操作。此外,可以观察到,第1-5行求 − n − 1 m o d r -n^{-1} \mod r −n−1modr 的过程仅和 n n n 与 r r r 有关,因此可以事先计算好,且在同一场景下(模数 n n n 相同时)重复使用。
回到问题的一开始,若要计算 a ⋅ b m o d n a \cdot b \mod n a⋅bmodn,则仅需调用 M o n t M u l \mathrm{MontMul} MontMul 4次即可。令 ρ = r 2 k m o d n \rho = r^{2k} \mod n ρ=r2kmodn:
- A = M o n t M u l ( a , ρ ) = a ⋅ r k m o d n A = \mathrm{MontMul}(a, \rho) = a \cdot r^k \mod n A=MontMul(a,ρ)=a⋅rkmodn,
- B = M o n t M u l ( b , ρ ) = b ⋅ r k m o d n B = \mathrm{MontMul}(b, \rho) = b \cdot r^k \mod n B=MontMul(b,ρ)=b⋅rkmodn,
- D = M o n t M u l ( A , B ) = a ⋅ b ⋅ r k m o d n D = \mathrm{MontMul}(A, B) = a \cdot b \cdot r^k \mod n D=MontMul(A,B)=a⋅b⋅rkmodn,
- d = M o n t M u l ( D , 1 ) = a ⋅ b m o d n d = \mathrm{MontMul}(D, 1) = a \cdot b \mod n d=MontMul(D,1)=a⋅bmodn .
第4步得到的 d d d 即为所需结果。其中,1、2步称之为进入蒙哥马利域,即将操作数 x x x 变为 x ⋅ r k m o d n x \cdot r^k \mod n x⋅rkmodn 的形式;第4步称之为退出蒙哥马利域,也叫做蒙哥马利约减,即由 x ⋅ r k m o d n x \cdot r^k \mod n x⋅rkmodn 的形式恢复出 x x x . 另外, ρ = r 2 ⋅ k m o d n \rho = r^{2\cdot k} \mod n ρ=r2⋅kmodn 是个仅与 r , n r, n r,n 有关的变量,可事先计算好,在同一场景下(模数 n n n 相同时)重复使用。
蒙哥马利模幂算法的实现
实际上,蒙哥马利模幂运算( M o n t E x p \mathrm{MontExp} MontExp)可直接由快速模幂算法 Algorithm 4 结合蒙哥马利模乘算法 Algorithm 5 ( M o n t M u l \mathrm{MontMul} MontMul) 稍作修改得到,伪代码表示为
Input: a , e = e k − 1 e k − 2 … e 0 a, e = e_{k-1}e_{k-2}\dots e_0 a,e=ek−1ek−2…e0
Output: t = a e m o d n t = a^e \mod n t=aemodn
Algorithm 6: M o n t E x p \mathrm{MontExp} MontExp
- T ← M o n t M u l ( 1 , ρ ) T \leftarrow \mathrm{MontMul}(1, \rho) \quad T←MontMul(1,ρ)\\ 进入蒙哥马利域: T = 1 ⋅ r k m o d n T = 1 \cdot r^k \mod n T=1⋅rkmodn
- A ← M o n t M u l ( a , ρ ) A \leftarrow \mathrm{MontMul}(a, \rho) \quad A←MontMul(a,ρ)\\ 进入蒙哥马利域: A = a ⋅ r k m o d n A = a \cdot r^k \mod n A=a⋅rkmodn
- for i = k − 1 i = k - 1 i=k−1 downto 0 0 0 step − 1 -1 −1 do
- T ← M o n t M u l ( T , T ) \qquad T \leftarrow \mathrm{MontMul}(T, T) T←MontMul(T,T)
- \qquad if e i = = 1 e_i == 1 ei==1 then
- T ← M o n t M u l ( T , A ) \qquad \qquad T \leftarrow \mathrm{MontMul}(T, A) T←MontMul(T,A)
- \qquad end
- end
- return M o n t M u l ( T , 1 ) \mathrm{MontMul}(T, 1) \quad MontMul(T,1)\\ 退出蒙哥马利域 .
以上,即是关于蒙哥马利算法的全部内容。