Montgomery Modular Multiplication

Montgomery 模乘法是美国数学家 Peter L. Montgomery 于 1985 年在论文 Modular Multiplication Without Trial Division 中提出的，本文将介绍它。

不过，单个模乘法并不是 Montgomery 算法的应用场景，它的主要应用场景是模幂（Modular exponentiation）计算，即已知 $x,e,m$ 时计算 $x^{e}modm$ 。 当模数比较大时，使用 Montgomery 算法可以大大加快模幂的计算。

设 $R>m,gcd(R,m)=1$ 。数 $x$ 在 “Montgomery 空间”的表示可记为 $\tilde{x}$ ，它的定义是： $$\tilde{x}\triangleq xRmodm$$

设 $m=17,R=100$ ，显然满足 $gcd(R,m)=1$ 。表 1 展示了 $3,5,7,15$ 分别在 Montgomery 空间的对应表示。

Table 1: 原空间和 Montgomery 空间

原空间	Montgomery 空间
$3$	$3\times 100mod17=11$
$5$	$5\times 100mod17=7$
$7$	$7\times 100mod17=3$
$15$	$15\times 100mod17=4$

下面我们看看在 “Montgomery 空间” 进行模加法运算有什么特点。

$$\begin{array}{rl}\tilde{a}+\tilde{b}modm& =aR+bRmodm\\ & =(a+b)Rmodm& \text{这是 }a+b\text{ 在 Montgomery 空间的表示}\end{array}$$

从上式可知，要计算 $a+bmodm$ 除了直接计算外，还有一种在 Montgomery 空间进行计算的方法：

1. 把 $a,b$ 转换到 Montgomery 空间，分别记为 $\tilde{a},\tilde{b}$ ；
   
2. 在 Montgomery 空间进行模加法运算 $\tilde{a}+\tilde{b}modm$ ，假设结果是 $\tilde{x}$ ；
   
3. 把得到结果 $\tilde{x}$ 转换回原空间。
   

下面演示一下这个过程。假设要计算 $7+15mod17$ ，我们先把 $7$ 和 $15$ 转换到“Montgomery 空间”的数，分别得到 $3$ 和 $4$ ，然后计算 $3+4mod17=7$ ，最后把结果 $7$ 转换到原空间，得到结果 $5$ （参考表 1）。容易验证直接计算 $7+15mod17$ 的结果就是 $5$ 。

下面我们再看看在 “Montgomery 空间” 进行普通模乘法运算有什么特点。

$$\begin{array}{rl}\tilde{a}\cdot \tilde{b}modm& =(aRmodm)(bRmodm)modm\\ & =(abR)Rmodm& \text{这是 }abR\text{ 在 Montgomery 空间的表示}\end{array}$$

可以发现，要计算 $abmodm$ 时，如果转换到 “Montgomery 空间”进行计算，在计算 $\tilde{a}\cdot \tilde{b}modm$ 后，再转换回原空间得到的是 $abR$ ，而不是 $ab$ 。为了解决这个问题，我们定义 “Montgomery 空间” 的模乘运算为： $$\text{Mont}(\tilde{a},\tilde{b})\triangleq \tilde{a}\cdot \tilde{b}\cdot R^{-1}modm$$

在这种新模乘运算定义下，有： $$\begin{array}{rl}\text{Mont}(\tilde{a},\tilde{b})& \triangleq \tilde{a}\cdot \tilde{b}\cdot R^{-1}modm\\ & =(aRmodm)(bRmodm)R^{-1}modm\\ & =(abR)modm& \text{这是 }ab\text{ 在 Montgomery 空间的表示}\end{array}$$

这样，计算 $abmodm$ 有两种方法了，一种是原空间中直接计算，另一种转换到 Montgomery 空间后利用 $\text{Mont}$ 进行计算：

1. 把 $a,b$ 转换到 Montgomery 空间，分别记为 $\tilde{a},\tilde{b}$ ；
   
2. 在 Montgomery 空间计算 $\text{Mont}(\tilde{a},\tilde{b})$ ，假设结果是 $\tilde{x}$ ；
   
3. 把得到结果 $\tilde{x}$ 从 Montgomery 空间转换回原空间。
   

干嘛这么折腾呢？ 直接计算 $abmodm$ 比较耗时（求余过程中会涉及多次除法运算，而除法运算比较耗时），后面（节 4）将介绍一种方法可以快速地计算 $\text{Mont}(\tilde{a},\tilde{b})$ ，这样在 Montgomery 空间进行计算会比原空间更快。

不过，需要说明的是，如果只是进行一次模乘运算，转换到 Montgomery 空间进行计算没有优势，因为计算前需要把数字转换到 Montgomery 空间，在 Montgomery 空间得到结果后，还需要从 Montgomery 空间转换回原空间。如果是进行模幂运算（相当于多次模乘），则转移到 Montgomery 空间进行计算有很大优势，因为进入 Montgomery 空间，和退出 Montgomery 空间只需要在计算开始和结束时分别进行一次。

这节介绍 Montgomery Reduction 算法，它是下一节 4 的基础。

Montgomery Reduction 是指已知 $m,R,T$ ，且它们满足关系 $R>m,gcd(R,m)=1,0\le T<mR$ ，求解 $T{R}^{-1}modm$。

下面以具体值 $R=2^5=64$ 来介绍一下 $T{R}^{-1}modm$ 的求解。

1. 假设 $T$ 能够被 $R=2^5$ 整除，则计算 $T{R}^{-1}modm$ 非常快，因为这时有 $T{R}^{-1}modm=T>>5$ ，也就是说把 $T$ 右移 5 位就行了，不用进行除法运算了。由于 $T<mR$ ， $T$ 右移 5 位后会小于 $m$ ，不用再取模了。
   
2. 假设 $T$ 不能被 $R=2^5$ 整除，那怎么办呢？思路如下，我们想办法找这样的 $U$ ，使得 $T+Um$ 可以被 $R=2^5$ 整除，这样的话 $(T+Um)R^{-1}modm$ 就很容易求解了，因为和前面介绍的一样，它会等于 $T+Um>>5$ ，右移运算不涉及除法，是很快的。而 $T{R}^{-1}modm$ 显然和 $(T+Um)R^{-1}modm$ 是相等的，所以问题的关键就是寻找合适的 $U$ 。前面提到对 $U$ 的要求是 $T+Um$ 可以被 $R=2^5$ 整除，写成式子就是 $T+UmmodR=0$ ，从而有 $U=-T{m}^{-1}modR$ ，如果定义新记号 $m^{\prime }\triangleq -m^{-1}modR$ （一般我们可以提前把 $m^{\prime }$ 计算出来），则 $U=T{m}^{\prime }modR$ 。对 $U$ 的计算也是很快的，由于 $R=2^5$ ，取 $T{m}^{\prime }$ 的最低 5 个二进制位就行。
   

function REDC is
    input: Integers R and m with gcd(R, m) = 1,
           Integer m' in [0, R − 1] such that mm' ≡ −1 mod R,
           Integer T in the range [0, Rm − 1].
    output: Integer S in the range [0, m − 1] such that S ≡ TR⁻¹ mod m

    U ← ((T mod R)m') mod R        # mod R 的性能很高，编程实现时 R 一般是 2ⁿ，相当于只留下最低 n 个二进制位
    t ← (T + Um) / R               # 这里 T + Um 一定可以被 R 整除（前面有分析），编程实现是往往为右移运算
    if t ≥ m then                  # 因为 T < mR, U < R ，从而有 t = (T + Um) / R < (mR+mR)/R = 2m，所以 t ≥ m 时，变为 t-m 即可
        return t − m
    else
        return t
    end if
end function

在节 2 介绍中介绍了“Montgomery 空间” 的模乘运算 $\text{Mont}(x,y)\triangleq xy{R}^{-1}modm$ ，这节将介绍它的一个高效算法。

节 2 中提到过，如果只是进行一次模乘运算，转换到 Montgomery 空间进行计算没有优势，因为计算前需要把数字转换到 Montgomery 空间，在 Montgomery 空间得到结果后，还需要从 Montgomery 空间转换回原空间。 Montgomery 算法真正的应用场景是计算模幂，即 $x^{e}modm$ 。 因为进入 Montgomery 空间，和退出 Montgomery 空间只需要在计算开始和结束时分别进行一次。

在介绍使用 Mount(x,y) 快速计算模幂前，先介绍一下“Left-to-right binary exponentiation”算法，如图 3 所示。

下面以 $g^{283}$ 为例，介绍一下 Left-to-right binary exponentiation 算法的应用过程。

先把十进制 $283$ 表示为二进制，即 $(100011011{)}_2$ ，从而 $t=8$ ，按照算法的 Step 2 进行迭代，如表 4 所示。

Table 4: Left-to-right binary exponentiation 实例

$i$	$e_i$	A
8	1	$g$
7	0	$g^2$
6	0	$g^4$
5	0	$g^8$
4	1	$g^{17}$
3	1	$g^{35}$
2	0	$g^{70}$
1	1	$g^{141}$
0	1	$g^{283}$

注：Left-to-right binary exponentiation 并不是乘法次数最少的方法。比如求 $g^{15}$ 时，按照 Left-to-right binary exponentiation 算法需要 6 次乘法。我们还有更快的方法：先使用 2 次乘法求得 $g^3$ ，然后用 $g^3$ 自乘一次得 $g^6$ ，再用 $g^6$ 自乘一次得 $g^{12}$ ，最后把之前的 $g^3$ 和 $g^{12}$ 相乘可得最终结果 $g^{15}$ ，一共只使用了 5 次乘法。不过，这种方法没有 Left-to-right binary exponentiation 算法简单。

前面介绍的算法都是教科书中算法，在编程实现时，往往还有其它一些优化，论文 Analyzing and Comparing Montgomery Multiplication Algorithms 中分析和比较了几种实现。

OpenSSL 中 BN_mod_exp/BN_mod_exp_mont 的实现用到了 Montgomery 算法来加快模幂计算，参考：https://github.com/openssl/openssl/blob/9c8d04dbec03172d6ffe4eaa38ea4b1ac2741f26/crypto/bn/bn_exp.c#L304

1. 本文主要总结于 Handbook of Applied Cryptography，Chapter 14 Efficient Implementation
   
2. https://en.wikipedia.org/wiki/Montgomery_modular_multiplication
   
3. 1985 年 Peter L. Montgomery 发表的论文 Modular Multiplication Without Trial Division
   
4. 论文 Comparison of three modular reduction functions 中对传统方法、Barrett Reduction 算法和 Montgomery 算法进行了性能比较
   
5. Efficient Software Implementations of Modular Exponentiation
   
6. Rethinking Modular Multi-Exponentiation in Real-World Applications