欢迎关注个人公众号摸鱼范式
转载自：
这篇文章为大家梳理一下整个蒙哥马利算法的本质，蒙哥马利算法并不是一个独立的算法，而是三个相互独立又相互联系的算法集合，其中包括
- 蒙哥马利乘模，是用来计算(xcdot y (mod N))
- 蒙哥马利约减，是用来计算(tcdot ho^{-1} (mod N))
- 蒙哥马利幂模，是用来计算(x^y (mod N))
其中蒙哥马利幂乘是RSA加密算法的核心部分。
基本概念
梳理几个概念，试想一个集合是整数模N之后得到的(Z_N=left{0,1,2,cdots,N-1 ight})
注：N在base-b进制下有(l_N)位。 比如10进制和100进制，都属于base-10进制，因为(100=10^2)，所以b=10。在10进制下，667的(l_N=3)
这样的集合叫做N的剩余类环，任何属于这个集合Z的x满足以下两个条件：
- 正整数
- 最大长度是(l_N)
这篇文章中讲到的蒙哥马利算法就是用来计算基于(Z_N)集合上的运算，简单讲一下原因，因为RSA是基于大数运算的，通常是1024bit或2048bit，而我们的计算机不可能存储完整的大数，因为占空间太大，而且也没必要。因此，这种基于大数运算的加密体系在计算的时候都是基于(Z_N)集合的，自然，蒙哥马利算法也是基于(Z_N)。
在剩余类环上，有两种重要的运算，一类是简单运算，也就是加法和减法，另一类复杂运算，也就是乘法。我们比较熟悉的是自然数集上的运算，下面看下怎么从自然数集的运算演变成剩余类环上的运算。
对于加法运算，如果计算(xpm y (mod N)),(0leqslant x,y)，试想自然数集上的 (xpm y)
我们可以简单的通过加减N来实现从自然数到剩余类集的转换
另外一类是乘法操作，也就是(xcdot y (mod N)),(0leqslant x,y)，那么
如果在自然数集下，令(t=xcdot y)，那么对于(mod N)我们需要计算
加减操作很简单，具体的算这里就不细说了，我们用(Z_N-ADD)来代表剩余类环上的加法操作。既然我们可以做加法操作，那么我们就可以扩展到乘法操作，算法如下
但是这并不是一个好的解决方案，因为通常来说，我们不会直接做w位乘w位的操作，这个后面会用蒙哥马利的乘法来代替解决。
对于取模操作，一般有以下几种方法
- 根据以下公式，来计算取模操作
- 整个计算过程是基于标准的数字表示
- 不需要预计算（也就是提前计算一些变量，以备使用）
- 涉及到一个除法操作，非常费时和复杂
- 用Barrett reduction算法，这篇文章不细说，但是有以下特征
- 基于标准的数字表示
- 不需要预计算
- 需要(2 cdot (l_N+1) cdot (l_N+1))次数乘运算
- 用蒙哥马利约减，也就是下面要讲的算法，有以下特征
- 不是基于标准的数字表示（后文中有提到，是基于蒙哥马利表示法）
- 需要预计算
- 需要(2 cdot (l_N) cdot (l_N))次数乘运算
蒙哥马利预备知识
在将蒙哥马利算法之前，先看一下在自然数下的乘法公式
计算(xcdot y)，想象一下我们常用的计算乘法的方法，用乘数的每一位乘上被乘数，然后把得到的结果相加，总结成公式，可以写成如下的形式。
尝试下面一个例子，10进制下（也就是b=10），y=456（也就是(l_n=3)），计算(xcdot y)，公式可演变如下：
最后一次演变其实就是用霍纳法则(Horner’s rule)所讲的规则，关于霍纳法则，可以自行百度。
这个计算过程写成代码实现的算法应该是这样的：
接下来我们来看下面这样的计算，计算(x⋅y)/1000(x⋅y)/1000，由前面可以知道
由此可以知道：
这个计算过程写成代码实现的算法是这样的：
接下来我们再来看在剩余类集合下的乘法操作 (xcdot y/1000 (mod 667))
我们知道剩余类集合(Z_{667}=left{0,1 cdots 666 ight})，是不存在小数的，而如果我们采用自然数集的计算方式的话，就会出现小数，比如前面的例子，除10就会有小数。
这个问题是这样的，我们知道(u·667 equiv 0 (mod 667))（≡表示取模相等），所以我们可以选择一个合适的u，用u乘667再加上r，使得和是一个可以除10没有小数，这样在mod 667之后依然是正确的结果。至于u怎么算出来，这篇文章会在后面的章节说明。
这个过程之后(xcdot y/1000 (mod 667)) 的计算算法可以写成如下的形式
至此，你可能还不明白上面说这一堆演变的原因，其实很简单，原来是一个((xcdot y) (mod 667))的运算，这个运算中的模操作，正常情况下是要通过除法实现的，而除法是一个特别复杂的运算，要涉及到很多乘法，所以在大数运算时，我们要尽量避免除法的出现。而通过以上几个步骤，我们发现((xcdot y)/1000 (mod 667))这个操作是不用除法的。等等，算法中明明有个除10的操作，你骗谁呢。不知道你有没有发现，除数其实是我们的进制数，除进制数在计算机中是怎么做呢，其实很简单，左移操作就ok了。所以这个计算方法是不涉及到除法操作的。
但是我们要计算的明明是((x_1cdot y_1) (mod 667))，怎么现在变成了((x_2cdot y_2)/1000 (mod 667))，所以在下一步，我们要思考的是怎么样让((x_1cdot y_1) (mod 667))转变成((x_2cdot y_2)/1000 (mod 667))这种形式。
考虑这样两个算法
- 第一个是输入x和y，计算(x cdot y cdot ho^{-1})
- 第二个算法，输入一个t，计算(t cdot ho^{-1})。
是不是变成了我们需要的((xcdot y)/1000 (mod 667))模式，而且这个转变过程是不是可以通过上面两个算法来实现，输入值如果是((xcdot1000))和((ycdot1000))，则通过第一个算法可以得到(((xcdot1000)cdot(ycdot1000)/1000))，把结果作为第二个算法的输入值，则通过第二个算法可以得到(((xcdot1000)cdot(ycdot1000)/1000)/1000)。
扯了一大顿，终于引出了今天文章的主角，前面讲到的两个算法，第一个就是蒙哥马利乘模，第二个就是蒙哥马利约减。下面我们来讲这两个算法的详解。
正如前面提到的蒙哥马利算法的三个特性之一是，不是基于普通的整数表示法，而是基于蒙哥马利表示法。什么是蒙哥马利表示法呢，其实也很简单，上面我们提到，要让((x_1cdot y_1) (mod 667))转变成((x_2cdot y_2)/1000 (mod 667))这种形式，我们需要将输入参数变成((xcdot1000))和((ycdot1000))，而不是x和y本身，而(((xcdot1000)cdot(ycdot1000)/1000)/1000) 其实就是蒙哥马利表示法。
所以我们先定义几个概念：
-
蒙哥马利参数
给定一个N，N在b进制（例如，二进制时，b=2）下共有l位，(gcd(N,b)=1)，先预计算以下几个值(这就是前面提到的特性之一，需要预计算）： -
( ho = b^k) 指定一个最小的k，使得(b^k>N)
(omega = -N^{-1} (mod ho))
这两个参数是做什么用的呢，你对照前面的演变过程可以猜到( ho)就是前面演变中的1000，而(omega)则是用来计算前面提到的u的。 -
蒙哥马利表示法
对于x，(0leqslant xleqslant N-1)，x的蒙哥马利表示法表示为(x=xcdot ho (mod N))
蒙哥马利约减
蒙哥马利约减的定义如下
给定一整数t，蒙哥马利约减的计算结果是(tcdot ho^{-1} (mod N))
蒙哥马利约减的算法可表示为
蒙哥马利约减可以算作是下面要说的蒙哥马利模乘当(x=1)时的一种特殊形式，。同时它又是蒙哥马利乘模要用到的一部分，这在下一部分讲蒙哥马利乘模的时候有讲到。
蒙哥马利约减可以用来计算某个值得取模操作，比如我们要计算(m(mod N))，只需要将m的蒙哥马利表示法(mcdot ho)作为参数，带入蒙哥马利约减，则计算结果就(m(mod N))。
蒙哥马利乘模
一个蒙哥马利乘模包括整数乘法和蒙哥马利约减，现在我们有蒙哥马利表示法：
它们相乘的结果是
最后，用一次蒙哥马利约减得到结果
上面我们可以看出，给出的输入参数是(hat{x}) 和(hat{y})， 得到的结果是((x cdot y) cdot ho (mod N))，所以蒙哥马利乘法也可以写成如下形式，已知输入参数x和y，蒙哥马利乘法是计算((x cdot y) cdot ho ^ {-1} (mod N))
举个例子：
b=10，也就是说在10进制下，N=667
让(b^k>N)的最小的k是3，所以( ho=b^k=10^3=1000)
(omega=-N^{-1} (mod ho)=-667^{-1} (mod ho)=997)
因为(x=421)，所以(hat{x}=xcdot ho(mod N)=421cdot1000(mod 667)=123)
因为(y=422)，所以(hat{y}=ycdot ho(mod N)=422cdot1000(mod 667)=456)
所以计算(hat{x}) 和(hat{y})蒙哥马利乘结果是
然后总结一下蒙哥马利约减和蒙哥马利乘法的伪代码实现，这个算法其实就是从蒙哥马利预备知识讲到的算法演变来的。
上面的例子用这个算法可以描述为
蒙哥马利算法是一套很完美的算法，为什么这么说呢，你看一开始已知(x)，我们要求(hat{x}=x cdot ho)，这个过程可以通过蒙哥马利乘法本身来计算，输入参数(x)和( ho^2)，计算结果就是(hat{x}=x cdot ho)。然后在最后，我们知道(hat{x}=x cdot ho)，要求得(x)的时候，同样可以通过蒙哥马利算法本身计算，输入参数(hat{x})和(1)，计算结果就是(x)。有没有一种因就是果，果就是因的感觉，这就是为什么说蒙哥马利算法是一套很完美的算法。
蒙哥马利幂模
最后，才说到我们最开始提到的RSA的核心幂模运算，先来看一下普通幂运算的算法是怎么得出来的。
以下资料来自于百度百科
针对快速模幂运算这一课题，西方现代数学家提出了大量的解决方案，通常都是先将幂模运算转化为乘模运算。
例如求D=C^15%N
由于：ab % 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 b)) ，其实在算法中，我们看到，(omega)仅仅被用来做(mod b)操作，所以事实上，我们只需要计算(mod b)即可。
尽管N有可能是合数（因为两个素数的乘积不一定是素数），但通常N和( ho)（也就是N和b）是互质的，也就是说(N^{phi(b)}=1(mob b))(费马定理)，(N^{phi(b)-1}=N^{-1}(mob b)),因为(b=2^omega)，所以(phi(b)=2^{(omega-1)})，写成算法是这样的
还有一个参数是( ho^2)，还记得前面说过( ho)是怎么得出来吗，选定一个最小的(k)，使得(b^k>N)，我们还知道(N)在(b)进制下是(l_N)位，所以当(k=l_N)的时候肯定是符合要求。
(b=2^{omega}) 所以( ho=b^k=({2^{omega}})^k)
( ho^2={({2^w})^k)}^2=2^{2cdot kcdot omega}=2^{2cdot l_Ncdot omega})，算法如下
至此整个蒙哥马利算法就全部说完了。通过这个算法，我们可以实现快速幂模。