Montgomery reduction——多精度模乘法运算算法

1. modular reduction模简化定义

modular reduction模简化的定义为：
若z为任意整数，则 z mod m 的结果在区间[0, m-1]，相当于z除以m的余数，该运算过程称为z对模m的模简化。

2. 模运算

有限域Z_m内的加减乘除运算，其中的m为多精度正整数，m称为模。

将正整数m，非负整数x,y (x<m y<m)以b进制格式表示：
m=(m_nm_n-1…m₁m₀)_b
x=(x_nx_n-1…x₁x₀)_b
y=(y_ny_n-1…y₁y₀)_b

模运算：

• 模加法：x + y mod m
• 模减法：x - y mod m
• 模乘法：x * y mod m
• 模取反：x^-1 mod m

2.1 模加法和模减法

对于常规的多精度运算来说，模加法和模减法为最简单的运算。

举例：
当0<x<m, 0<y<m时，则有：
1）x+y < 2m;
2) 若x >= y，则0 =< x-y < m;
3) 若x <y，则有0 =< x+m-y < m.

因此对于模加法和模减法，可分别采用参考资料1中的$14.7算法（当x+y>=m时，需额外增加一步-m的操作）和$14.9算法（当x >= y）。

2.2 乘法运算

以b进制表示：
x=(x_nx_n-1…x₁x₀)_b 【n+1位】
y=(y_ty_t-1…y₁y₀)_b 【t+1位】
则 x*y 的乘积结果以b进制表示，最多有(n+t+2)位。

小学笔算级别求乘积的算法如下：

举例：

2.3 除法运算

除法运算为基础多精度加减乘除运算中最复杂且最耗时的运算。
以下算法为计算x除以y的商q和余数r，即：
x=y*q+r, 0=<r<y

2.4 模乘法

模乘法既需要2.2中的乘法运算，也需要第一节中提到的模简化运算。

2.4.1 模乘法classical算法

传统直观的求解模乘法x*y mod m的算法可分解为：
1）用2.2乘法运算求得xy；
2）用2.3除法运算球的xy除以m的余数r;
3）返回r值。

该算法简单直接，但效率不高。

2.4.2 模乘法Montgomery reduction算法

• Montgomery reduction定义： 模m为正整数，R，T为整数，且R>m，gcd(m,R) = 1, 0 =< T < mR，则 TR^-1 mod m 称为相对于R的模m的T Montgomery reduction运算。

通过选择合适的R值，Montgomery reduction用于模乘法运算效率将更高。
对于0 =< T < mR，T可取值T=aR，其中0 =< a < m。

对于 0 =< x,y < m，假设
X = xR mod m
Y = yR mod m
则XY的Montgomery reduction为 XYR^-1 mod m = xyR mod m。

对于任意的1 =< x < m，若需求 x⁵ mod m 的值，可转换为：
1）求X = xR mod m；
2）求X²的Montgomery reduction A = X²R^-1 mod m；
3）求A² 的Montgomery reduction B = A²R^-1 mod m = X⁴ R^-3 mod m；
4）求BX的Montgomery reduction C = BXR^-1 mod m = X⁴R^-3XR^-1 mod m = X⁵R^-4 mod m = x⁵R mod m;
5) C*(R^-1 mod m) mod m = x⁵ mod m.

通过以上五步可看出，通过Montgomery reduction进行模运算比2.4.1的classical算法效率更高。

实际操作时，一般地，若模m值以b进制表示长度为n，则通常地，R的取值为bⁿ，由此 R > m条件肯定满足。但是gcd(R,m)=1 条件成立的前提需为gcd(b,m)=1。对于素数有限域内的运算，m为素数，所以gcd(b,m)=1 亦成立。

Montgomery的一些特征：

• Fact(Montgomery reduction)：若m,R为整数，且gcd(m,R)=1。设定 m’ = -m^-1 mod R，T为整数且 0 =< T < mR。若U = Tm’ mod R，则 (T+Um)/R 为整数，且 (T+Um)/R=TR^-1 mod m.

以上Fact的证明如下：
∵ T+Um = T mod m
∴ (T+Um)R^-1 = TR^-1 mod m
∵ m’ = -m^-1 mod R
∴ mm’ = -1 + nR 其中n为整数
∵ U = Tm’ mod R
∴ U = Tm’ + kR 其中k为整数
∴ (T+Um)/R = (T+(Tm’+kR)m)/R = (T+Tmm’+mkR)/R = (T+T(-1+nR)+mkR)/R = nT+km
∴ (T+Um)/R为整数。

具体的算法细节为：

参考资料：
[1] The Montgomery reduction here is based on Algorithm 14.32 in Handbook of Applied Cryptography chap14 http://cacr.uwaterloo.ca/hac/about/chap14.pdf.