蒙哥马利算法

蒙哥马利模乘运算（Montgomery Modular Multiplication）[1]与蒙哥马利幂模运算（Montgomery power module）和蒙哥马利约减运算（Montgomery model reduction）统称蒙哥马利算法（Montgomery Algorithm）。

• 蒙哥马利模乘运算：$abmodN$
• 蒙哥马利幂模运算：$a^{b}modN$
• 蒙哥马利约减运算：$a{b}^{-1}modN$

蒙哥马利模乘#

模乘是计算$ab(\mathrm{mod}N)$。在普通计算模N时，利用的是带余除法，除法运算需要太多次乘法，计算复杂度较高，「蒙哥马利模乘」的思想就是利用进制表示简化除法运算，转化成位运算。

• 蒙哥马利形式(Montgomery form, Montgomery domain)

为了计算$ab(\mathrm{mod}N)$，需要找到一个$R$，然后使得${\mathrm{a}}^{\prime }\equiv \mathrm{aR}(\mathrm{mod}\mathrm{N}),{\mathrm{b}}^{\prime }\equiv \mathrm{bR}(\mathrm{mod}\mathrm{N})$。

这个$R$不是随便取的，需要满足两个条件：

1. $R=2^k>N$ ，其中$k$是满足条件的最小正数，这就能保证除以$R$就相当于右移$k$位
2. $gcd(R,N)=1$，这就可以求出$m$

在计算$ab(\mathrm{mod}N)$时需要利用「蒙哥马利形式」。令$X=a^{\prime }{b}^{\prime }$，可以设计一个函数$REDC(X)=X{R}^{-1}(\mathrm{mod}N)$，计算结果为$X_1=REDC(X)=a^{\prime }{b}^{\prime }{R}^{-1}(\mathrm{mod}N)=abR(\mathrm{mod}N)$，这样再调用一次函数计算$REDC(X_1)=X_1{R}^{-1}(\mathrm{mod}N)$就能得到结果$ab(\mathrm{mod}N)$。

这个函数$REDC()$就是蒙哥马利约减算法，即求$X{R}^{-1}(\mathrm{mod}N)$。

所以蒙哥马利模乘可以分三步进行计算：

1. 将输入$a{R}^2,b{R}^2$转成蒙哥马利形式，即$aR=REDC(a{R}^2),bR=REDC(b{R}^2)$；
2. 做一次标准乘法得$ab{R}^2=aR\ast bR$；
3. 最后做一次REDC得$abR=REDC(ab{R}^2)$；
4. 注意上面三个步骤返回的是蒙哥马利形式的$ab$，即$abR$。若需要转换成正常形式的$ab$，需要再做一次$REDC$得$ab=REDC(abR)$；

代码实现#

给出一个蒙哥马利模乘的Python实现计算 (23456789*12345678)%123456789，注意程序里面取$R=2^{64}$， 所以$modR$相当于取低64bit，$/R$相当于右移64bit。

import math

class MontMul:
    """docstring for ClassName"""
    def __init__(self, R, N):
        self.N = N
        self.R = R
        self.logR = int(math.log(R, 2)) #log_2 ^ R
        N_inv = MontMul.modinv(N, R)
        self.N_inv_neg = R - N_inv    #N_inv_neg=R-N^{-1}
        self.R2 = (R*R)%N

    @staticmethod        
    def egcd(a, b):
        if a == 0:
            return (b, 0, 1)
        else:
            g, y, x = MontMul.egcd(b % a, a)
            return (g, x - (b // a) * y, y) #g=gcd(a,b)

    @staticmethod
    def modinv(a, m):
        g, x, y = MontMul.egcd(a, m)
        if g != 1:
            raise Exception('modular inverse does not exist')
        else:
            return x % m

    def REDC(self, T):
        N, R, logR, N_inv_neg = self.N, self.R, self.logR, self.N_inv_neg

        m = ((T&int('1'*logR, 2)) * N_inv_neg)&int('1'*logR, 2) # m = (T%R * N_inv_neg)%R        
        t = (T+m*N) >> logR # t = int((T+m*N)/R)
        if t >= N:
            return t-N
        else:
            return t

    def ModMul(self, a, b):
        if a >= self.N or b >= self.N:
            raise Exception('input integer must be smaller than the modulus N')

        R2 = self.R2
        aR = self.REDC(a*R2) # convert a to Montgomery form
        bR = self.REDC(b*R2) # convert b to Montgomery form
        T = aR*bR # standard multiplication
        abR = self.REDC(T) # Montgomery reduction
        return self.REDC(abR) # covnert abR to normal ab

if __name__ == '__main__':
    N = 123456789
    R = 2**64 # assume here we are working on 64-bit integer multiplication
    g, x, y = MontMul.egcd(N,R)
    if R<=N or g !=1: 
        raise Exception('N must be larger than R and gcd(N,R) == 1')
    inst = MontMul(R, N)

    input1, input2 = 23456789, 12345678
    mul = inst.ModMul(input1, input2)
    if mul == (input1*input2)%N:
        print ('({input1}*{input2})%{N} is {mul}'.format(input1 = input1, input2 = input2, N = N, mul = mul))

蒙哥马利模约简#

蒙哥马利模约简（REDC）是蒙哥马利模乘最重要的部分，主要是计算 $T{R}^{-1}(\mathrm{mod}N)\leftarrow REDC(T)$，算法描述如下：

一般做模约减运算$T{R}^{-1}modN$，相当于$\frac{T}{R}(\mathrm{mod}N)$，需要做一次除运算，如何避免除法呢？

由于$R=2^k$，所以$\frac{T}{R}$相当于T右移k位，即$T>>k$。但右移k位可能会抹掉T低位中的一些1，如$7÷4=0b111>>2=0b1=1$，这个不是精确计算，而是向下取整的除法，当且仅当T是R的整数倍时，$T/R=T>>k$。所以实际上就是找一个$m$，使得$T+mN$是$R$的倍数，这样计算$\frac{T+mN}{R}$就相当于$(T+mN)>>k$ 。

• $m$如何找？

由于$gcd(R,N)=1$，根据扩展欧几里得算法得：有$R{R}^{\prime }-N{N}^{\prime }=1$，且有$1<N^{\prime }<R,0<R^{\prime }<N<R$。（这里是$-N{N}^{\prime }$，所以$N^{\prime }=-N^{-1}\mathrm{mod}R$）

扩展欧几里得:

若$gcd(a,b)=1$，则必存在整数$x,y$，使得$ax+by=gcd(a,b)=1$。

$\begin{array}{c}\mathrm{T}+\mathrm{mN}\equiv 0(\mathrm{mod}R)\\ {\mathrm{TN}}^{\prime }+{\mathrm{mNN}}^{\prime }\equiv 0(\mathrm{mod}R)\\ {\mathrm{TN}}^{\prime }+\mathrm{m}({\mathrm{RR}}^{\prime }-1)\equiv 0(\mathrm{mod}R)\\ {\mathrm{TN}}^{\prime }\equiv \mathrm{m}(\mathrm{mod}\mathrm{R})\end{array}$

这样就求出了$m=T{N}^{\prime }(\mathrm{mod}R)$。

所以在已知$a,b,N$，并计算出$a^{\prime },b^{\prime },R,T$下，蒙哥马利约减算法$REDC(T)=T{R}^{-1}(\mathrm{mod}N)$如下：

1. 计算$N^{\prime }=-N^{-1}(\mathrm{mod}R),m=T{N}^{\prime }(\mathrm{mod}R)$；
2. 计算$y=\frac{T+mN}{R}$，即将$T+mN$右移$k$位；
3. 若$y>N$，则$y=y-N$，由于$\mathrm{X}<{\mathrm{N}}^2,\mathrm{m}<\mathrm{R},\mathrm{N}<\mathrm{R}$，所以$0<y<2N$，故：$\frac{\mathrm{X}+\mathrm{m}\mathrm{N}}{\mathrm{R}}<\frac{{\mathrm{N}}^2+\mathrm{R}\cdot \mathrm{N}}{\mathrm{R}}<\frac{\mathrm{R}\mathrm{N}+\mathrm{R}\cdot \mathrm{N}}{\mathrm{R}}=2\mathrm{N}$
4. 返回$y=T{R}^{-1}(\mathrm{mod}N)$；

蒙哥马利幂模运算#

蒙哥马利幂模运算是快速计算$a^{b}modN$的一种算法，是RSA加密算法的核心之一。

蒙哥马利幂模的优点在于减少了取模的次数（在大数的条件下）以及简化了除法的复杂度（在2的k次幂的进制下除法仅需要进行左移操作）。

计算方式是将幂模运算转化为模乘运算

例如：求D=C^15 % N

由于：ab % n = (a % n)(b % n) % n

所以令：

C = C%N

C1 =C*C % N =C^2 % N

C2 =C1*C % N =C^3 % N

C3 =C2*C2 % N =C^6 % N

C4 =C3*C % N =C^7 % N

C5 =C4*C4 % N =C^14 % N

C6 =C5*C % N =C^15 % N

即：对于E=15的幂模运算可分解为6 个乘模运算。

归纳分析以上方法可以发现：对于任意指数E，都可采用以下算法计算 D=C^E % N：

D=1

WHILE E>0

  IF E%2=0

    C=C*C % N

    E=E/2

  ELSE

    D=D*C % N

    E=E-1

RETURN D

继续分析会发现，要知道E何时能整除 2，并不需要反复进行减一或除二的操作，只需验证E 的二进制各位是0还是1就可以了，从左至右或从右至左验证都可以，从左至右会更简洁。

D=1

FOR i=n TO 0

  D=D*D % N

  IF e=1	//e是E的最后一位【判断E是否为奇数】

    D=D*C % N
  
RETURN D

这样，模幂运算就转化成了一系列的模乘运算。

/*例如求D=C^15%N 
由于：C*k % n = (C % n)*(k % n) % n 
所以令: 
【奇数】 C1 = C*C % N =C^2 % N       1   15 
			  C2 = C1*C % N =C^3 % N      3   7 
【奇数】C3 = C2*C2 % N =C^6 % N 
				C4 = C3*C % N =C^7 % N      7   3 
【奇数】C5 = C4*C4 % N =C^14 % N 
				C6 = C5*C % N =C^15 % N     15  1 
				
蒙哥马利幂模运算*/   
#include <iostream> 
using namespace std; 
typedef long long  __int64;
  
__int64 Montgomery(__int64 base,__int64 exp,__int64 mod)  
{  
    __int64 res = 1;  
    while(exp)  
    {  
        if ( exp&1 )  //取exp的最后一位为1（奇数）
            res = (res*base) % mod;  
        exp >>= 1;    //exp/2
        base = (base*base) % mod;  
    }  
    return res;  
}  
int main()  
{  
    //base 底数，exponential　指数，mod 模  
    __int64 base,exp,mod;                   //base^exp % mod
    base=12,exp=15,mod=99;
    cout << base<<"^"<<exp<<"%"<<mod<<"="<<Montgomery(base,exp,mod) << endl;  
    return 0;  
}  

复杂度分析#

• 将模除运算转换为移位运算；

• 当出现大量模乘运算时，可以通过并行运算进行预计算，节省时间；

参考#

1. https://en.wikipedia.org/wiki/Montgomery_modular_multiplication
2. 密码学中的模乘算法——蒙哥马利模乘（Montgomery Modular Multiplication）
3. 蒙哥马利约减算法
4. 蒙哥马利模乘
5. 蒙哥马利算法(Montgomery Algorithm)|蒙哥马利约简、模乘、模幂