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