RSA 的大数运算

俺曾经查阅了网上找得到的各种用于实现RSA 的大数运算库，然而最终还是决定自己动手写一个。因为凡是效率高速度快的代码（crypto+、miracl、freelip、rsaref等），要么使用的数据结构过于复杂，要么编码风格杂乱无章，俺的水平和耐心都实在是有限，以至于无法读懂这些东西。而俺读得懂的一些代码，其实现方式却又过于幼稚，效率极低速度一塌糊涂。俺觉得像俺这样的人不在少数，于是决心写一个清晰易懂，效率也过得去的东西奉献给大家。 这个函数库刚做好的时候，生成1024位的随机密钥耗时大约5 分钟，俺认为是可以接受的。但后来找到一个叫tE! 的老外用miracl库写的RsaTools，发现其生成1024位的密钥耗时不超过三秒钟！于是俺针对俺的代码开始了艰苦的优化工作，希望能达到甚至超过这一水平?一周之后1024位密钥的平均生成时间已经降至5 秒左右，但是单单依靠优化代码来进一步提高速度也非常困难了。于是俺开始借助金山词霸来查阅能够通过google找到的一切与RSA 算法相关的论文，但是网上关于RSA算法的论述绝大多数都是用于硬件实现的，将其算法流程用软件设计语言来实现极其繁琐?而且俺发现这样做下去俺只会离自己的初衷越来越远: 俺的代码将不再清晰易懂?所以俺一度准备放弃? 准备放弃之后，心态平静了许多，再回头去看那些原来不太能够理解的RSA 算法原理，却发现其实也不是那么高深莫测，不急不躁地慢慢看，慢慢想，突然就一下子全明白了。一番改进之后，现在这个版本的函数库同样具有非常简单而清晰的结构，速度也不算慢，生成1024位的密钥在俺PIII 900的笔记本上平均耗时不超过两秒。程序使用C+ 编写，可在VC6.0 下直接编译通过，希望大家喜欢。如果发现Bug 或者有好的修改建议，俺将非常感谢您能够给俺一个Mail。 最后，感谢看雪论坛，感谢雪兄多次热心相助，俺在此学到了很多知识，当然还要乘机拍拍马屁，感谢俺家甜甜的支持！ afanty RSA 原理:选取两个不同的大素数p、q，并计算N=p*q选取小素数d，并计算e，使d*e % (p-1)(q-1)=1对于任意AN：若B=A*d % N则A=B*e % N可见d、e形成了非对称秘钥关系，加密者用公钥d加密，解密者可用私钥e解密，第三者即使拦截了密文B、公钥d和N，在不知道p、q的前提下，无法推算出e，从而无法获得明文A。当N取非常大的值时，将其因式分解成p、q是非常困难的，例如当N为1024 bit时，据分析，需动用价值数千万美金的大型计算机系统并耗费一年的时间?RSA 密钥的选取和加解密过程都非常简洁，在算法上主要要实现四个问题：1 ?如何处理大数运算2、如何求解同余方程 XY % M = 13 ?如何快速进行模幂运算4 ?如何获取大素数实际上，在实现RSA 算法的过程中大家会发现后三个问题不是各自独立的，它们互有关联，环环相套，相信届时你会意识到：RSA算法是一种“优美”的算法！ RSA 依赖大数运算，目前主流RSA 算法都建立在1024位的大数运算之上。而大多数的编译器只能支持到64位的整数运算，即我们在运算中所使用的整数必须小于等于64位，即：0xffffffffffffffff，也就是18446744073709551615，这远远达不到RSA 的需要，于是需要专门建立大数运算库来解决这一问题。 最简单的办法是将大数当作数组进行处理，数组的各元素也就是大数每一位上的数字，通常采用最容易理解的十进制数字09。然后对“数字数组”编写加减乘除函数。但是这样做效率很低，因为二进制为1024位的大数在十进制下也有三百多位，对于任何一种运算，都需要在两个有数百个元素的数组空间上多次重循环，还需要许多额外的空间存放计算的进退位标志及中间结果。另外，对于某些特殊的运算而言，采用二进制会使计算过程大大简化，而这种大数表示方法转化成二进制显然非常麻烦，所以在某些实例中则干脆采用了二进制数组的方法来记录大数，当然这样效率就更低了? 一个有效的改进方法是将大数表示为一个n 进制数组，对于目前的32位系统而言n 可以取值为2 的32次方，即 0x100000000，假如将一个二进制为1024位的大数转化成0x10000000进制，就变成了32位，而每一位的取值范围不再是二进制的01或十进制的09，而是0-0xffffffff，我们正好可以用一个32位的DWORD （如：无符号长整数，unsigned long） 类型来表示该值。所以1024位的大数就变成一个含有32个元素的 DWORD数组，而针对 DWORD数组进行各种运算所需的循环规模至多32次而已。而且0x100000000 进制与二进制，对于计算机来说，几乎是一回事，转换非常容易? 例如大数18446744073709551615，等于 0xffffffff ffffffff，就相当于十进制的99：有两位，每位都是0xffffffff。而18446744073709551616等于0x0000000100000000 00000000，就相当于十进制的100：有三位，第一位是1 ，其它两位都是0 ，如此等等。在实际应用中，“数字数组”的排列顺序采用低位在前高位在后的方式，这样，大数A 就可以方便地用数学表达式来表示其值：A=Sumi=0 to n(Ai*r*i)，r=0x100000000，0=Ai=B：A=Sumi=0 to p(Ai*r*i)B=Sumi=0 to q(Bi*r*i)C=Sumi=0 to n(Ci*r*i)r=0x100000000，0=Ai,Bi,Ci=q则当C=A+B、C=A-B、C=A*B时，我们都可以通过计算出Ci来获得C：C=A+B，显然Ci不总是等于Ai+Bi，因为Ai+Bi可能0xffffffff，而Ci必须=0xffffffff，这时就需要进位，当然在计算Ci-1时也可能产生了进位，所以计算Ci时还要加上上次的进位值。如果用一个64位变量result来记录和，另一个32位变量carry来记录进位，则有： carry = 0 FOR i FROM 0 TO p result=Ai+Bi+carry Ci=result%0x100000000 carry=result/0x100000000 IF carry=0 n=p Else: n = p + 1C=A-B，同理Ci不总是等于Ai-Bi，因为Ai-Bi可能=0，这时就需要借位，同样在计算Ci-1时也可能产生了借位，所以计算Ci时还要减去上次的借位值: carry = 0 FOR i FROM 0 TO p IF Ai-Bi-carry=0 Ci=Ai-Bi-carry carry = 0 Else Ci=0x100000000+Ai-Bi-carry carry = 1 n = p WHILE Cn=0 n=n-1C=A*B，首先我们需要观察日常做乘法所用的“竖式计算”过程： A3 A2 A1 A0* B2 B1 B0-= A3B0 A2B0 A1B0 A0B0+ A3B1 A2B1 A1B1 A0B1+ A3B2 A2B2 A1B2 A0B2-= C5 C4 C3 C2 C1 C0可以归纳出：Ci=Sumj=0 to q(Ai-j*Bj)，其中i-j必须=0且=p。当然这一结论没有考虑进位，虽然计算Ai-j*Bj和Sum的时候都可能发生进位，但显然这两种原因产生的进位可以累加成一个进位值?最终可用如下算法完成乘法: n = p + q - 1 carry = 0 For i FROM 0 TO n Result = carry For j FROM 0 TO q IF 0=i-j=p result=result+Ai-j*Bj Ci=result%0x100000000 carry=result/0x100000000 IF carry!=0 n = n + 1 c n = carry对于C=A/B，由于无法将B 对A“试商”，我们只能转换成Bq对Ap的试商来得到一个近似值，所以无法直接通过计算Ci来获得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） 11 x - 5 y = 1 11%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=1 While a! = 0 i = y y = X - y * a / b X = i i = a a=b%a b = i IF x= 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=Sumi=0 to n(Ei*2*i)，0=Ei=1，则： d = 1 For i = n To 0 D=D*D % N IF Ei=1 D=D*C % N RETURN D 这样，模幂运算就转化成了一系列的模乘运算。 对于乘模运算 A*B%N，如果A、B都是1024位的大数，先计算A*B，再% N，就会产生2048位的中间结果，如果不采用动态内存分配技术就必须将大数定义中的数组空间增加一倍，这样会造成大量的浪费，因为在绝大多数情况下不会用到那额外的一倍空间，而采用动态内存分配技术会使大数存储失去连续性而使运算过程中的循环操作变得非常繁琐。所以模乘运算的首要原则就是要避免直接计算A*B。 设A=Sumi=0 to k(Ai*r*i)，r=0x10000000，0=Air，则： C = A*B = Sumi=0 to n(A*Bi*r*i) %N可以用一个循环来处理: C=0; For i = n To 0 c = c * r C=C+A*Bi RETURN 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 0 C=C*r %N C=C+A*Bi %N RETURN C 这样产生的最大中间结果是A*Bi 或C*r，都不超过1056位，空间代价会小得多，但是时间代价却加大了，因为求模的过程由一次变成了多次。对于孤立的乘模运算而言这种时间换空间的交易还是值得的，但是对于反复循环的乘模运算，这种代价就无法承受，必须另寻出路。 由于RSA 的核心算法是模幂运算，模幂运算又相当于模乘运算的循环，要提高RSA 算法的效率，首要问题在于提高模乘运算的效率。不难发现，模乘过程中复杂度最高的环节是求模运算，因为一次除法实际上包含了多次加法、减法和乘法，如果在算法中能够尽量减少除法甚至避免除法，则算法的效率会大大提高。 设A=Sumi=0 to k(Ai*2*i)，0=Ai=N C=C-N RETURN C 由于在RSA中A、B总是小于N，又0=Ai,C0=1，所以： c = (C+Ai*B+C0*N)/2 c (C+2N)/22 c C+2N c 2N 既然C 总是小于2N，所以求C %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(Ei*2*i)： a =A*R %N X = a FOR i FROM n TO 0 X = X * X * r %N IF Ei=1 X=X*A*R %N X = X * 1 * r %N RETURN X最初: X = A*R %N，开始循环时: X = X * X * r %N = A*R*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 = Sumi=0 to k(Ai*r*i) 0=Ai=N C=C-N RETURN C 因为在循环中每次C =C/r 时，都可能有余数被舍弃。假如我们能够找到一个系数 q，使得(C + Ai*B + q*N) %r =0，并将算法修改为： c =0 FOR i FROM 0 TO k c =C+Ai*B+q*N c =C/r IF C=N C=C-N RETURN C 则C 的最终返回值就是A*B*R %N的精确值，所以关键在于求q。由于： (C + Ai*B + q*N) %r =0= (C %r + Ai*B %r + q*N %r) %r =0= (C0 + Ai*B0 + q*N0) %r =0 若令N0*N0 %r =1，q=(C0+Ai*B0)*(r-N0) %r，则： (C0 + Ai*B0 + q*N0) %r = (C0+Ai*B0 - (C0+Ai*B0)*N0*N0) %r) %r = 0于是我们可以得出r为任何值的蒙哥马利算法: m=r-N0 c =0 FOR i FROM 0 TO k q=(C0+Ai*B0)*m %r c =(C+Ai*B+q*N)/r IF C=N C=C-N RETURN C 如果令 r=0x100000000，则 %r 和 /r 运算都会变得非常