第二十三個知識點：寫一個實作蒙哥馬利算法的C程式

這次部落格我将通過對蒙哥馬利算法的一個實際的實作，來補充我們上周蒙哥馬利算法的理論方面。這個用C語言實作的蒙哥馬利算法，是為一個位數為64的計算機編寫的。模數\(m\)是以能和\(2^{64}-1\)一樣大，\(a\)和\(b\)能和\(m-1\)一樣大。我們采用\(r = 2^{64}\)。在之前的部落格裡，給出的大部分資訊都來自于[1],是以請參考這裡的資訊。

在讀過上次部落格後，你知道我們需要四個步驟。為了我們的目的，我們将這些分為三個階段。

1.The GCD Operation

這個函數用了二進制擴充歐幾裡得算法，找出一個\(r^{-1}\)和\(m^{'}\)使得\(rr^{-1} = 1 + mm^{'}\)。這些整數在後面的算法中需要使用到。算法用\(r^{-1}\)和\(m^{'}\)計算出\(m,m^{'}\),這個部落格的目的不是介紹這個二進制擴充歐幾裡得算法。你想知道的更多可以看連結[1]和[2]。

2.Transform the Multipliers

第二階段是計算兩個值\(abar = ar \mod m\)和\(bbar = br \mod m\).因為\(r = 2^{64}\),這裡隻需要右移64位。就是輸出128位，前64位是\(a,b\)的值，後64位都是0。然後計算\(m\)的模數。這個函數接受64位的\(x\)，同時接受低64位的\(y\)和一個\(m\)的值。之後傳回一個64位的值。

uint64 modul64(uint64 x, uint64 y, uint64 z);

uint64

是這樣定義的：

typedef unsigned long long uint64;

3.Montgomery Multiplication

這個函數定義成接受64位的abar,bbar,m和mprime。然後傳回64位的值。

首先計算\(t = abar*bbar\)。這個得到一個128位的整數。

然後計算\(u = (t + ( tm^{'} \mod r)*m)/r\)。\(t\)是一個128位的整數。這裡就可以計算了。定義如下:

tm = tlo*mprime;
mulul64(tm,m,&tmmhi,%tmmlo);
           

然後計算:

ulo = tlo + tmmlo;
uhi = thi + tmmhi;
if (ulo < tlo) uhi = uhi +1; // test for overflow from ulo and add if necessary to uhi
ov = (uhi < thi) | ((uhi == thi) & (ulo < tlo)); // check for carry
           

最後一步約減到\(m\)。

ulo = uhi;
uhi = 0;
if(ov > 0 || ulo >= m) // test if there was overflow or ulo is higher that
             ulo = ulo – m;
return ulo;
           

4.The Inverse Transformation

最後計算\(a*b \mod m = ur^{-1} \mod m\)

調用之前的函數。

mulul64(p, rinv, &phi, &plo); // performs multiplication and returns two 64 bit values phi and plo
p = modul64(phi, plo, m); // returns value of 128bit input mod m
           

這裡\(p\)就是蒙哥馬利算法的結果了。

[1] http://www.hackersdelight.org/MontgomeryMultiplication.pdf

[2] http://www.ucl.ac.uk/~ucahcjm/combopt/ext_gcd_python_programs.pdf