蒙哥馬利（Montgomery）冪模運(yùn)算是快速計(jì)算a^b%k的一種算法，是RSA加密算法的核心之一。

　　蒙哥馬利模乘的優(yōu)點(diǎn)在于減少了取模的次數(shù)（在大數(shù)的條件下）以及簡(jiǎn)化了除法的復(fù)雜度（在2的k次冪的進(jìn)制下除法僅需要進(jìn)行左移操作）。模冪運(yùn)算是RSA 的核心算法，最直接地決定了RSA 算法的性能。

　　針對(duì)快速模冪運(yùn)算這一課題，西方現(xiàn)代數(shù)學(xué)家提出了大量的解決方案，通常都是先將冪模運(yùn)算轉(zhuǎn)化為乘模運(yùn)算。

　　這里為大家梳理一下整個(gè)蒙哥馬利算法的本質(zhì)，蒙哥馬利算法并不是一個(gè)獨(dú)立的算法，而是三個(gè)相互獨(dú)立又相互聯(lián)系的算法集合，其中包括

　　蒙哥馬利乘模，是用來(lái)計(jì)算x?y （mod N）

　　蒙哥馬利約減，是用來(lái)計(jì)算t?ρ?1 （mod N）

　　蒙哥馬利冪模，是用來(lái)計(jì)算xy （mod N）

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

　　基本概念

　　梳理幾個(gè)概念，試想一個(gè)集合是整數(shù)模N之后得到的

　　ZN={0，1，2，?，N?1}

　　注：N在base-b進(jìn)制下有l(wèi)N位。 比如10進(jìn)制和100進(jìn)制，都屬于base-10進(jìn)制，因?yàn)?00=102，所以b=10。在10進(jìn)制下，667的lN=3這樣的集合叫做N的剩余類(lèi)環(huán)，任何屬于這個(gè)集合Z的x滿足以下兩個(gè)條件：

　　1. 正整數(shù)

　　2. 最大長(zhǎng)度是lN

　　文中講到的蒙哥馬利算法就是用來(lái)計(jì)算基于ZN集合上的運(yùn)算，簡(jiǎn)單講一下原因，因?yàn)镽SA是基于大數(shù)運(yùn)算的，通常是1024bit或2018bit，而我們的計(jì)算機(jī)不可能存儲(chǔ)完整的大數(shù)，因?yàn)檎伎臻g太大，而且也沒(méi)必要。因此，這種基于大數(shù)運(yùn)算的加密體系在計(jì)算的時(shí)候都是基于ZN集合的，自然，蒙哥馬利算法也是基于ZN。

　　在剩余類(lèi)環(huán)上，有兩種重要的運(yùn)算，一類(lèi)是簡(jiǎn)單運(yùn)算，也就是加法和減法，另一類(lèi)復(fù)雜運(yùn)算，也就是乘法。我們比較熟悉的是自然數(shù)集上的運(yùn)算，下面看下怎么從自然數(shù)集的運(yùn)算演變成剩余類(lèi)環(huán)上的運(yùn)算。

　　對(duì)于加法運(yùn)算，如果計(jì)算x±y （mod N） （0≤x，y＜N），試想自然數(shù)集上的 x±y

　　0≤x+y≤2?（N?1）

　　?（N?1）≤x?y≤（N?1）我們可以簡(jiǎn)單的通過(guò)加減N來(lái)實(shí)現(xiàn)從自然數(shù)到剩余類(lèi)集的轉(zhuǎn)換

　　另外一類(lèi)是乘法操作，也就是x?y （mod N）（0≤x，y＜N），那么

　　0≤x?y≤（N?1）2如果在自然數(shù)集下，令t=x?y，那么對(duì)于modN我們需要計(jì)算

　　t?（N??t/N?）加減操作很簡(jiǎn)單，具體的算這里就不細(xì)說(shuō)了，我們用ZN?ADD 來(lái)代表剩余類(lèi)環(huán)上的加法操作。既然我們可以做加法操作，那么我們就可以擴(kuò)展到乘法操作，算法如下

　　但是這并不是一個(gè)好的解決方案，因?yàn)橥ǔ?lái)說(shuō)，我們不會(huì)直接做w位乘w位的操作，這個(gè)后面會(huì)用蒙哥馬利的乘法來(lái)代替解決。

　　對(duì)于取模操作，一般有以下幾種方法

　　1，根據(jù)以下公式，來(lái)計(jì)算取模操作

　　t?（N??t/N?）

　　這種解法有以下特征

　　整個(gè)計(jì)算過(guò)程是基于標(biāo)準(zhǔn)的數(shù)字表示

　　不需要預(yù)計(jì)算（也就是提前計(jì)算一些變量，以備使用）

　　涉及到一個(gè)除法操作，非常費(fèi)時(shí)和復(fù)雜

　　2，用Barrett reduction算法，這篇文章不細(xì)說(shuō)，但是有以下特征

　　基于標(biāo)準(zhǔn)的數(shù)字表示

　　不需要預(yù)計(jì)算

　　需要2?（lN+1）?（lN+1） 次數(shù)乘運(yùn)算

　　3，用蒙哥馬利約減，也就是下面要講的算法，有以下特征

　　不是基于標(biāo)準(zhǔn)的數(shù)字表示（后文中有提到，是基于蒙哥馬利表示法）

　　需要預(yù)計(jì)算

　　需要2?（lN）?（lN） 次數(shù)乘運(yùn)算

　　蒙哥馬利預(yù)備知識(shí)

　　在將蒙哥馬利算法之前，先看一下在自然數(shù)下的乘法公式

　　計(jì)算x?y，想象一下我們常用的計(jì)算乘法的方法，用乘數(shù)的每一位乘上被乘數(shù)，然后把得到的結(jié)果相加，總結(jié)成公式，可以寫(xiě)成如下的形式。

　　x?y=x?sumly?1i=0yi?bi

　　=sumly?1i=0yi?x?bi

　　嘗試下面一個(gè)例子，10進(jìn)制下（也就是b=10），y=456（也就是ln=3），計(jì)算x?y，公式可演變?nèi)缦拢?/p>

　　x?y=（y0?x?100）+（y1?x?101）+（y2?x?102）

　　=（y0?x?0）+（y1?x?10）+（y2?x?100）

　　=（y0?x）+10?（y1?x+10?（y2?x?+10?（0）））

　　最后一次演變其實(shí)就是用霍納法則（Horner’s rule）所講的規(guī)則，關(guān)于霍納法則，可以自行百度。

　　這個(gè)計(jì)算過(guò)程寫(xiě)成代碼實(shí)現(xiàn)的算法應(yīng)該是這樣的：

　　接下來(lái)我們來(lái)看下面這樣的計(jì)算，計(jì)算（x?y）/1000，由前面可以知道

　　x?y=（y_0?x）+10?（y_1?x+10?（y_2?x?+10?（0）））

　　由此可以知道：

　　這個(gè)計(jì)算過(guò)程寫(xiě)成代碼實(shí)現(xiàn)的算法是這樣的：

　　接下來(lái)我們?cè)賮?lái)看在剩余類(lèi)集合下的乘法操作 x?y/1000 （mod 667）

　　我們知道剩余類(lèi)集合Z667={0，1?666}，是不存在小數(shù)的，而如果我們采用自然數(shù)集的計(jì)算方式的話，就會(huì)出現(xiàn)小數(shù)，比如前面的例子，除10就會(huì)有小數(shù)。

　　這個(gè)問(wèn)題是這樣的，我們知道u?667≡0（mod667）（≡表示取模相等），所以我們可以選擇一個(gè)合適的u，用u乘667再加上r，使得和是一個(gè)可以除10沒(méi)有小數(shù)，這樣在mod 667之后依然是正確的結(jié)果。至于u怎么算出來(lái)，這篇文章會(huì)在后面的章節(jié)說(shuō)明。

　　這個(gè)過(guò)程之后x?y/1000 （mod 667） 的計(jì)算算法可以寫(xiě)成如下的形式

　　至此，你可能還不明白上面說(shuō)這一堆演變的原因，其實(shí)很簡(jiǎn)單，原來(lái)是一個(gè)（x?y） （mod 667）的運(yùn)算，這個(gè)運(yùn)算中的模操作，正常情況下是要通過(guò)除法實(shí)現(xiàn)的，而除法是一個(gè)特別復(fù)雜的運(yùn)算，要涉及到很多乘法，所以在大數(shù)運(yùn)算時(shí)，我們要盡量避免除法的出現(xiàn)。而通過(guò)以上幾個(gè)步驟，我們發(fā)現(xiàn)（x?y）/1000 （mod 667）這個(gè)操作是不用除法的。等等，算法中明明有個(gè)除10的操作，你騙誰(shuí)呢。不知道你有沒(méi)有發(fā)現(xiàn)，除數(shù)其實(shí)是我們的進(jìn)制數(shù)，除進(jìn)制數(shù)在計(jì)算機(jī)中是怎么做呢，其實(shí)很簡(jiǎn)單，左移操作就ok了。所以這個(gè)計(jì)算方法是不涉及到除法操作的。

　　但是我們要計(jì)算的明明是（x1?y1） （mod 667），怎么現(xiàn)在變成了（x2?y2）/1000 （mod 667），所以在下一步，我們要思考的是怎么樣讓?zhuān)▁1?y1） （mod 667）轉(zhuǎn)變成（x2?y2）/1000 （mod 667）這種形式。

　　考慮這樣兩個(gè)算法

　　- 第一個(gè)是輸入x和y，計(jì)算x?y?ρ?1

　　- 第二個(gè)算法，輸入一個(gè)t，計(jì)算t?ρ?1。

　　x?y （mod 667）=（（x?1000）?（y?1000）/1000）/1000 （mod 667）

　　是不是變成了我們需要的（x?y）/1000 （mod 667）模式，而且這個(gè)轉(zhuǎn)變過(guò)程是不是可以通過(guò)上面兩個(gè)算法來(lái)實(shí)現(xiàn)，輸入值如果是（x?1000）和（y?1000），則通過(guò)第一個(gè)算法可以得到（（x?1000）?（y?1000）/1000），把結(jié)果作為第二個(gè)算法的輸入值，則通過(guò)第二個(gè)算法可以得到（（x?1000）?（y?1000）/1000）/1000。

? ? ? ? 算法的詳解

　　扯了一大頓，終于引出了今天文章的主角，前面講到的兩個(gè)算法，第一個(gè)就是蒙哥馬利乘模，第二個(gè)就是蒙哥馬利約減。下面我們來(lái)講這兩個(gè)算法的詳解。

　　正如前面提到的蒙哥馬利算法的三個(gè)特性之一是，不是基于普通的整數(shù)表示法，而是基于蒙哥馬利表示法。什么是蒙哥馬利表示法呢，其實(shí)也很簡(jiǎn)單，上面我們提到，要讓?zhuān)▁1?y1） （mod 667）轉(zhuǎn)變成（x2?y2）/1000 （mod 667）這種形式，我們需要將輸入參數(shù)變成（x?1000）和（y?1000），而不是x和y本身，而（x?1000）和（y?1000） 其實(shí)就是蒙哥馬利表示法。

　　所以我們先定義幾個(gè)概念：

　　蒙哥馬利參數(shù)

　　給定一個(gè)N，N在b進(jìn)制（例如，二進(jìn)制時(shí)，b=2）下共有l(wèi)位，gcd（N，b）=1，先預(yù)計(jì)算以下幾個(gè)值（這就是前面提到的特性之一，需要預(yù)計(jì)算）：

　　ρ=bk 指定一個(gè)最小的k，使得bk》N

　　ω=?N?1（mod ρ）

　　這兩個(gè)參數(shù)是做什么用的呢，你對(duì)照前面的演變過(guò)程可以猜到ρ 就是前面演變中的1000，而ω 則是用來(lái)計(jì)算前面提到的u的。

　　蒙哥馬利表示法

　　對(duì)于x，0≤x≤N?1，x的蒙哥馬利表示法表示為x=x?ρ （mod N）

　　蒙哥馬利約減

　　蒙哥馬利約減的定義如下

　　給定一些整數(shù)t，蒙哥馬利約減的計(jì)算結(jié)果是t?ρ?1 （mod N）

　　蒙哥馬利約減的算法可表示為

　　蒙哥馬利約減可以算作是下面要說(shuō)的蒙哥馬利模乘當(dāng)x=1時(shí)的一種特殊形式，。同時(shí)它又是蒙哥馬利乘模要用到的一部分，這在下一部分講蒙哥馬利乘模的時(shí)候有講到。

　　蒙哥馬利約減可以用來(lái)計(jì)算某個(gè)值得取模操作，比如我們要計(jì)算m（mod N），只需要將m

　　的蒙哥馬利表示法m?ρ作為參數(shù)，帶入蒙哥馬利約減，則計(jì)算結(jié)果就是m（mod N）。

　　蒙哥馬利乘模

　　一個(gè)蒙哥馬利乘模包括整數(shù)乘法和蒙哥馬利約減，現(xiàn)在我們有蒙哥馬利表示法：

　　x^=x?ρ （mod N）

　　y^=y?ρ （mod N）

　　它們相乘的結(jié)果是

　　t=x^?y^

　　=（x?ρ）?（y?ρ）

　　=（x?y）?ρ2

　　最后，用一次蒙哥馬利約減得到結(jié)果

　　t^=（x?y）?ρ （mod N）

　　上面我們可以看出，給出的輸入?yún)?shù)是x^ 和y^， 得到的結(jié)果是（x?y）?ρ （mod N），所以蒙哥馬利乘法也可以寫(xiě)成如下形式，已知輸入?yún)?shù)x和y，蒙哥馬利乘法是計(jì)算（x?y）?ρ?1 （mod N）

　　舉個(gè)例子：

　　b=10，也就是說(shuō)在10進(jìn)制下，N=667

　　讓bk》N的最小的k是3，所以ρ=bk=103=1000

　　ω=?N?1 （mod ρ）=?667?1 （mod ρ）=997

　　因?yàn)閤=421，所以x^=x?ρ（mod N）=421?1000（mod 667）=123

　　因?yàn)閥=422，所以y^=y?ρ（mod N）=422?1000（mod 667）=456

　　所以計(jì)算x^和y^蒙哥馬利乘結(jié)果是

　　x^?y^?ρ?1=（421?1000?422?1000）?1000?1 （mod 667）

　　（421?422）?1000 （mod 667）

　　（240）?1000 （mod 667）

　　547 （mod 667）

　　然后總結(jié)一下蒙哥馬利約減和蒙哥馬利乘法的偽代碼實(shí)現(xiàn)，這個(gè)算法其實(shí)就是從蒙哥馬利預(yù)備知識(shí)講到的算法演變來(lái)的。

　　上面的例子用這個(gè)算法可以描述為

　　蒙哥馬利算法是一套很完美的算法，為什么這么說(shuō)呢，你看一開(kāi)始已知x，我們要求x^=x?ρ，這個(gè)過(guò)程可以通過(guò)蒙哥馬利乘法本身來(lái)計(jì)算，輸入?yún)?shù)x和ρ2，計(jì)算結(jié)果就是x^=x?ρ。然后在最后，我們知道x^=x?ρ，要求得x的時(shí)候，同樣可以通過(guò)蒙哥馬利算法本身計(jì)算，輸入?yún)?shù)x^和1，計(jì)算結(jié)果就是x。有沒(méi)有一種因就是果，果就是因的感覺(jué)，這就是為什么說(shuō)蒙哥馬利算法是一套很完美的算法。

　　蒙哥馬利冪模

　　最后，才說(shuō)到我們最開(kāi)始提到的RSA的核心冪模運(yùn)算，先來(lái)看一下普通冪運(yùn)算的算法是怎么得出來(lái)的。

　　以下資料來(lái)自于百度百科快速模冪運(yùn)算

　　針對(duì)快速模冪運(yùn)算這一課題，西方現(xiàn)代數(shù)學(xué)家提出了大量的解決方案，通常都是先將冪模運(yùn)算轉(zhuǎn)化為乘模運(yùn)算。

　　例如求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

　　即：對(duì)于E=15的冪模運(yùn)算可分解為6 個(gè)乘模運(yùn)算，歸納分析以上方法可以發(fā)現(xiàn)：

　　對(duì)于任意指數(shù)E，都可采用以下算法計(jì)算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

　　繼續(xù)分析會(huì)發(fā)現(xiàn)，要知道E 何時(shí)能整除 2，并不需要反復(fù)進(jìn)行減一或除二的操作，只需驗(yàn)證E 的二進(jìn)制各位是0 還是1 就可以了，從左至右或從右至左驗(yàn)證都可以，從左至右會(huì)更簡(jiǎn)潔，

　　設(shè)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這樣，模冪運(yùn)算就轉(zhuǎn)化成了一系列的模乘運(yùn)算。

　　算法可以寫(xiě)成如下的形式

　　如果我們現(xiàn)在用蒙哥馬利樣式稍作改變，就可以變成如下的形式：

　　以上就是蒙哥馬利算法的全部，通過(guò)蒙哥馬利算法中的約減運(yùn)算，我們將大數(shù)運(yùn)算中的模運(yùn)算變成了移位操作，極大地提高了大數(shù)模乘的效率。

　　但是在以上的算法，可以發(fā)現(xiàn)還有兩個(gè)變量的計(jì)算方式不是很清楚，一個(gè)是ω，前面說(shuō)過(guò)ω=?N?1（modN） ，其實(shí)在算法中，我們看到，omega僅僅被用來(lái)做modb操作，所以事實(shí)上，我們只需要計(jì)算modb即可。

　　盡管N有可能是合數(shù)（因?yàn)閮蓚€(gè)素?cái)?shù)的乘積不一定是素?cái)?shù)），但通常N和ρ（也就是N和b）是互質(zhì)的，也就是說(shuō)N?（b）=1（mob b）（費(fèi)馬定理），N?（b）?1=N?1（mob b）

　　因?yàn)閎=2ω，所以?（b）=2（ω?1），寫(xiě)成算法是這樣的

　　還有一個(gè)參數(shù)是ρ2，還記得前面說(shuō)過(guò)ρ是怎么得出來(lái)嗎，選定一個(gè)最小的k，使得bk》N，我們還知道N在b進(jìn)制下是lN位，所以當(dāng)k=lN的時(shí)候肯定是符合要求。

　　b=2ω 所以ρ=bk=（2ω）kρ2=（2w）k）2=22?k?ω=22?lN?ω，算法如下

　　至此整個(gè)蒙哥馬利算法就全部說(shuō)完了。通過(guò)這個(gè)算法，我們可以實(shí)現(xiàn)快速冪模。

　　C++實(shí)現(xiàn)

　　inline unsigned __int64 MulMod（unsigned __int64 a， unsigned __int64 b， unsigned __int64 n）

　　{

　　return a * b % n;

　　}

　　/*

　　模冪運(yùn)算，返回值 x=base^pow mod n

　　*/

　　unsigned __int64 PowMod（unsigned __int64 base， unsigned __int64 pow， unsigned __int64 &n）

　　{

　　unsigned __int64 a=base， b=pow， c=1;

　　while（b）

　　{

　　while（！（b & 1））

　　{

　　b》》=1; //a=a * a % n; //

　　a=MulMod（a， a， n）;

　　} b--; //c=a * c % n; //

　　c=MulMod（a， c， n）;

　　} return c;

　　}

? ? ? ?數(shù)論學(xué)家利用費(fèi)馬小定理研究出了多種素?cái)?shù)測(cè)試方法，目前最快的算法是拉賓

　　米勒測(cè)試算法，其過(guò)程如下：

　　（1）計(jì)算奇數(shù)M，使得N=（2**r）*M+1

　　（2）選擇隨機(jī)數(shù)A《N

　　（3）對(duì)于任意i《r，若A**（（2**i）*M） % N = N-1，則N通過(guò)隨機(jī)數(shù)A的測(cè)試

　　（4）或者，若A**M % N = 1，則N通過(guò)隨機(jī)數(shù)A的測(cè)試

　　（5）讓A取不同的值對(duì)N進(jìn)行5次測(cè)試，若全部通過(guò)則判定N為素?cái)?shù)

　　若N 通過(guò)一次測(cè)試，則N 不是素?cái)?shù)的概率為 25%，若N 通過(guò)t 次測(cè)試，則N 不是素?cái)?shù)的概率為1/4**t。事實(shí)上取t 為5 時(shí)，N 不是素?cái)?shù)的概率為 1/128，N 為素?cái)?shù)的概率已經(jīng)大于99.99%。

　　在實(shí)際應(yīng)用中，可首先用300—500個(gè)小素?cái)?shù)對(duì)N 進(jìn)行測(cè)試，以提高拉賓米勒測(cè)試通過(guò)的概率，從而提高整體測(cè)試速度。而在生成隨機(jī)素?cái)?shù)時(shí)，選取的隨機(jī)數(shù)最好讓r=0，則可省去步驟（3） 的測(cè)試，進(jìn)一步提高測(cè)試速度。

　　素?cái)?shù)測(cè)試是RSA 選取秘鑰的第一步，奇妙的是其核心運(yùn)算與加解密時(shí)所需的運(yùn)算完全一致：都是模冪運(yùn)算。而模冪運(yùn)算過(guò)程中中所需求解的歐幾里德方程又恰恰正是選取密鑰第二步所用的運(yùn)算。可見(jiàn)整個(gè)RSA 算法具有一種整體的和諧。