第二十三个知识点：写一个实现蒙哥马利算法的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