蒙哥馬利算法

歡迎關注個人公眾號摸魚范式

轉載自：

蒙哥馬利算法

這篇文章為大家梳理一下整個蒙哥馬利算法的本質，蒙哥馬利算法並不是一個獨立的算法，而是三個相互獨立又相互聯系的算法集合，其中包括

• 蒙哥馬利乘模，是用來計算\(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或2048bit，而我們的計算機不可能存儲完整的大數，因為占空間太大，而且也沒必要。因此，這種基於大數運算的加密體系在計算的時候都是基於\(Z_N\)集合的，自然，蒙哥馬利算法也是基於\(Z_N\)。

在剩余類環上，有兩種重要的運算，一類是簡單運算，也就是加法和減法，另一類復雜運算，也就是乘法。我們比較熟悉的是自然數集上的運算，下面看下怎么從自然數集的運算演變成剩余類環上的運算。

對於加法運算，如果計算\(x\pm y\ (mod\ N)\),\(0\leqslant x,y\)，試想自然數集上的 \(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\)，那么

\[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⋅y)/1000(x⋅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)\)（≡表示取模相等），所以我們可以選擇一個合適的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)\cdot(y\cdot1000)/1000)/1000\) 其實就是蒙哥馬利表示法。

所以我們先定義幾個概念：

• 蒙哥馬利參數
  給定一個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\qquad=547\ (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
由於：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和\(\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}\)，算法如下

至此整個蒙哥馬利算法就全部說完了。通過這個算法，我們可以實現快速冪模。