学习模幂运算－平方和乘算法
姓 名：李锦 学 院：计算机学院 班 级：计算机二班 学 号：07095203 制 作 日 期：2010.5.2
随着各种加解密算法密钥长度的逐步增加，在一些具有安全性需求的芯片设计中，大规格数据运算的硬件实现已成为硬件设计的主要考虑因素和设计难点．比如RSA等基于大数分解的公钥密码算法，虽然目前密钥长度已达1024位，但是仍然不能避免将被破解的厄运，致使密钥还需进一步增加．这种运算规格的增长不仅使加解密运算速度降低，而且增加了硬件实现的难度．
1 CSTU 安全芯片体系结构简介
随着人们对安全需求的不断增加，采用固定或单一加解密算法的产品已经无法满足人们的需求，目前的安全产品需要经常更换加解密算法甚至改变整个安全策略．适应这种需求常用的方法是在基本运算器之上，使用软件编程的方式灵活的实现算法的转换．但是面对不断升级的软件破解技术的挑战，以及软件方式的低速率性，各种加解密算法也由软件实现向硬件电路实现过渡．为解决这一矛盾．可支持多种加解密算法的硬件安全产品就应运而生，其中基于可重组方式设计的安全芯片无疑又具有领先优势．
CSTU保密终端安全芯片采用了可重组设计思想，综合分析了当前大量使用的DES，AES，IDEA，RSA，MD5等十余种加解密算法的实现过程，支持对称、公钥、摘要密码算法及用户隐秘算法，提供这些算法实现所需的IP平台，不同的用户可以根据自己的需要在平台上进行二次开发，形成自己定义的安全算法及策略．
CSTU安全芯片可用于保密电话、安全卡证或移动安全终端等产品中，这些产品的共同特点是对规模要求比较严格，对公钥密码算法的速度要求不高．为提供对公钥密码算法和数字签名算法的支持，大数运算器成为CSTU安全体系中关键的核心IP．根据实际需求，本设计在满足硬件规模尽可能小同时支持尽可能多的运算功能和多种规格的数据运算的条件下，最终保证整个系统的灵活性．
2 模幂运算
模幂运算是RSA 的核心算法，最直接地决定了RSA 算法的性能。针对快速模幂运算这一课题，西方现代数学家提出了大量的解决方案，通常都是先将幂模运算转化为乘模运算。
例如求D=C**15 % N，由于：a*b % n = (a % n)*(b % n) % n，所以：
C1 =C*C % N =C**2 % N C2 =C1*C % N =C**3 % N C3 =C2*C2 % N =C**6 % N C4 =C3*C % N =C**7 % N C5 =C4*C4 % N =C**14 % N C6 =C5*C % N =C**15 % N
即：对于E=15的幂模运算可分解为6 个乘模运算，归纳分析以上方法可以发现对于任意E，都可采用以下算法计算D=C**E % N：
D=1 WHILE E>=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=Sum[i=0 to n](E*2**i)，0
D=1 FOR i=n TO 0 D=D*D % N IF E=1 D=D*C % N RETURN D
这样，模幂运算就转化成了一系列的模乘运算。 3.模乘运算
对于乘模运算 A*B%N，如果A、B都是1024位的大数，先计算A*B，再% N，就会产生2048位的中间结果，如果不采用动态内存分配技术就必须将大数定义中的数组空间增加一倍，这样会造成大量的浪费，因为在绝大多数情况下不会用到那额外的一倍空间，而采用动态内存分配技术会使大数存储失去连续性而使运算过程中的循环操作变得非常繁琐。所以模乘运算的首要原则就是要避免直接计算A*B。
设A=Sum[i=0 to k](A*r**i)，r=0x10000000，0
这样产生的最大中间结果是A*B 或C*r，都不超过1056位，空间代价会小得多，但是时间代价却加大了，因为求模的过程由一次变成了多次。对于孤立的乘模运算而言这种时间换空间的交易还是值得的，但是对于反复循环的乘模运算，这种代价就无法承受，必须另寻出路。 蒙哥马利模乘
由于RSA 的核心算法是模幂运算，模幂运算又相当于模乘运算的循环，要提高RSA 算法的效率，首要问题在于提高模乘运算的效率。不难发现，模乘过程中复杂度最高的环节是求模运算，因为一次除法实际上包含了多次加法、减法和乘法，如果在算法中能够尽量减少除法甚至避免除法，则算法的效率会大大提高。
设A=Sum[i=0 to k](A*2**i)，0
通过这一算法求A*B*2**(-k)是不精确的，因为在循环中每次除以2都可能有余数被舍弃了，但是可以通过这一算法求A*B*2**(-k) %N的精确值，方法是在对C\'除2之前，让C\'加上C\'[0]*N。由于在RSA中N是两个素数的积，总是奇数，所以当C\'是奇数时，C\'[0]=1，C\'+C\'[0]*N 就是偶数，而当C\'为偶数时C\'[0]=0，C\'+C\'[0]*N还是偶数，这样C\'/2 就不会有余数被舍弃。又因为C\'+N %N = C\' %N，所以在计算过程中加若干次N，并不会影响结果的正确性。可以将算法整理如下： C\'=0 FOR i FROM 0 TO k C\'=C\'+A*B C\'=C\'+C\'[0]*N C\'=C\'/2 IF C\'>=N C\'=C\'-N RETURN C\' 由于在RSA中A、B总是小于N，又0
既然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=Sum[i=0 to n](E*2**i)： A\'=A*R %N X=A\' FOR i FROM n TO 0 X=X*X*R\' %N IF E=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 = Sum[i=0 to k](A*r**i) 0=N C\'=C\'-N RETURN C\'
因为在循环中每次C\'=C\'/r 时，都可能有余数被舍弃。假如我们能够找到一个系数 q，使得(C\' + A*B + q*N) %r =0，并将算法修改为： C\'=0 FOR i FROM 0 TO k C\'=C\'+A*B+q*N C\'=C\'/r IF C\'>=N C\'=C\'-N RETURN C\' 则C\'的最终返回值就是A*B*R\' %N的精确值，所以关键在于求q。由于： (C\' + A*B + q*N) %r =0 ==> (C\' %r + A*B %r + q*N %r) %r =0 ==> (C\'[0] + A*B[0] + q*N[0]) %r =0 若令N[0]*N[0]\' %r =1，q=(C\'[0]+A*B[0])*(r-N[0]\') %r，则： (C\'[0] + A*B[0] + q*N[0]) %r = (C\'[0]+A*B[0] - (C\'[0]+A*B[0])*N[0]\'*N[0]) %r) %r = 0 于是我们可以得出r为任何值的蒙哥马利算法： m=r-N[0]\' C\'=0 FOR i FROM 0 TO k q=(C\'[0]+A*B[0])*m %r C\'=(C\'+A*B+q*N)/r IF C\'>=N C\'=C\'-N RETURN C\'
如果令 r=0x100000000，则 %r 和 /r 运算都会变得非常容易，在1024位的运算中，循环次数k 不大于32，整个运算过程中最大的中间变量C\'=(C\'+A*B+q*N)
3 利用平方-乘算法进行大数计算之C语言代码
在RSA算法中，往往要计算ma mod r的值（很抱歉这里用的符号可能与教科书上并不一致），由于ma的值可能会很大而产生溢出从而导致错误。以前我曾经写过一个逐步取余数的小程序来解决这个问题，回头看来，感觉效率不是很高。下面再给一个利用平方-乘算法解决这个问题的小程序，仅供参考。
程序要求：输入m,a,b ma mod b；输出ma mod b ，a的二进制位数为128位，m和b的二进制位数为8位。
实验代码： #include #include int pfcheng(int m,int a,int b) { int d[100],length=0; int c=1; do { d[length++]=a%2; a=a/2; } while(a!=0); while(length>=0) { c=(c*c)%b; if(d[length]==1) { c=(c*m)%b; } length--; } return c; }
void main() {
int m,a,b; cout>m>>a>>b; cout
