蒙哥马利算法(Montgomery Algorithm)|蒙哥马利约简、模乘、模幂

Montgomery Algorithm(蒙哥马利算法)

蒙哥马利算法分为3种，蒙哥马利模乘，蒙哥马利约简，蒙哥马利模幂

1、从蒙哥马利模乘说起

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

• 蒙哥马利形式：为了计算$ab(\mathrm{m}\mathrm{o}\mathrm{d}N)$，找一个$R$，然后使得$a^{\prime }\equiv aR(\mathrm{m}\mathrm{o}\mathrm{d}N),b^{\prime }\equiv bR(\mathrm{m}\mathrm{o}\mathrm{d}N)$，这就是蒙哥马利形式。
• 这个$R$不是随便一个数，需要满足两个条件：1:$R=2^k>N$,其中$k$是满足条件的最小的正数，这个取法能够保证除以$R$就相当于右移k位，避免除法运算；2:$gcd(R,N)=1$，这使得能够求出下面的m。
• 蒙哥马利形式有什么用？往下看

前面只是铺垫，下面才真正开始。为了计算一开始的$ab(\mathrm{m}\mathrm{o}\mathrm{d}N)$，需要用到上面的蒙哥马利形式。令$X=a^{\prime }{b}^{\prime }$，我们可以设计一个函数来计算$X{R}^{-1}(\mathrm{m}\mathrm{o}\mathrm{d}N)$，简单计算发现这个函数的计算结果为$X_1\equiv X{R}^{-1}\equiv a^{\prime }{b}^{\prime }{R}^{-1}\equiv abR(\mathrm{m}\mathrm{o}\mathrm{d}N)$，这样在调用一遍函数计算$X_1{R}^{-1}(\mathrm{m}\mathrm{o}\mathrm{d}N)$就得到我们最终需要的结果$ab(\mathrm{m}\mathrm{o}\mathrm{d}N)$了。我们称这个算法叫蒙哥马利约简算法。所以说，蒙哥马利约简的产生是为了蒙哥马利模乘计算服务的。

2、蒙哥马利约简

根据上述分析可知蒙哥马利算法的核心在于蒙哥马利约简。而且前面提到，蒙哥马利算法的主要思想是把取模运算变得简单，到底是怎么做到的呢？

主要思想：蒙哥马利约简是计算$X{R}^{-1}(\mathrm{m}\mathrm{o}\mathrm{d}N)$，这相当于$\frac{X}{R}(\mathrm{m}\mathrm{o}\mathrm{d}N)$，如何才能避免除法？从前面的$R$的定义中我们知道，$R=2^k$，所以$\frac{X}{R}=X>>k$ （X右移k位）。但是新的问题出现了，右移k位可能会抹掉X的低位中的一些1，如$7÷4=0b111>>2=0b1=1$，这个不是精确计算，而是向下取整的除法。当且仅当X是R的整数倍时，$X÷R$严格等于$X>>k$。所以我们实际上是找一个$m$，使得$X+mN$是$R$的倍数，这样计算$\frac{X+mN}{R}(\mathrm{m}\mathrm{o}\mathrm{d}N)$就可以了。

m如何找

根据R的定义，$gcd(R,N)=1$，根据扩展的欧几里得算法，有$R{R}^{\prime }-N{N}^{\prime }=1$并且有$0<N^{\prime }<R,0<R^{\prime }<N<R$。（注意这里是$-N{N}^{\prime }$，所以有$N^{\prime }=-N^{-1}(\mathrm{m}\mathrm{o}\mathrm{d}R)$）
$$\begin{gathered} X+mN\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}R) \\ X{N}^{\prime }+mN{N}^{\prime }\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}R) \\ X{N}^{\prime }+m(R{R}^{\prime }-1)\equiv 0(\mathrm{m}\mathrm{o}\mathrm{d}R) \\ X{N}^{\prime }\equiv m(\mathrm{m}\mathrm{o}\mathrm{d}R) \end{gathered}$$
这样就求出了m。

约简算法总流程

提前工作：已知$a,b,N$，并计算出蒙哥马利形式中的$a^{\prime },b^{\prime },R$以及$X$。

Montgomery reduction$(X,R,N)$：

1. 计算$N^{\prime }=-N^{-1}(\mathrm{m}\mathrm{o}\mathrm{d}R)$，计算$m=X{N}^{\prime }(\mathrm{m}\mathrm{o}\mathrm{d}R)$；

2. 计算$y=\frac{X+mN}{R}$:将$X+mN$右移$k$位；

3. 若$y>N$，则$y=y-N$；

   ​ 这时的$y$满足：$0<y<2N$。因为
   $$X<N^2,m<R,N<R$$
   ​ 所以
   $$\frac{X+mN}{R}<\frac{N^2+R\cdot N}{R}<\frac{RN+R\cdot N}{R}=2N$$

4. 返回y。

3、蒙哥马利模乘流程

Montgomery Multiply$(a,b,N)$:

1. 计算$a^{\prime }\equiv aR(\mathrm{m}\mathrm{o}\mathrm{d}N),b^{\prime }\equiv bR(\mathrm{m}\mathrm{o}\mathrm{d}N)，X=a^{\prime }{b}^{\prime }$；
2. 调用蒙哥马利约简算法：$X_1$=Montgomery reduction$(X,R,N)=X/R=abR(\mathrm{m}\mathrm{o}\mathrm{d}N)$；
3. 再调用蒙哥马利约简算法：y=Montgomery reduction$(X_1,R,N)=X_1/R=ab(\mathrm{m}\mathrm{o}\mathrm{d}N)$；
4. 返回y。

4、复杂度分析

蒙哥马利算法是为了简化模$N$的复杂度，当$N$是比较大的数时，模N需要多次的加、减、乘运算。从上述过程看来，蒙哥马利约简的复杂度确实降低了，因为只有模$R$的移位运算。但是看蒙哥马利模乘的流程，在第一步中进行蒙哥马利表示时就计算了两次模N，（$a^{\prime }\equiv aR(\mathrm{m}\mathrm{o}\mathrm{d}N),b^{\prime }\equiv bR(\mathrm{m}\mathrm{o}\mathrm{d}N)，X=a^{\prime }{b}^{\prime }$），总的来看复杂度也没有降低啊？实际上，第一步可以看成是蒙哥马利的预先计算，在硬件实现中，先把预先计算的算好，在后面运行就会快很多。尤其是当出现大量的模乘运算时，可以通过并行运算进行预计算，这会大量节省运行时间。

5、蒙哥马利模幂

当进行模幂运算时，也可以利用蒙哥马利算法。如计算$a^e(\mathrm{m}\mathrm{o}\mathrm{d}N)$。根据重复平方-乘方法，模幂运算中就有很多步的模乘运算，这恰恰可以发挥蒙哥马利模乘的优势。具体的，在重复平方法的每一次模乘中都利用蒙哥马利模乘进行计算，把需要的参数提前计算好就行。