蒙哥马利算法简介

1、蒙哥马利( Montgomery )算法简介俺曾经查阅了网上找得到的各种用于实现 RSA 的大数运算库，然而最终还是决 定自己动手写一个。 因为凡是效率高速度快的代码 ( crypto+ 、 miracl 、 freelip 、 rsaref 等)，要么使用的数据结构过于复杂， 要么编码风格杂 乱无章，俺的水平和 耐心都实在是有限，以至于无法读懂这些东西。而俺读得懂 的一些代码，其实现方 式却又过于幼稚，效率极低速度一塌糊涂。俺觉得像俺这样 的人不在少数，于是决 心写一个清晰易懂，效率也过得去的东西奉献给大家。 这个函数库刚做好的时候，生成 1024 位的随机密钥耗时大 约 5 分钟，俺认为是

2、 可以接受的。但后来找到一个叫 tE! 的老外用 miracl 库写的 RsaTools ，发现其生成 1024 位的密钥耗时不超过三秒钟！于是俺针对俺的代码开 始了艰苦的优化工作，希 望能达到甚至超过这一水平。一周之后 1024 位密钥的平均 生成时间已经降至 5 秒左右，但是单单依靠优化代码来进一步提高速度也非常困难了。 于是俺开始借助金山词霸来查阅能够通过 google 找到的一切与 RSA 算法相关的论文，但是网上关于 RSA 算法的论述绝大多数都是用于硬件实现的，将其算法流程用 软件设计语言来实现极 其繁琐。而且俺发现这样做下去俺只会离自己的初衷越来越 远：俺的代码将不再清晰易懂。所

3、以俺一度准备放弃。 准备放弃之后，心态平静了许多，再回头去看那些原来不太 能够理解的 RSA算 法原理，却发现其实也不是那么高深莫测，不急不躁地慢慢 看，慢慢想，突然就一 下子全明白了。一番改进之后，现在这个版本的函数库同样 具有非常简单而清晰的结构，速度也不算慢，生成 1024 位的密钥在俺 PIII 900 的笔记本上平均耗时不超过两秒。程序使用 C+ 编写，可在 VC6.0 下直接编译通过， 希望大家喜欢。如果发现Bug 或者有好的修改建议，俺将非常感谢您能够给俺一个 Mail 最后，感谢看雪论坛，感谢雪兄多次热心相助，俺在此学到 了很多知识，当然 还要乘机拍拍马屁，感谢俺家甜甜的支持！

4、 afanty原理介绍RSA 原理：选取两个不同的大素数p、q，并计算N=p*q选取小素数d，并计算e,使d*e % (p-1)(q-1)=1 对于任意 A<N ：若 B=A*d % N则 A=B*e % N可见d、e形成了非对称秘钥关系，加密者用公钥d加密，解密者可用私钥 e 解密，第三者即使拦截了密文B、公钥d和N，在不知道p、q的前提下，无法推算出e，从而无法获得明文A。当N取非常大的值时，将其因式分解成p、q 是非常困难的，例如当 N为 1024bit 时，据分析， 需动用价值数千万美金的大型计算机系统并 耗费一年的时间。RSA 密钥的选取和加解密过程都非常简洁，在算法上主要要实

5、现四个问题：1、如何处理大数运算2、如何求解同余方程 XY % M = 13、如何快速进行模幂运算4、如何获取大素数 实际上，在实现 RSA 算法的过程中大家会发现后三个问题不是各自独立的，它们 互有关联，环环相套，相信届时你会意识到： RSA 算法是一种 “优美”的算法！大数存储RSA 依赖大数运算，目前主流 RSA 算法都建立在 1024 位的大数运算之上。而大 多数的编译器只能支持到 64 位的整数运算，即我们在运算 中所使用的整数必须小于等于 64 位，即： 0xffffffffffffffff ，也就是 18446744073709551615 ，这远远达不 到 RSA 的需要，于是

6、需要专门建立大数运算库来解决这一问题。 最简单的办法是将大数当作数组进行处理，数组的各元素也 就是大数每一位上 的数字，通常采用最容易理解的十进制数字 0 9。然后对“数 字数组”编写加减乘 除函数。但是这样做效率很低，因为二进制为 1024 位的大 数在十进制下也有三百多 位，对于任何一种运算，都需要在两个有数百个元素的数组 空间上多次重循环，还 需要许多额外的空间存放计算的进退位标志及中间结果。另 外，对于某些特殊的运 算而言，采用二进制会使计算过程大大简化，而这种大数表 示方法转化成二进制显 然非常麻烦，所以在某些实例中则干脆采用了二进制数组的 方法来记录大数，当然这样效率就更低了。 一

7、个有效的改进方法是将大数表示为一个n 进制数组， 对于目前的 32 位系统而言 n 可以取值为 2 的 32 次方，即 0x100000000 ，假如将一个二进制为 1024 位的大数 转化成 0x10000000 进制，就变成了 32 位，而每一位的取 值范围不再是二进制的 0 1 或十进制的 0 9，而是 0-0xffffffff ，我们正好可以用一个 32 位的 DWORD（如：无符号长整数， unsigned long ） 类型来表示该值。所以 1024位的大数就变成一个含有 32 个元素的 DWORD 数组，而针对DWORD 数组进行各种运算所需的循环规模至多 32 次而已。而且 0

8、x100000000 进制与二进制，对于计算机来说，几乎是一回事，转换 非常容易。例如大数 18446744073709551615 ，等于 0xffffffffffffffff ，就相当于十进制的 99 ：有两位，每位都是 0xffffffff 。而18446744073709551616 等于 0x000000010000000000000000 ，就相当于十进制的 100 ：有三位，第一位是 1 ， 其它两位都是0 ，如此等等。在实际应用中，“数字数组”的排列顺序采用低 位在前高位在后的方式，这样，大数 A 就可以方便地用数学表达式来表示其值： A=Sumi=0 to n(A*r*i)

9、， r=0x100000000 ， 0<=A<r 其中 Sum 表示求和， A 表示用以记录 A 的数组的第 i 个元素， * 表示乘方。 任何整数运算最终都能分解成数字与数字之间的运算，在0x100000000进制下其“数字”最大达至到）xffffffff，其数字与数字之间的运算，结果也必然超出了目前 32 位系统的字长。在 VC+ 中，存在一个 _int64 类型可以处理 64 位的整数，所以 不用担心这一问题，而在其它编译系统中如果不存在 64 位 整形，就需要采用更小的进制方式来存储大数， 例如 16 位的 WORD 类型可以用来表示 0x10000进制。但效率更高的办法还

10、是采用 32 位的 DWORD 类型，只不过将0x100000000进制改成 0x40000000 进制，这样两个数字进行四则运算的最大结果为 0x3fffffff *0x3fffffff ，小于0xffffffffffffff ，可以用一个双精度浮点类型（ double ， 52 位有 效数字）来储存这一中间结果，只是不能简单地用高位低位来将中间结果拆 分成两个“数字”。加减乘除设有大数 A、B 和 C ，其中 A>=B ：A=Sumi=0 to p（A*r*i）B=Sumi=0 to q(B*r*i)C=Sumi=0 ton(C*r*i) r=0x100000000 ， 0<=

11、A,B,C<r ，p>=q 则当 C=A+B 、 C=A-B 、 C=A*B 时，我们都可以通过计算出 C 来获得 C ：C=A+B ，显然 C 不总是等于 A+B ，因为 A+B 可能 >0xffffffff ， 而C必须<=Oxffffffff，这时就需要进位，当然在计算Ci-1时也可能产生了进位，所以计算 C 时还要加上上次的进位值。如果用一个 64 位变量 result 来记录和，另一个 32 位变量 carry 来记录进位，则有：carry=0FOR i FROM 0 TOpresult=A+B+carryC=result%0x100000000 carry=

12、result/0x100000000IF carry=0 n=pELSE n=p+1C=A-B，同理C不总是等于 A-B，因为A-B可能<O ,而C 必须 >=0 ，这时就需要借位，同样在计算 Ci-1 时也可能产生了借位， 所以计算 C 时还要减去上次的借位值：carry=0FOR i FROM 0 TO pIFA-B-carry>=0C=A-B-carrycarry=0ELSEC=0x100000000+A-B-carrycarry=1n=pWHILECn=0 n=n-1C=A*B ，首先我们需要观察日常做乘法所用的“竖式计算”过 程：A3 A2 A1 A0* B2 B1

13、 B0= A3B0 A2B0A1B0 A0B0+ A3B1 A2B1 A1B1 A0B1+ A3B2 A2B2 A1B2A0B2= C5 C4 C3 C2 C1 C0可以归纳出： C=Sumj=0 toq(Ai-j*Bj)，其中 i-j 必须 >=O 且<=p。当然这 一结论没有考虑进位，虽然计算 Ai-j*Bj 和 Sum 的时候都 可能发生进位，但显然这两种原因产生的进位可以累加成一个进位值。最终可用如 下算法完成乘法：n=p+q-1carry=0For i FROM 0 TO nresult=carryFor j FROM 0 TOqIF 0<=i-j<=pres

14、ult=result+Ai-j*BjC=result%0x100000000carry=result/0x100000000IFcarry!=0n=n+1Cn=carry对于 C=A/B ，由于无法将 B对A “试商”，我们只能转换成Bq对Ap的试商来得到 一个近似值，所以无法直接通过计算C来获得C ,只能步逼近 C 。由于：B*(Ap/Bq-1)*0x100000000*(p-q)<C令：X=0，重复 A=A-X*B，X=X+(Ap/Bq-1)*0x100000000*(p-q) ，直到 A<B 则： X=A/B ，且此时的 A=A%B的数注意对于任意大数 A*0x1000000

15、00*k ，都等价于将 A组中的各元素左移 k位，不必计算；同样， A/0x100000000*k 则等价于右移。 欧几里得方程在RSA算法中，往往要在已知 A、N的情况下，求 B,使得 （A*B）%N=1 。即相当于求解 B、 M 都是未知数的二元一次不定方程 A*B-N*M=1 的最小整数解。而针对不定方程 ax-by=c 的最小整数解，古今中外都进行过详尽的研究，西方 有著名的欧几里德算法， 即辗转相除法， 中国有秦九韶的“大 衍求一术”。事实上 二元一次不定方程有整数解的充要条件是 c 为 a、 b 的最大 公约数。即当 c=1 时， a、 b必须互质。而在 RSA 算法里由于公钥 d

16、 为素数，素数与任 何正整数互质，所以可以通过欧几里德方程来求解私钥 e 。 欧几里德算法是一种递归算法，比较容易理解： 例如： 11x-49y=1 ，求 x（a ） 11 x - 49 y = 1 49%11=5 ->（b ） 11 x - 5 y = 111%5 =1 ->（c ） x - 5 y = 1令 y=0 代入（ c ）得 x=1令 x=1 代入（ b ）得 y=2令 y=2代入（ a ）得 x=9 同理可使用递归算法求得任意 ax-by=1 （ a、b 互质）的解 实际上通过分析归 纳将递归算法转换成非递归算法就变成了大衍求一术： x=0,y=1WHILE a!=0

17、 i=y y=x-y*a/b x=i i=a a=b%a b=i IF x<0 x=x+b RETURN x 模幂运算 模幂运算是 RSA 的核心算法，最直接地决定了 RSA 算法的性能。针对快速模幂 运算这一课题，西方现代数学家提出了大量的解决方案，通 常都是先将幂模运算转 化为乘模运算。例如求 D=C*15 % N ，由于： a*b % n = (a % n)*(b % n) %n，所以：C1 =C*C % N =C*2 % NC2 =C1*C % N =C*3 % NC3 =C2*C2 % N =C*6 %NC4 =C3*C % N =C*7 % NC5 =C4*C4 % N =C

18、*14 % NC6 =C5*C % N =C*15 %N即：对于 E=15 的幂模运算可分解为 6 个乘模运算，归纳分 析以上方法可以发现对于任意E，都可采用以下算法计算D=C*E % N :D=1WHILE E>=0IF E%2=0C=C*C % NE=E/2ELSED=D*C %NE=E-1RETURN D继续分析会发现，要知道 E 何时能整除 2，并不需要反复 进行减一或除二的操作，只需验证 E 的二进制各位是 0 还是 1 就可以了，从左至右或从右至左验证都可 以，从左至右会更简洁，设 E=Sumi=0 to n(E*2*i) ， 0<=E<=1 ，则：D=1FOR

19、i=n TO 0D=D*D % NIF E=1 D=D*C % NRETURN D 这样，模幂运算就转化成了一系列的模乘运算。模乘运算对于乘模运算 A*B%N ，如果 A、B 都是 1024 位的大数， 先 计算A*B，再%N ，就会产生 2048 位的中间结果，如果不采用动态内存分配技术就 必须将大数定义中的数组 空间增加一倍，这样会造成大量的浪费，因为在绝大多数情 况下不会用到那额外的 一倍空间，而采用动态内存分配技术会使大数存储失去连续 性而使运算过程中的循 环操作变得非常繁琐。所以模乘运算的首要原则就是要避免 直接计算 A*B 。设 A=Sumi=0 to k(A*r*i) ， r=0

20、x10000000 ， 0<=A<r ， 则：C = A*B = Sumi=0 to n(A*B*r*i) %N 可以用一个循环来处理：C=0;FOR i=n to 0C=C*rC=C+A*BRETURN C 这样将一个多位乘法转换成了一系列单位乘法和加法，由于： a*b %n = (a%n * b%n) %n a+b %n = (a%n + b%n) %n所以， 对于求 C=A*B %N ，我们完全可以在求 C 的过程中完 成：C=0;FOR i=n to 0C=C*r %NC=C+A*B %NRETURN C 这样产生的最大中间结果是 A*B或 C*r ，都不超过 1056 位

21、，空间代价会小得 多，但是时间代价却加大了，因为求模的过程由一次变成了 多次。对于孤立的乘模运算而言这种时间换空间的交易还是值得的，但是对于反复 循环的乘模运算，这种代价就无法承受，必须另寻出路。蒙哥马利模乘由于 RSA 的核心算法是模幂运算，模幂运算又相当于模乘 运算的循环，要提高RSA算法的效率， 首要问题在于提高模乘运算的效率。 不难发现， 模乘过程中复杂度最高的环节是求模运算，因为一次除法实际上包含了多次 加法、减法和乘法，如果在算法中能够尽量减少除法甚至避免除法，则算法的效率 会大大提高。设 A=Sumi=0 to k(A*2*i) ， 0<=A<=1 ，则：C= A*B

22、 = Sumi=0 to k(A*B*2*i)可用循环处理为：C=0FOR i FROM k TO 0C=C+A*BRETURN C若令 C= A*B*2*(-k) ，则：C= Sumi=0 to k(A*B*2*(i-k) 用循环处理即：C=0FOR i FROM 0 TO kC=C+A*BC=C/2RETURN C通过这一算法求 A*B*2*(-k) 是不精确的， 因为在循环中每次 除以 2 都可能有余数被舍弃了，但是可以通过这一算法求 A*B*2*(-k) %N 的精确值，方法是在对 C 除2 之前，让 C 加上 C0*N 。由于在 RSA 中 N 是两个素数的 积，总是奇数，所以当 C

23、 是奇数时， C0=1 ， C+C0*N就是偶数，而当 C为偶数时C0=0 , C+CO*N还是偶数， 这样 C/2 就不会有余数被舍弃。 又因为 C+N %N= C%N ，所以在计算过程中加若干次 N，并不会影响结果的正确性。可以将算法整理如下：C=0FOR i FROM 0 TO kC=C+A*BC=C+C0*NC=C/2IFC>=N C=C-NRETURN C由于在 RSA 中 A、B 总是小于 N，又 O<=A,CO<=1 ， 所以：C = (C+A*B+C0*N)/2C < (C+2N)/22C < C+2NC <2N既然C总是小于2N，所以求C

24、%N就可以很简单地在结束 循环后用一次减法来完成，即在求 A*B*2*(-k)%N 的过程中不用反复求模，达到了我们避免做除法的目的。当然，这一算法求得的是 A*B*2*(-k) %N ，而不是我们最初需要的 A*B%N 。但是利用 A*B*2*(-k) 我们同样可以求得 A*E %N 。设 R=2*k %N ，R=2*(-k) %N ，E=Sumi=0 to n(E*2*i)A=A*R %NX=AFOR i FROM n TO 0X=X*X*R %NIF E=1 X=X*A*R%NX=X*1*R %NRETURN X最初：X = A*R %N ，开始循环时：X = X*X*R %N= A*R

25、*A*R*R %N= A*2*R %N反复循环之后：X = A*E*R %N最后：X = X*1*R %N= A*E*R*R %N= A*E %N 如此，我们最终实现了不含除法的模幂算法，这就是著名的 蒙哥马利算法，而X*Y*R %N 则被称为“蒙哥马利模乘”。以上讨论的是蒙哥马利模乘最简 单，最容 易理解的二进制形式。蒙哥马利算法的核心思想在于将求 A*B%N 转化为不需要反复取模的 A*B*R %N ，但是利用二进制算法求 1024 位的A*B*R%N ，需要循环 1024 次之 多，我么必然希望找到更有效的计算 A*B*R %N 的算法。 考虑将 A 表示为任意的 r 进制：A = Su

26、mi=0 to k(A*r*i) 0<=A<=r 我们需要得到的蒙哥马利乘积为：C= A*B*R %N R=r*(-k) 则以下算法只能得到 C 的近似值C=0FOR i FROM 0 TO kC=C+A*BC=C/rIF C>=NC=C-NRETURN C因为在循环中每次 C=C/r 时，都可能有余数被舍弃。假如 我们能够找到一个系数 q ，使得 (C + A*B + q*N) %r=0，并将算法修改为：C=0FOR i FROM 0 TO kC=C+A*B+q*NC=C/rIF C>=NC=C-NRETURN C则C的最终返回值就是 A*B*R %N的精确值，所以关键在 于求 q 。由于：(C + A*B + q*N) %r =0=> (C %r + A*B %r + q*N %r) %r =0=>(C0 + A*B0 + q*N0) %r =0若令 N0*N0 %r =1 ， q=(C0+A*B0)*(r-N0) %r ，则：(C0 + A*B0 + q*N0) %r= (C0+A*B0 - (C0+A*B0)*N0*N0)%r) %r= 0于