Montgomery算法原理及Python实现

论文原文：Modular Multiplication Without Trial Division

用处：当a，b，N很大时，减速计算a*b(modN)

一、步骤

将a，b映射到Montgomery域内：

$$a'=aR(modN),b'=bR(modN)$$

此处的$R=2^k>N$，令$X=a'b'=abR^2(modN)$

Montgomery约简：

$$Input:X\\Output:X/R$$

约简1次：$X=abR^2(modN)$ -> $X/R=abR(modN)$

约简2次：$X=abR(modN)$ -> $X/R=ab(modN)$

X/R是不能够间接除的，因为X的低k位很可能不都为0。因而心愿有一个数m使得

$$X+mN=0(modR)$$

由拓展的欧几里得算法能够解出以下方程的一组特解：

$$RR'-NN'=gcd(R,N)$$

对上述两个公式进行推导，就能够失去m的计算方法：

$$XN'+mNN'=0(modR)\\XN'+m(RR'-gcd(R,N))=0(modR)\\XN'=m*gcd(R,N)(modR)\\m=XN'/gcd(R,N)(modR)$$

蒙哥马利算法限度$gcd(R,N)=1$，接着

$$X+mX=X+XNN'(modR)=X+(XN'(modR))N\\(X+mX)/R=(X+(XN'(modR))N)/R$$

二、减速的起因：变乘除为移位

以下运算将被大大简化：

$$aR=a<<k\\bR=b<<k\\(X+mN)/R=(X+mN)>>k$$

1次蒙哥马利算法中，进行了2次对N取模，2次对R取模。

三、Python实现

申明：上面的代码来自GitHub，原作者是LivingLeopold

import mathdef get_gcd(a,b):    if a%b == 0:        return b    else :        return get_gcd(b,a%b)def get_(a, b):    if b == 0:        return 1, 0    else:        k = a // b        remainder = a % b        x1, y1 = get_(b, remainder)        x, y = y1, x1 - k * y1    return x, ydef Multiplicative_inverse_modulo( a , b ):    if b < 0:        m = abs(b)    else:        m = b    flag = get_gcd(a, b)    if flag == 1:        x, y = get_(a, b)        x0 = x % m        return x0    else:        print("Error!")        exit()print('To calculate a * b ( mod N ),')a , b , N = eval(input('Please input a , b , N (Like 68,57,109)：'))R = 2**(int(math.log2(N)+1))N1 = Multiplicative_inverse_modulo( N , R )N2 = R-N1a_ =( a * R ) % Nb_ =( b * R ) % Nc_ =((( a_ * b_ ) + ((( a_ * b_ ) * N2 ) % R ) * N )/ R ) % Nc = int(((c_ + ( c_* N2 ) % R * N ))/R)print( str(a) + ' * ' + str(b) + ' ( mod ' + str(N) +' ) = ' + str(c))

参考文献

[1] 高效幂模算法探索：Montgomery算法解析
[2] 蒙哥马利算法详解
[3] 蒙哥马利算法约简、模乘、幂模
[4] 欧几里得算法&辗转相除法，百度百科
[5] 拓展的欧几里得算法
[6] RSA大数运算实现（1024位n）（5）蒙哥马利模幂
[7] Python代码实现蒙哥马利算法