蒙哥马利模乘算法简介

标签（空格分隔）： 蒙哥马利模乘

解决的问题：

快速求解大整数的模乘： 已 知 x < N , y < N ， 求 x y m o d    N 已知 x < N, y < N，求xy mod N 已知x

维基百科：https://en.wikipedia.org/wiki/Montgomery_modular_multiplication
csdn： https://blog.csdn.net/zgzczzw/article/details/52712980
上面两篇文章从原理上讲述了蒙哥马利模乘，本文仅在数学上进行推导。

欧拉函数

m>1时， φ ( m ) varphi(m) φ(m)为[1,m-1]中与m互质的数的个数；
φ ( 1 ) = 1 varphi(1)=1 φ(1)=1
对 于 > 1 的 质 数 p ， φ ( p ) = p − 1 , 对于>1的质数p， varphi(p)=p-1, 对于>1的质数p，φ(p)=p−1,
若 m = m 1 × m 2 , 且 g c d ( m 1 , m 2 ) = 1 m=m1times m2, 且 gcd(m1, m2)=1 m=m1×m2,且gcd(m1,m2)=1 则
φ ( m ) = φ ( m 1 ) × φ ( m 2 ) varphi(m) = varphi(m1)timesvarphi(m2) φ(m)=φ(m1)×φ(m2)

费马小定理

如果 p p p是一个质数，而整数 a a a不是p的倍数，则有 a φ ( p ) = a ( p − 1 ) ≡ 1 ( m o d p ) a^{varphi(p)} = a^{(p-1)} ≡1 pmod p aφ(p)=a(p−1)≡1(modp)
（证明参见维基百科）

模逆元

如果 a × a ′ ≡ 1 ( m o d N ) atimes a' equiv 1 pmod N a×a′≡1(modN)，则称 a ′ a' a′为 a a a 的模 N N N逆元

模运算法则

a + b ≡ ( a m o d    n ) + b ≡ a + ( b m o d    n ) ≡ ( ( a m o d    n ) + ( b m o d    n ) ) ( m o d n ) a+bequiv (a mod n) + b equiv a + (b mod n)equiv bigl((a mod n)+(bmod n)bigr) pmod n a+b≡(amodn)+b≡a+(bmodn)≡((amodn)+(bmodn))(modn)
a − b ≡ ( a m o d    n ) − b ≡ a − ( b m o d    n ) ≡ ( ( a m o d    n ) − ( b m o d    n ) ) ( m o d n ) a-bequiv (a mod n) - b equiv a - (b mod n)equiv bigl((a mod n)-(bmod n)bigr) pmod n a−b≡(amodn)−b≡a−(bmodn)≡((amodn)−(bmodn))(modn)
a b ≡ ( a m o d    n ) b ≡ a ( b m o d    n ) ≡ ( a m o d    n ) ( b m o d    n ) ( m o d n ) ab equiv (a mod n)b equiv a (bmod n) equiv (a mod n)(bmod n) pmod n ab≡(amodn)b≡a(bmodn)≡(amodn)(bmodn)(modn)

其他

已知： N ∈ Z 且 N > 1 ； R ∈ Z 且 R > N , ; g c d ( R , N ) = 1 N in Z 且N>1；R in Z且R>N,; gcd(R, N)=1 N∈Z且N>1；R∈Z且R>N,;gcd(R,N)=1
则 ∃ R ′ , N ′ exists R', N' ∃R′,N′,使下列成立：

R R ′ ≡ 1 ( m o d N ) RR'equiv 1 pmod N RR′≡1(modN)
R R ′ − N N ′ = 1 RR' - NN' = 1 RR′−NN′=1
R R ′ − N N ′ ≡ 1 ( m o d R ) RR'- NN' equiv 1 pmod R RR′−NN′≡1(modR)
− N N ′ ≡ 1 ( m o d R ) -NN' equiv 1 pmod R −NN′≡1(modR)
0 < R ′ < N 0 < R' < N 0 0 < N ′ < R 0 < N' < R 0 证明见文末尾证明

要计算 x y   m o d N xy mod N xy modN，需要先计算出 x y xy xy，再进行模运算。
传统的取模运算，需要通过除法来求得余数，而除法比较耗时。也可以通过反复做减法来实现，但是当商也是一个很大的数时，反复减法的次数过多，也会很耗时。
蒙哥马利模乘约减的思路是通过变换，将需要取模的数控制到很小的范围（由 [ 0 , N 2 − 2 N + 1 ] 变 为 [ 0 ， 2 N − 1 ] [0, N^2 - 2N + 1]变为 [0，2N-1] [0,N2−2N+1]变为[0，2N−1]，这样只需要通过最多一次减法即可完成取模运算；同时在变换过程中，通过选择除数，将除法的开销降到最小（比如在计算机中，除以2的幂通过移位即可实现，开销大大降低）。

考察如下性质：
令 T ∈ Z 且 T < R × N Tin Z 且 T< Rtimes N T∈Z且T 有： T = T × 1 = T ( R R ′ − T N N ′ ) = T R R ′ − T N N ′ T=T times1=T(RR' -TNN') = TRR'-TNN' T=T×1=T(RR′−TNN′)=TRR′−TNN′
∴ T + T N N ′ = T R R ′ therefore T+TNN'=TRR' ∴T+TNN′=TRR′
令 m = T N ′ m=T N' m=TN′
∴ T + m N = T R R ′ therefore T+ mN=TRR' ∴T+mN=TRR′
∴ ( T + m N ) / R = T R ′ therefore (T+ mN)/R=TR' ∴(T+mN)/R=TR′
∴ 对 任 意 T < N R , ∃ m 使 ( T + m N ) / R 是 整 数 therefore 对任意 T < NR, exist m使 (T+ mN)/R是整数 ∴对任意T

∴ T R ′ m o d    N therefore TR' mod N ∴TR′modN

= T + m N R m o d    N = frac{T+mN}R mod N =RT+mN​modN

= T + T N ′ N R m o d    N = frac{T + TN ' N}R mod N =RT+TN′N​modN

= T + ( T N ′ + k R ) N R m o d    N = frac{T + (TN ' + kR) N}R mod N =RT+(TN′+kR)N​modN

= T + ( ( ( T m o d    R ) N ′ ) m o d    R ) N R m o d    N =frac{T + Big(big((T mod R) N'big) mod RBig) N}R mod N =RT+(((TmodR)N′)modR)N​modN
即可以不计算 m = T N ′ m=TN' m=TN′，而是计算 m = ( ( T m o d    R ) N ′ ) m o d    R m=big((T mod R) N'big) mod R m=((TmodR)N′)modR

∴ ⇒ T R ′ m o d    N = T + m N R m o d    N ; m = ( ( T m o d    R ) N ′ ) m o d    R therefore Rightarrow TR' mod N = frac{T+mN}R mod N;m=big((T mod R) N'big) mod R ∴⇒TR′modN=RT+mN​modN;m=((TmodR)N′)modR

又 ∵ T < N R , m < R 又because T < NR, m < R 又∵T ∴ T + m N < N R + R N = 2 N R therefore T+mN < NR + RN = 2NR ∴T+mN ∴ T + m N R = T + ( ( ( T m o d    R ) N ′ ) m o d    R ) N R < 2 N therefore frac{T +mN}R= frac{T + Big(big((T mod R) N'big) mod RBig) N}R < 2 N ∴RT+mN​=RT+(((TmodR)N′)modR)N​<2N

∴ 最 后 求 m o d    N 的 计 算 量 很 小 ： therefore 最后求 mod N的计算量很小： ∴最后求modN的计算量很小：
令 t = T + m N R = T + ( ( ( T m o d    R ) N ′ ) m o d    R ) N R 令 t= frac{T +mN}R=frac{T + Big(big((T mod R) N'big) mod RBig) N}R 令t=RT+mN​=RT+(((TmodR)N′)modR)N​

if t >= N
	T mod N = t - N;
else 
	T mod N = t

同时，R选取满足gcd(R,N)=1且是2的幂的整数，除以R的运算转变为移位运算，对R求模的运算变为与运算。

所以通过选择 R R R以及上述变换，可以快速求取 T R ′ m o d N TR' mod N TR′modN。

上述变换，即为蒙哥马利约减算法Montgomery Reduction，记为REDC
即： R E D C ( T ) = T R ′ m o d    N REDC(T)=TR'mod N REDC(T)=TR′modN

REDC的算法伪代码见文末尾。

下面来看如何通过蒙哥马利约减计算蒙哥马利模乘。

已知： x ∈ Z 且 x < R N x in Z且x < RN x∈Z且x

求 ： z = a b m o d    N 求：z = abmod N 求：z=abmodN

z = a b m o d    N z=abmod N z=abmodN
= a b R R ′ m o d    N = abRR'mod N =abRR′modN
= R E D C ( a b R ) = REDC(abR) =REDC(abR)
= R E D C ( a b R R R ′ ) = REDC(abRRR') =REDC(abRRR′)
= R E D C ( R E D C ( a b R R ) ) = REDCbig(REDC(abRR)big) =REDC(REDC(abRR))
= R E D C ( R E D C ( ( a R m o d    N ) ( b R m o d    N ) ) ) = REDCBig(REDCbig((aRmod N) (bRmod N)big)Big) =REDC(REDC((aRmodN)(bRmodN)))
= R E D C ( R E D C ( ( a R R R ′ m o d    N ) ( b R R R ′ m o d    N ) ) ) = REDCBig(REDCbig((aRRR'mod N) (bRRR'mod N)big)Big) =REDC(REDC((aRRR′modN)(bRRR′modN)))
= R E D C ( R E D C ( R E D C ( a R R ) R E D C ( b R R ) ) ) =REDCBig(REDCbig(REDC(aRR) REDC(bRR)big)Big) =REDC(REDC(REDC(aRR)REDC(bRR)))
∵ R E D C ( a R R ) = a R R R ′ m o d    N = a ( R R m o d    N ) R ′ m o d    N = R E D C ( a ( R R m o d    N ) ) because REDC(aRR)=aRRR'mod N=a(RR mod N)R'mod N=REDC(a (RR mod N)) ∵REDC(aRR)=aRRR′modN=a(RRmodN)R′modN=REDC(a(RRmodN))
∴ z = R E D C ( R E D C ( R E D C ( a ( R R m o d    N ) ) R E D C ( b ( R R m o d    N ) ) ) ) therefore z=REDCBigg(REDCBig(REDCbig(a(RR mod N)big) REDCbig(b(RR mod N)big)Big)Bigg) ∴z=REDC(REDC(REDC(a(RRmodN))REDC(b(RRmodN))))

约定： x x x的蒙哥马利形式为： x ‾ = x R m o d    N overline x = xRmod N x=xRmodN
则上述过程可以简记为：

z = a b m o d    N = R E D C ( R E D C ( a ‾ b ‾ ) ) z=abmod N=REDCbig(REDC(overline a overline b)big) z=abmodN=REDC(REDC(ab))
a ‾ = R E D C ( a ( R R m o d    N ) ) overline a=REDCbig( a (RR mod N)big) a=REDC(a(RRmodN))
b ‾ = R E D C ( b ( R R m o d    N ) ) overline b=REDCbig( b (RR mod N)big) b=REDC(b(RRmodN))

可以看出，整个蒙哥马利模乘的过程中，只有求 R ′ , N ′ 以 及 R R m o d    N R', N' 以及RR mod N R′,N′以及RRmodN的计算无法简化。而通常 R 、 N R、N R、N都是预先选取好的常数，所以通常可以预计算 R ′ , N ′ , R R m o d    N R', N', RR mod N R′,N′,RRmodN

在RSA、ECC算法中，通常N都为一个大的奇数（更通常的为质数），所以R一般可以直接选取大于N的最小的2的幂。

前述简单描述了蒙哥马利模乘算法的原理和过程，但是在具体实现中，由于蒙哥马利约减和大整数乘法可以有多种组合方式；同时对于大整数的乘法和加法，通常需要分解为多个单字节的乘法、加法以及进位的循环，蒙哥马利约减与循环的组合也有多种方式。

论文《analyzing and comparing montgomery multiplication algorithms》详细论述了目前主流的5种算法的实现以及空时开销对比。

由于个人还未完全弄明白这几种算法的实现，所以本章内容暂时空缺。

(1)证明 ∃ R ′ , 使 R R ′ ≡ 1 ( m o d N ) 且 R ′ < N exists R', 使 RR'equiv 1 pmod N且R'< N ∃R′,使RR′≡1(modN)且R′

由费马小定理得知，如果 g c d ( R , N ) = 1 , 则 R φ ( N ) ≡ 1 ( m o d N ) gcd(R, N) = 1, 则 R^{varphi(N)} equiv 1pmod N gcd(R,N)=1,则Rφ(N)≡1(modN)
∴ R R φ ( N ) − 1 ≡ 1 ( m o d N ) therefore RR^{varphi(N) - 1} equiv 1pmod N ∴RRφ(N)−1≡1(modN)
∴ R 存 在 关 于 N 的 模 拟 元 ， 即 ∃ R ′ , 使 R R ′ ≡ 1 ( m o d N ) . therefore R存在关于N的模拟元，即exists R', 使 RR'equiv 1 pmod N. ∴R存在关于N的模拟元，即∃R′,使RR′≡1(modN).
可知R’不唯一，假设R’ >= N, 则：
R R ′ ≡ R ( k N + R ′ − k N ) ≡ R k N + R ( R ′ − k N ) ≡ R ( R ′ − k N ) ( m o d N ) RR' equiv R(kN+R'-kN) equiv RkN + R(R'-kN)equiv R(R'-kN)pmod N RR′≡R(kN+R′−kN)≡RkN+R(R′−kN)≡R(R′−kN)(modN)
∴ ∃ R ′ , 使 R R ′ ≡ 1 ( m o d N ) 且 R ′ < N therefore exists R', 使 RR'equiv 1 pmod N且R'< N ∴∃R′,使RR′≡1(modN)且R′

(2)证明 ∃ N ′ , 使 R R ′ − N N ′ = 1 , 且 N ′ < R exists N', 使 RR'-NN'=1,且N'

∵ R R ′ ≡ 1 ( m o d N ) because RR' equiv 1 pmod N ∵RR′≡1(modN)
∴ R R ′ = m N + 1 therefore RR' = mN + 1 ∴RR′=mN+1
∴ R R ′ − m N = 1 therefore RR' - mN =1 ∴RR′−mN=1
∴ ∃ N ′ = m 使 R R ′ − N N ′ = 1 ， 且 R R ′ − N N ′ ≡ 1 ( m o d R ) therefore exists N'= m 使 RR'-NN'=1 ，且 RR'-NN'equiv 1 pmod R ∴∃N′=m使RR′−NN′=1，且RR′−NN′≡1(modR)
∵ R R ′ = m N + 1 , R ′ < N because RR'=mN + 1, R' < N ∵RR′=mN+1,R′ ∴ R R ′ = m N + 1 < R N therefore RR' = mN + 1 < RN ∴RR′=mN+1 ∴ m N < R N therefore mN < RN ∴mN ∵ R > 0 because R > 0 ∵R>0
∴ N ′ = m < R therefore N'=m < R ∴N′=m

(3)证明 − N N ′ ≡ 1 ( m o d R ) -NN' equiv 1 pmod R −NN′≡1(modR)

∵ R R ′ − N N ′ ≡ 1 ( m o d R ) because RR'-NN' equiv 1 pmod R ∵RR′−NN′≡1(modR)
∴ − N N ′ ≡ 1 ( m o d R ) therefore -NN' equiv 1 pmod R ∴−NN′≡1(modR)

摘录于维基百科

function REDC is
    input: Integers R and N with gcd(R, N) = 1,
           Integer N′ in [0, R − 1] such that NN′ ≡ −1 mod R,
           Integer T in the range [0, RN − 1].
    output: Integer S in the range [0, N − 1] such that S ≡ TR−1 mod N

    m ← ((T mod R)N′) mod R
    t ← (T + mN) / R
    if t ≥ N then
        return t − N
    else
        return t
    end if
end function