Montgomery Modular Multiplication
Table of Contents
1. 简介
Montgomery 模乘法是美国数学家 Peter L. Montgomery 于 1985 年在论文 Modular Multiplication Without Trial Division 中提出的,本文将介绍它。
不过,单个模乘法并不是 Montgomery 算法的应用场景,它的主要应用场景是模幂(Modular exponentiation)计算,即已知 \(x,e,m\) 时计算 \(x^e \bmod m\) 。 当模数比较大时,使用 Montgomery 算法可以大大加快模幂的计算。
2. Montgomery 空间
设 \(R > m, \gcd(R,m) = 1\) 。数 \(x\) 在 “Montgomery 空间”的表示可记为 \(\widetilde{x}\) ,它的定义是: \[\widetilde{x} \triangleq xR \bmod m\]
设 \(m = 17, R = 100\) ,显然满足 \(\gcd(R,m)=1\) 。表 1 展示了 \(3,5,7,15\) 分别在 Montgomery 空间的对应表示。
原空间 | Montgomery 空间 |
---|---|
\(3\) | \(3 \times 100 \bmod 17 = 11\) |
\(5\) | \(5 \times 100 \bmod 17 = 7\) |
\(7\) | \(7 \times 100 \bmod 17 = 3\) |
\(15\) | \(15 \times 100 \bmod 17 = 4\) |
下面我们看看在 “Montgomery 空间” 进行模加法运算有什么特点。
\[\begin{aligned} \widetilde{a} + \widetilde{b} \bmod m &= a R + b R \bmod m \\
&= (a+b) R \bmod m & \text{这是 $a+b$ 在 Montgomery 空间的表示} \end{aligned}\]
从上式可知,要计算 \(a+b \bmod m\) 除了直接计算外,还有一种在 Montgomery 空间进行计算的方法:
- 把 \(a,b\) 转换到 Montgomery 空间,分别记为 \(\widetilde{a},\widetilde{b}\) ;
- 在 Montgomery 空间进行模加法运算 \(\widetilde{a} + \widetilde{b} \bmod m\) ,假设结果是 \(\widetilde{x}\) ;
- 把得到结果 \(\widetilde{x}\) 转换回原空间。
下面演示一下这个过程。假设要计算 \(7+15 \bmod 17\) ,我们先把 \(7\) 和 \(15\) 转换到“Montgomery 空间”的数,分别得到 \(3\) 和 \(4\) ,然后计算 \(3+4 \bmod 17 = 7\) ,最后把结果 \(7\) 转换到原空间,得到结果 \(5\) (参考表 1)。容易验证直接计算 \(7+15 \bmod 17\) 的结果就是 \(5\) 。
下面我们再看看在 “Montgomery 空间” 进行普通模乘法运算有什么特点。
\[\begin{aligned} \widetilde{a} \cdot \widetilde{b} \bmod m &= (a R \bmod m) (b R \bmod m) \bmod m \\
&= (abR)R \bmod m & \text{这是 $abR$ 在 Montgomery 空间的表示} \end{aligned}\]
可以发现,要计算 \(ab \bmod m\) 时,如果转换到 “Montgomery 空间”进行计算,在计算 \(\widetilde{a} \cdot \widetilde{b} \bmod m\) 后,再转换回原空间得到的是 \(abR\) ,而不是 \(ab\) 。为了解决这个问题,我们定义 “Montgomery 空间” 的模乘运算为: \[\text{Mont}(\widetilde{a},\widetilde{b}) \triangleq \widetilde{a} \cdot \widetilde{b} \cdot R^{-1}\bmod m\]
在这种新模乘运算定义下,有: \[\begin{aligned} \text{Mont}(\widetilde{a},\widetilde{b}) &\triangleq \widetilde{a} \cdot \widetilde{b} \cdot R^{-1}\bmod m \\
&=(a R \bmod m) (b R \bmod m) R^{-1} \bmod m \\
&=(abR) \bmod m & \text{这是 $ab$ 在 Montgomery 空间的表示} \end{aligned}\]
这样,计算 \(ab \bmod m\) 有两种方法了,一种是原空间中直接计算,另一种转换到 Montgomery 空间后利用 \(\text{Mont}\) 进行计算:
- 把 \(a,b\) 转换到 Montgomery 空间,分别记为 \(\widetilde{a},\widetilde{b}\) ;
- 在 Montgomery 空间计算 \(\text{Mont}(\widetilde{a},\widetilde{b})\) ,假设结果是 \(\widetilde{x}\) ;
- 把得到结果 \(\widetilde{x}\) 从 Montgomery 空间转换回原空间。
干嘛这么折腾呢? 直接计算 \(ab \bmod m\) 比较耗时(求余过程中会涉及多次除法运算,而除法运算比较耗时),后面(节 4)将介绍一种方法可以快速地计算 \(\text{Mont}(\widetilde{a},\widetilde{b})\) ,这样在 Montgomery 空间进行计算会比原空间更快。
不过,需要说明的是,如果只是进行一次模乘运算,转换到 Montgomery 空间进行计算没有优势,因为计算前需要把数字转换到 Montgomery 空间,在 Montgomery 空间得到结果后,还需要从 Montgomery 空间转换回原空间。如果是进行模幂运算(相当于多次模乘),则转移到 Montgomery 空间进行计算有很大优势,因为进入 Montgomery 空间,和退出 Montgomery 空间只需要在计算开始和结束时分别进行一次。
3. Montgomery Reduction
这节介绍 Montgomery Reduction 算法,它是下一节 4 的基础。
Montgomery Reduction 是指已知 \(m, R, T\) ,且它们满足关系 \(R > m, \gcd(R,m) = 1, 0 \le T < mR\) ,求解 \(TR^{-1} \bmod m\)。
下面以具体值 \(R = 2^5 = 64\) 来介绍一下 \(TR^{-1} \bmod m\) 的求解。
- 假设 \(T\) 能够被 \(R=2^5\) 整除,则计算 \(TR^{-1} \bmod m\) 非常快,因为这时有 \(TR^{-1} \bmod m = T >> 5\) ,也就是说把 \(T\) 右移 5 位就行了,不用进行除法运算了。由于 \(T < mR\) , \(T\) 右移 5 位后会小于 \(m\) ,不用再取模了。
- 假设 \(T\) 不能被 \(R=2^5\) 整除,那怎么办呢?思路如下,我们想办法找这样的 \(U\) ,使得 \(T + Um\) 可以被 \(R=2^5\) 整除,这样的话 \((T+Um) R^{-1} \bmod m\) 就很容易求解了,因为和前面介绍的一样,它会等于 \(T+Um >> 5\) ,右移运算不涉及除法,是很快的。而 \(TR^{-1} \bmod m\) 显然和 \((T+Um) R^{-1} \bmod m\) 是相等的,所以问题的关键就是寻找合适的 \(U\) 。前面提到对 \(U\) 的要求是 \(T + Um\) 可以被 \(R=2^5\) 整除,写成式子就是 \(T+Um \bmod R = 0\) ,从而有 \(U = -T m^{-1} \bmod R\) ,如果定义新记号 \(m' \triangleq - m^{-1} \bmod R\) (一般我们可以提前把 \(m'\) 计算出来),则 \(U = Tm' \bmod R\) 。对 \(U\) 的计算也是很快的,由于 \(R=2^5\) ,取 \(Tm'\) 的最低 5 个二进制位就行。
3.1. REDC 算法
基于上面的分析,我们可总结出 Montgomery Reduction 算法(简记为 REDC 算法,摘自:The REDC algorithm):
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
3.1.1. R 的选择
前面提到需要满足 \(R>m, \gcd(R,m) = 1\) ,下面介绍一下 \(R\) 的常用选择方法。当模数 \(m\) 确定后,假设 \(m\) 可表示为 \(n\) 位 \(b\) 进制数,即 \(m = (m_{n-1}\cdots m_1 m_0)_b\) ,则取 \(R=b^n\) 时显然会有 \(R>m\) 。而对于条件 \(\gcd(R,m) = 1\) ,比较容易满足,例如当模数 \(m\) 是奇数时,只要 \(b\) 是 2 的次方(即二进制、四进制、八进制等等),则 \(\gcd(R,m)=1\) 会成立。
假设 \(m\) 是 32 个二进制位组成(即 \(b=2\) ),则可取 \(R=2^{32}\) ;假设 \(m\) 是 5 个十进制位组成(即 \(b=10\) ),则可取 \(R=10^5\) 。在编程实现时 \(b\) 往往是 2 的次方(以方便移位运算),而在书本上介绍算法时 \(b\) 常常是 10(以方便直接验证算法的正确性)。
3.1.2. REDC 实例
设 \(m=72639, R = 10^5, T=7118368\) ,容易计算出 \(m' = -m^{-1} \bmod R = 46241\) ,下面演示一下 \(TR^{-1} \bmod m\) 的计算过程。
根据 REDC 算法,有: \[\begin{aligned} U &= (T \bmod R) m' \bmod R = (7118368 \bmod 10^5) \times 46241 \bmod 10^5 = 54688 \\ t &= (T + Um) / R = (7118368 + 54688 \times 72639) / R = 3979600000 / 10^5 = 39796 \end{aligned}\] 由于 \(t < m\) ,所以 \(TR^{-1} \bmod m\) 结果就是 \(t\) 即 \(39796\) 。
3.2. REDC(Multiprecision Version)
密码学库一般能处理上百比特位,甚至上千比特位的大数。这样的数字太大无法直接表示为计算机的单个 Machine Word,大数运算往往需要组合好几个小数运算。
节 3.1 中介绍的 REDC 算法中,需要计算 \(\bmod R\) ,如果 \(R\) 可以表示为某个数的次方(如 \(R=b^n\) );则有下面的 REDC 变种算法,不再需要 \(\bmod R\) 运算了,只需要计算 \(\bmod b\) ,如图 1 所示。
Figure 1: Montgomery Reduction
图 1 所示算法和节 3.1 中算法的核心思想是一样的,它们的区别在于,图 1 所示算法的 \(m'=-m^{-1} \bmod b\) ,而节 3.1 中则是 \(m'=-m^{-1} \bmod R\) 。这个差别导致了算法主体有些差别。后文会通过例子演示图 1 所示算法。
3.2.1. REDC 实例
设 \(m=72639, b=10, R = 10^5, T=7118368\) ,有 \(n=5\) ,容易计算出 \(m' = -m^{-1} \bmod b = 1\) ,下面演示一下使用图 1 所示算法计算 \(TR^{-1} \bmod m\) 的过程。
\(i\) | \(u_i = a_i m' \bmod 10\) | \(u_i m b^i\) | \(A = (a_{2n-1}\cdots a_3 a_2 a_1 a_0)_{10}\) |
---|---|---|---|
- | - | - | \(7118368\) |
0 | \(a_0 \times 1 \bmod 10 = 8\) | \(581112\) | \(7118368 + 581112 = 7699480\) |
1 | \(a_1 \times 1 \bmod 10 = 8\) | \(5811120\) | \(7699480 + 5811120 = 13510600\) |
2 | \(a_2 \times 1 \bmod 10 = 6\) | \(43583400\) | \(13510600 + 43583400 = 57094000\) |
3 | \(a_3 \times 1 \bmod 10 = 4\) | \(290556000\) | \(57094000 + 290556000 = 347650000\) |
4 | \(a_4 \times 1 \bmod 10 = 5\) | \(3631950000\) | \(347650000 + 3631950000 = 3979600000\) |
Step 2 结束后,得到了 \(A=3979600000\) 。再执行 Step 3,即 \(A = 3979600000/10^5 = 39796\) ;由于它小于 \(m\) ,这就是 \(TR^{-1} \bmod m\) 的结果。这个结果和节 3.1.2 中通过原始 REDC 算法的计算结果是一样的。
4. Montgomery Multiplication
在节 2 介绍中介绍了“Montgomery 空间” 的模乘运算 \(\text{Mont}(x, y) \triangleq xyR^{-1} \bmod m\) ,这节将介绍它的一个高效算法。
4.1. Mont(x,y) 算法
4.2. Mont(x,y) 实例
设 \(m=72639, b=10, R = 10^5, x=5792, y=1229\) ,有 \(n=5\) ,容易计算出 \(m' = -m^{-1} \bmod b = 1\) ,表 3 演示了使用图 2 所示算法计算 \(\text{Mont}(x, y)\) 的过程。
\(i\) | \(x_i\) | \(x_i y_0\) | \(u_i = (a_0 + x_i y_0)m' \bmod 10\) | \(x_iy\) | \(u_i m\) | \(A = (a_{n}\cdots a_3 a_2 a_1 a_0)_{10}\) |
- | - | - | - | - | - | 0 |
0 | 2 | 18 | \((0 + 18) \bmod 10 = 8\) | 2458 | 581112 | 58357 |
1 | 9 | 81 | \((7 + 81) \bmod 10 = 8\) | 11061 | 581112 | 65053 |
2 | 7 | 63 | \((3 + 63) \bmod 10 = 6\) | 8603 | 435834 | 50949 |
3 | 5 | 45 | \((9 + 45) \bmod 10 = 4\) | 6145 | 290556 | 34765 |
4 | 0 | 0 | \((5 + 0) \bmod 10 = 5\) | 0 | 363195 | 39796 |
由于得到的 \(A=39796\) 小于 \(m\) ,所以 \(\text{Mont}(x, y)\) 的最终结果就是 \(39796\) 。
5. 使用 Mont(x,y) 快速计算模幂
节 2 中提到过,如果只是进行一次模乘运算,转换到 Montgomery 空间进行计算没有优势,因为计算前需要把数字转换到 Montgomery 空间,在 Montgomery 空间得到结果后,还需要从 Montgomery 空间转换回原空间。 Montgomery 算法真正的应用场景是计算模幂,即 \(x^e \bmod m\) 。 因为进入 Montgomery 空间,和退出 Montgomery 空间只需要在计算开始和结束时分别进行一次。
在介绍使用 Mount(x,y) 快速计算模幂前,先介绍一下“Left-to-right binary exponentiation”算法,如图 3 所示。
Figure 3: Left-to-right binary exponentiation
下面以 \(g^{283}\) 为例,介绍一下 Left-to-right binary exponentiation 算法的应用过程。
先把十进制 \(283\) 表示为二进制,即 \((100011011)_2\) ,从而 \(t=8\) ,按照算法的 Step 2 进行迭代,如表 4 所示。
\(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 算法简单。
5.1. Montgomery Exponentiation 算法
结合节 4.1 中求解 \(\text{Mont}(x,y)\) 的算法和前面介绍的 Left-to-right binary exponentiation 算法,可以得到图 4 所示的 Montgomery Exponentiation 算法。
Figure 4: Montgomery Exponentiation
其中算法 Step 1 中的 \(\widetilde{x} \leftarrow \text{Mont}(x, R^2 \bmod m)\) 就是把 \(x\) 转换到 Montgomery 空间。因为 \(\text{Mont}(x, R^2 \bmod m) = x R^2 R^{-1} \bmod m = x R \bmod m\) ,这就是数 \(x\) 在 Montgomery 空间的定义。
5.1.1. 求模幂实例
设 \(x, m, R\) 满足图 4 算法中指定的要求,设 \(e=11=(1011)_2\) ,从而 \(t=3\) ,下面演示一下利用图 4 算法求解 \(x^e \bmod m\) 的过程。
对这个例子应用算法 4 中 Step 2,其过程如表 5 所示。
\(i\) | \(e_i\) | \(A \bmod m\) |
---|---|---|
3 | 1 | \(\widetilde{x}\) |
2 | 0 | \(\widetilde{x}^2R^{-1}\) |
1 | 1 | \(\widetilde{x}^5R^{-4}\) |
0 | 1 | \(\widetilde{x}^{11}R^{-10}\) |
Step 2 结束后,得到的 \(A=\widetilde{x}^{11}R^{-10}\) ,执行 Step 3,有 \(A = \widetilde{x}^{11}R^{-10}R^{-1} = (xR^{11})R^{-10}R^{-1} = x^{11} \bmod m\) ,显然这就是待求的 \(x^e \bmod m\) 。
6. Montgomery 的不同实现
前面介绍的算法都是教科书中算法,在编程实现时,往往还有其它一些优化,论文 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
7. 参考
- 本文主要总结于 Handbook of Applied Cryptography,Chapter 14 Efficient Implementation
- https://en.wikipedia.org/wiki/Montgomery_modular_multiplication
- 1985 年 Peter L. Montgomery 发表的论文 Modular Multiplication Without Trial Division
- 论文 Comparison of three modular reduction functions 中对传统方法、Barrett Reduction 算法和 Montgomery 算法进行了性能比较
- Efficient Software Implementations of Modular Exponentiation
- Rethinking Modular Multi-Exponentiation in Real-World Applications