蒙哥马利算法

Posted 2023-05-11 Pam

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

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

蒙哥马利模乘

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

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

为了计算\\(ab\\pmodN\\)，需要找到一个\\(R\\)，然后使得\\(\\mathrma^\\prime\\equiv\\mathrmaR\\pmod\\mathrmN,\\mathrmb^\\prime\\equiv\\mathrmbR\\pmod\\mathrmN\\)。

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

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

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

这个函数\\(REDC()\\)就是蒙哥马利约减算法，即求\\(XR^-1\\pmodN\\)。

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

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 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）是蒙哥马利模乘最重要的部分，主要是计算 \\(TR^-1\\pmodN \\gets REDC(T)\\)，算法描述如下：

一般做模约减运算\\(TR^-1 mod N\\)，相当于$\\fracTR\\pmodN $，需要做一次除运算，如何避免除法呢？

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

• \\(m\\)如何找？

由于\\(gcd(R,N)=1\\)，根据扩展欧几里得算法得：有\\(RR\'-NN\'=1\\)，且有\\(1<N\'<R,0<R\'<N<R\\)。（这里是\\(-NN\'\\)，所以\\(N\'=-N^-1\\modR\\)）

扩展欧几里得:

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

\\(\\begingathered \\mathrmT+mN\\equiv0\\pmodR \\\\ \\mathrmTN\'+mNN\'\\equiv0\\pmodR \\\\ \\mathrmTN^\\prime+\\mathrmm(\\mathrmRR^\\prime-1)\\equiv0\\pmodR \\\\ \\mathrmTN^\\prime\\equiv\\mathrmm\\pmod\\mathrmR \\endgathered\\)

这样就求出了\\(m=TN\'\\pmodR\\)。

所以在已知\\(a,b,N\\)，并计算出\\(a\',b\',R,T\\)下，蒙哥马利约减算法\\(REDC(T)=TR^-1\\pmodN\\)如下：

1. 计算\\(N\'=-N^-1\\pmodR,m=TN\'\\pmodR\\)；
2. 计算\\(y=\\fracT+mNR\\)，即将\\(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=TR^-1\\pmodN\\)；

蒙哥马利幂模运算

蒙哥马利幂模运算是快速计算\\(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=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)|蒙哥马利约简、模乘、模幂

蒙哥马利算法详解

这篇文章为大家梳理一下整个蒙哥马利算法的本质，蒙哥马利算法并不是一个独立的算法，而是三个相互独立又相互联系的算法集合，其中包括

• 蒙哥马利乘模，是用来计算$x\cdot y\ (mod\ N)$
• 蒙哥马利约减，是用来计算$t\cdot \rho^{-1}\ (mod\ N)$
• 蒙哥马利幂模，是用来计算$x^y\ (mod\ N)$

其中蒙哥马利幂乘是RSA加密算法的核心部分。

基本概念

梳理几个概念，试想一个集合是整数模N之后得到的
$Z_N=\left\{0,1,2,\cdots,N-1\right\}$

注：N在base-b进制下有$l_N$位。 比如10进制和100进制，都属于base-10进制，因为$100=10^2$，所以b=10。在10进制下，667的$l_N=3$

这样的集合叫做N的剩余类环，任何属于这个集合Z的x满足以下两个条件：
1. 正整数
2. 最大长度是$l_N$

这篇文章中讲到的蒙哥马利算法就是用来计算基于$Z_N$集合上的运算，简单讲一下原因，因为RSA是基于大数运算的，通常是1024bit或2018bit，而我们的计算机不可能存储完整的大数，因为占空间太大，而且也没必要。因此，这种基于大数运算的加密体系在计算的时候都是基于$Z_N$集合的，自然，蒙哥马利算法也是基于$Z_N$。

在剩余类环上，有两种重要的运算，一类是简单运算，也就是加法和减法，另一类复杂运算，也就是乘法。我们比较熟悉的是自然数集上的运算，下面看下怎么从自然数集的运算演变成剩余类环上的运算。

对于加法运算，如果计算$x\pm y\ (mod\ N)$ ($0\leqslant x,y<N$)，试想自然数集上的 $x\pm y$

$\qquad 0\leqslant x+y\leqslant 2\cdot(N-1)$

$-(N-1)\leqslant x-y\leqslant (N-1)$

我们可以简单的通过加减N来实现从自然数到剩余类集的转换

另外一类是乘法操作，也就是$x\cdot y\ (mod\ N)$($0\leqslant x,y<N$)，那么

$0\leqslant x\cdot y\leqslant (N-1)^2$

如果在自然数集下，令$t=x\cdot y$，那么对于$\mod N$我们需要计算

$t-（N\cdot \lfloor\frac{t}{N}\rfloor）$

加减操作很简单，具体的算这里就不细说了，我们用$Z_N-ADD$ 来代表剩余类环上的加法操作。既然我们可以做加法操作，那么我们就可以扩展到乘法操作，算法如下

但是这并不是一个好的解决方案，因为通常来说，我们不会直接做w位乘w位的操作，这个后面会用蒙哥马利的乘法来代替解决。

对于取模操作，一般有以下几种方法

1，根据以下公式，来计算取模操作

$t-（N\cdot \lfloor\frac{t}{N}\rfloor）$

这种解法有以下特征

• 整个计算过程是基于标准的数字表示
• 不需要预计算（也就是提前计算一些变量，以备使用）
• 涉及到一个除法操作，非常费时和复杂

2，用Barrett reduction算法，这篇文章不细说，但是有以下特征

• 基于标准的数字表示
• 不需要预计算
• 需要$2 \cdot (l_N+1) \cdot (l_N+1)$ 次数乘运算

3，用蒙哥马利约减，也就是下面要讲的算法，有以下特征

• 不是基于标准的数字表示（后文中有提到，是基于蒙哥马利表示法）
• 需要预计算
• 需要$2 \cdot (l_N) \cdot (l_N)$ 次数乘运算

蒙哥马利预备知识

在将蒙哥马利算法之前，先看一下在自然数下的乘法公式

计算$x\cdot y$，想象一下我们常用的计算乘法的方法，用乘数的每一位乘上被乘数，然后把得到的结果相加，总结成公式，可以写成如下的形式。

$x\cdot y=x\cdot sum_{i=0}^{l_y-1}y_i \cdot b^i$

$\qquad=sum_{i=0}^{l_y-1}y_i \cdot x \cdot b^i$

尝试下面一个例子，10进制下（也就是b=10），y=456（也就是$l_n=3$），计算$x\cdot y$，公式可演变如下：

$x\cdot y=(y_{0}\cdot x\cdot 10^{0})+(y_{1}\cdot x\cdot 10^{1})+(y_{2}\cdot x\cdot 10^{2})$
$\qquad=(y_{0}\cdot x\cdot 0)+(y_{1}\cdot x\cdot 10)+(y_{2}\cdot x\cdot 100)$
$\qquad=(y_{0}\cdot x)+10\cdot(y_{1}\cdot x+10\cdot(y_{2}\cdot x\cdot +10\cdot(0)))$

最后一次演变其实就是用霍纳法则(Horner’s rule)所讲的规则，关于霍纳法则，可以自行百度。

这个计算过程写成代码实现的算法应该是这样的：

接下来我们来看下面这样的计算，计算$(x\cdot y)/1000$，由前面可以知道

$x\cdot y=(y\_{0}\cdot x)+10\cdot(y\_{1}\cdot x+10\cdot(y\_{2}\cdot x\cdot +10\cdot(0)))$

由此可以知道：

$\frac{x\cdot y}{1000}=\frac{(y_{0}\cdot x\cdot 10^{0})+(y_{1}\cdot x\cdot 10^{1})+(y_{2}\cdot x\cdot 10^{2})}{1000}$

$\qquad=\frac{(y_{0}\cdot x\cdot 0)+(y_{1}\cdot x\cdot 10)+(y_{2}\cdot x\cdot 100)}{1000}$

$\qquad=\frac{(y_{0}\cdot x)}{1000}+\frac{(y_{1}\cdot x)}{100}+\frac{(y_{2}\cdot x)}{10}$

$\qquad=(((((y_0\cdot x)/10)+y_1\cdot x)/10)+y_2\cdot x)/10$

这个计算过程写成代码实现的算法是这样的：

接下来我们再来看在剩余类集合下的乘法操作 $x\cdot y/1000\ (mod\ 667)$

我们知道剩余类集合$Z_{667}=\left\{0,1 \cdots 666\right\}$，是不存在小数的，而如果我们采用自然数集的计算方式的话，就会出现小数，比如前面的例子，除10就会有小数。

这个问题是这样的，我们知道$u·667 \equiv 0 (mod 667)$（$\equiv$表示取模相等），所以我们可以选择一个合适的u，用u乘667再加上r，使得和是一个可以除10没有小数，这样在mod 667之后依然是正确的结果。至于u怎么算出来，这篇文章会在后面的章节说明。

这个过程之后$x\cdot y/1000\ (mod\ 667)$ 的计算算法可以写成如下的形式

至此，你可能还不明白上面说这一堆演变的原因，其实很简单，原来是一个$(x\cdot y)\ (mod\ 667)$的运算，这个运算中的模操作，正常情况下是要通过除法实现的，而除法是一个特别复杂的运算，要涉及到很多乘法，所以在大数运算时，我们要尽量避免除法的出现。而通过以上几个步骤，我们发现$(x\cdot y)/1000\ (mod\ 667)$这个操作是不用除法的。等等，算法中明明有个除10的操作，你骗谁呢。不知道你有没有发现，除数其实是我们的进制数，除进制数在计算机中是怎么做呢，其实很简单，左移操作就ok了。所以这个计算方法是不涉及到除法操作的。

但是我们要计算的明明是$(x_1\cdot y_1)\ (mod\ 667)$，怎么现在变成了$(x_2\cdot y_2)/1000\ (mod\ 667)$，所以在下一步，我们要思考的是怎么样让$(x_1\cdot y_1)\ (mod\ 667)$转变成$(x_2\cdot y_2)/1000\ (mod\ 667)$这种形式。

考虑这样两个算法
- 第一个是输入x和y，计算$x \cdot y \cdot \rho^{-1}$
- 第二个算法，输入一个t，计算$t \cdot \rho^{-1}$。

$x\cdot y\ (mod\ 667)=((x\cdot1000)\cdot(y\cdot1000)/1000)/1000\ (mod\ 667)$

是不是变成了我们需要的$(x\cdot y)/1000\ (mod\ 667)$模式，而且这个转变过程是不是可以通过上面两个算法来实现，输入值如果是$(x\cdot1000)$和$(y\cdot1000)$，则通过第一个算法可以得到$((x\cdot1000)\cdot(y\cdot1000)/1000)$，把结果作为第二个算法的输入值，则通过第二个算法可以得到$((x\cdot1000)\cdot(y\cdot1000)/1000)/1000$。

扯了一大顿，终于引出了今天文章的主角，前面讲到的两个算法，第一个就是蒙哥马利乘模，第二个就是蒙哥马利约减。下面我们来讲这两个算法的详解。

正如前面提到的蒙哥马利算法的三个特性之一是，不是基于普通的整数表示法，而是基于蒙哥马利表示法。什么是蒙哥马利表示法呢，其实也很简单，上面我们提到，要让$(x_1\cdot y_1)\ (mod\ 667)$转变成$(x_2\cdot y_2)/1000\ (mod\ 667)$这种形式，我们需要将输入参数变成$(x\cdot1000)$和$(y\cdot1000)$，而不是x和y本身，而$(x\cdot1000)$和$(y\cdot1000)$ 其实就是蒙哥马利表示法。

所以我们先定义几个概念：

• 蒙哥马利参数
  给定一个N，N在b进制（例如，二进制时，b=2）下共有l位，$gcd(N,b)=1$，先预计算以下几个值(这就是前面提到的特性之一，需要预计算）：
• • $\rho = b^k$ 指定一个最小的k，使得$b^k>N$
• • $\omega = -N^{-1} (mod\ \rho)$
    这两个参数是做什么用的呢，你对照前面的演变过程可以猜到$\rho$ 就是前面演变中的1000，而$\omega$ 则是用来计算前面提到的u的。

• 蒙哥马利表示法
  对于x，$0\leqslant x\leqslant N-1$，x的蒙哥马利表示法表示为$x=x\cdot \rho\ (mod\ N)$

蒙哥马利约减

蒙哥马利约减的定义如下
给定一些整数t，蒙哥马利约减的计算结果是$t\cdot \rho^{-1}\ (mod\ N)$

蒙哥马利约减的算法可表示为

蒙哥马利约减可以算作是下面要说的蒙哥马利模乘当$x=1$时的一种特殊形式，。同时它又是蒙哥马利乘模要用到的一部分，这在下一部分讲蒙哥马利乘模的时候有讲到。

蒙哥马利约减可以用来计算某个值得取模操作，比如我们要计算$m(mod\ N)$，只需要将m
的蒙哥马利表示法$m\cdot \rho$作为参数，带入蒙哥马利约减，则计算结果就是$m(mod\ N)$。

蒙哥马利乘模

一个蒙哥马利乘模包括整数乘法和蒙哥马利约减，现在我们有蒙哥马利表示法：

$\hat{x}=x\cdot\rho\ (mod\ N)$
$\hat{y}=y\cdot\rho\ (mod\ N)$

它们相乘的结果是

$t=\hat{x}\cdot\hat{y}$
$\ =(x\cdot\rho)\cdot(y\cdot\rho)$
$\ =(x\cdot y)\cdot\rho^2$

最后，用一次蒙哥马利约减得到结果

$\hat{t}=(x \cdot y) \cdot \rho\ (mod\ N)$

上面我们可以看出，给出的输入参数是$\hat{x}$ 和$\hat{y}$， 得到的结果是$(x \cdot y) \cdot \rho\ (mod\ N)$，所以蒙哥马利乘法也可以写成如下形式，已知输入参数x和y，蒙哥马利乘法是计算$(x \cdot y) \cdot \rho ^ {-1}\ (mod\ N)$

举个例子：
b=10，也就是说在10进制下，N=667
让$b^k>N$的最小的k是3，所以$\rho=b^k=10^3=1000$
$\omega=-N^{-1}\ (mod\ \rho)=-667^{-1}\ (mod\ \rho)=997$

因为$x=421$，所以$\hat{x}=x\cdot\rho(mod\ N)=421\cdot1000(mod\ 667)=123$
因为$y=422$，所以$\hat{y}=y\cdot\rho(mod\ N)=422\cdot1000(mod\ 667)=456$

所以计算$\hat{x}$和$\hat{y}$蒙哥马利乘结果是

$\hat{x}\cdot\hat{y}\cdot\rho^{-1}=(421\cdot1000\cdot422\cdot1000)\cdot1000^{-1}\ (mod\ 667)$
$\qquad\qquad(421\cdot422)\cdot1000\ (mod\ 667)$
$\qquad\qquad(240)\cdot1000\ (mod\ 667)$
$\qquad\qquad547\ (mod\ 667)$

然后总结一下蒙哥马利约减和蒙哥马利乘法的伪代码实现，这个算法其实就是从蒙哥马利预备知识讲到的算法演变来的。

上面的例子用这个算法可以描述为

蒙哥马利算法是一套很完美的算法，为什么这么说呢，你看一开始已知$x$，我们要求$\hat{x}=x \cdot \rho$，这个过程可以通过蒙哥马利乘法本身来计算，输入参数$x$和$\rho^2$，计算结果就是$\hat{x}=x \cdot \rho$。然后在最后，我们知道$\hat{x}=x \cdot \rho$，要求得$x$的时候，同样可以通过蒙哥马利算法本身计算，输入参数$\hat{x}$和$1$，计算结果就是$x$。有没有一种因就是果，果就是因的感觉，这就是为什么说蒙哥马利算法是一套很完美的算法。

蒙哥马利幂模

最后，才说到我们最开始提到的RSA的核心幂模运算，先来看一下普通幂运算的算法是怎么得出来的。

以下资料来自于百度百科快速模幂运算
针对快速模幂运算这一课题，西方现代数学家提出了大量的解决方案，通常都是先将幂模运算转化为乘模运算。
例如求D=C^15%N
由于：a*b % n = (a % n)*(b % n) % 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 就可以了，从左至右或从右至左验证都可以，从左至右会更简洁，
设E=Sum[i=0 to n](E*2**i)，0<=E<=1
则：
D=1
FOR i=n TO 0
D=D*D % N
IF E=1
D=D*C % N
RETURN D这样，模幂运算就转化成了一系列的模乘运算。

算法可以写成如下的形式

如果我们现在用蒙哥马利样式稍作改变，就可以变成如下的形式：

以上就是蒙哥马利算法的全部，通过蒙哥马利算法中的约减运算，我们将大数运算中的模运算变成了移位操作，极大地提高了大数模乘的效率。

但是在以上的算法，可以发现还有两个变量的计算方式不是很清楚，一个是$\omega$，前面说过$\omega = -N^{-1} (mod N)$ ，其实在算法中，我们看到，$omega$仅仅被用来做$\mod b$操作，所以事实上，我们只需要计算$\mod b$即可。

尽管N有可能是合数（因为两个素数的乘积不一定是素数），但通常N和$\rho$（也就是N和b）是互质的，也就是说$N^{\phi(b)}=1(mob\ b)$(费马定理)，$N^{\phi(b)-1}=N^{-1}(mob\ b)$
因为$b=2^\omega$，所以$\phi(b)=2^{(\omega-1)}$，写成算法是这样的

还有一个参数是$\rho^2$，还记得前面说过$\rho$是怎么得出来吗，选定一个最小的$k$，使得$b^k>N$，我们还知道$N$在$b$进制下是$l_N$位，所以当$k=l_N$的时候肯定是符合要求。

$b=2^{\omega}$ 所以$\rho=b^k=({2^{\omega}})^k$

$\rho^2={({2^w})^k)}^2=2^{2\cdot k\cdot \omega}=2^{2\cdot l_N\cdot \omega}$，算法如下

至此整个蒙哥马利算法就全部说完了。通过这个算法，我们可以实现快速幂模。