用途：一种高速的取模算法，当$a$，$b$，$N$很大时，加速计算$a*b(modN)$。

一、原理

1. 将$a$，$b$映射到Montgomery域内：

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

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

2. $Montgomery$约简：

$$
\begin{align}
&Input:X \
&Output:X/R
\end{align}
$$

约简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)
$$
这样，约简函数的输入输出就变为：
$$
\begin{align}
&Input:X \
&Output:(X+mN)/R
\end{align}
$$

由拓展的欧几里得算法可以解出以下方程的一组特解：
$$
RR’-NN’=gcd(R,N)
$$
对上述两个公式进行推导，就可以得到$m$的计算方法：
$$
\begin{align}
&XN’+mNN’=0(modR)\
&XN’+m(RR’-gcd(R,N))=0(modR)\
&XN’=m*gcd(R,N)(modR)\
&m=XN’/gcd(R,N)(modR)
\end{align}
$$

通过阅读蒙哥马利算法原文可知，该算法限制$gcd(R,N)=1$，
所以接着推演可以得到
$$
\begin{align}
&X+mN=X+XNN’(modR)=X+(XN’(modR))N\
&(X+mN)/R=(X+(XN’(modR))N)/R
\end{align}
$$

二、优化加速

1. 分析

1次蒙哥马利算法中，进行了2次对$N$取模（$a->a’，b->b’$），2次对$R$取模（$reduction$），这就是我们需要优化的地方。对于Montgomery算法本身，它是把除法转化为移位算法，以下运算得到了简化：
$$
\begin{align}
&aR=a<<k\
&bR=b<<k\
&(X+mN)/R=(X+mN)>>k
\end{align}
$$

2. 对$N$取模

第1步计算$a’=aR(modN),b’=bR(modN)$时，若按照原来计算方法，一条转换的运算量可能就超过了$a*b(modN)$，所以通过模运算法则简化为：
$$
\begin{align}
&a’=(a(modN)*R(modN))(modN)\
&b’=(b(modN)*R(modN))(modN)
\end{align}
$$

3. 对$R$取模

在reduction函数里面，进行了对R取模的运算，具体为$XN’(modR)$。因为$R=2^k$，所以$XN’(modR)$相当于取$XN’$的后$k-1$位。
$$
XN’MODR = XN’[0:k-1]
$$
此处，$[0:k-1]$为XN’的第$0~(k-1)$位。

三、Python实现

声明：下面的代码来自GitHub，原作者是LivingLeopold

import math

def 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, y

def 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-N1

a_ =( a * R ) % N
b_ =( b * R ) % N
c_ =((( a_ * b_ ) + ((( a_ * b_ ) * N2 ) % R ) * N )/ R ) % N

c = 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代码实现蒙哥马利算法
[8] Markdown中输入多行并列的公式