最新要闻

广告

手机

iphone11大小尺寸是多少?苹果iPhone11和iPhone13的区别是什么?

iphone11大小尺寸是多少?苹果iPhone11和iPhone13的区别是什么?

警方通报辅警执法直播中被撞飞:犯罪嫌疑人已投案

警方通报辅警执法直播中被撞飞:犯罪嫌疑人已投案

家电

世界最资讯丨蒙哥马利算法

来源:博客园

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


(资料图片仅供参考)

  • 蒙哥马利模乘运算:\(ab mod N\)
  • 蒙哥马利幂模运算:\(a^b mod N\)
  • 蒙哥马利约减运算:\(ab^{-1} mod N\)

蒙哥马利模乘

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

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

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

这个\(R\)不是随便取的,需要满足两个条件:

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

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

这个函数\(REDC()\)就是蒙哥马利约减算法,即求\(XR^{-1}\pmod{N}\)。

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

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

代码实现

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

import mathclass 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 abif __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)是蒙哥马利模乘最重要的部分,主要是计算 \(TR^{-1}\pmod{N} \gets REDC(T)\),算法描述如下:

一般做模约减运算\(TR^{-1} mod N\),相当于$\frac{T}{R}\pmod{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 + m N\)是\(R\)的倍数,这样计算\(\frac{T+mN}{R}\)就相当于\((T+mN) >>k\) 。

  • \(m\)如何找?

由于\(gcd(R,N)=1\),根据扩展欧几里得算法得:有\(RR"-NN"=1\),且有\(1

扩展欧几里得:

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

\(\begin{gathered}\mathrm{T+mN}\equiv0{\pmod{R}} \\\mathrm{TN"+mNN"\equiv0}\pmod{R} \\\mathrm{TN}^{\prime}+\mathrm{m}(\mathrm{RR}^{\prime}-1)\equiv0{\pmod{R}} \\\mathrm{TN}^{\prime}\equiv\mathrm{m}\pmod{\mathrm{R}} \end{gathered}\)

这样就求出了\(m=TN"\pmod{R}\)。

所以在已知\(a,b,N\),并计算出\(a",b",R,T\)下,蒙哥马利约减算法\(REDC(T)=TR^{-1}\pmod{N}\)如下:

  1. 计算\(N"=-N^{-1}\pmod{R},m=TN"\pmod{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
  4. 返回\(y=TR^{-1}\pmod{N}\);

蒙哥马利幂模运算

蒙哥马利幂模运算是快速计算\(a^b mod N\)的一种算法,是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=1WHILE E>0  IF E%2=0    C=C*C % N    E=E/2  ELSE    D=D*C % N    E=E-1RETURN D

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

D=1FOR 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  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<<"^"<

复杂度分析

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

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

参考

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

关键词: