1、大整数乘法的实现与分析I摘 要随着计算机信息安全要求的不断提高，密码学被大量应用到生活中。RSA、 ElGamal、DSA、ECC 等公钥密码算法和数字签名算法都建立在大整数运算的基础上，比较耗时的大整数乘法、除法、模乘、幂运算、幂乘等运算却被上述算法大量使用，它们的运算速度对这些算法的高效实现起着重要的作用，如何快速实现上述几种运算是公钥密码领域普遍关注的热点问题。本文基于 32 位的系统，首先采用模块化的思想建立大整数运算库的基础框架，在实现一些辅助函数后在此框架上讨论并实现多精度大整数的基本加法、减法、乘法、除法、平方算法、缩减、模乘、模幂乘等算法。所用程序均采用 C/C+语言编写，所采
2、用的优化也均建立在 C/C+语言这一层面上，在保证算法有足够高的效率的同时力求代码清晰易懂，函数接口简单明了，具有可移植性和稳定性。关键词：多精度大整数，Comba，Montgomery，二分查找，笔算 注：本设计（论文）题目来源于企业项目。IIAbstractNowadays, as computer information security requirements improve continuously, the cryptology has been widely applied to life. Public key cryptographic algorithms and dig
3、ital signature algorithms such as RSA, ElGamal, DSA, ECC are all base on multiple precision arithmetic. Multiple precision multiplication, Division, modular multiplication ,exponen- tiation, modular exponentiation which need more working time is used by public key cryptographic algorithms widely, th
4、eir speed is very important to the implementations of those algorithms. How to fast implement those arithmetic above is the hot topic in the public key cryptographic field.This paper is based on the 32 bit system. First of all, we found the modular foundation of multiple precision arithmetic library
5、; After some auxiliary function is formed, we discuss and implement the multiple precision integer basic addition, Subtraction,multiplication, , kinds of square algorithms,division,reduction, and some relational function. All the algorithm discuss in this paper is implement entirely in portable ISO
6、C/C+and the optimization of those algorithms implementations is built on the C/C+ language level. the algorithm has high enough to ensure the efficiency of the code at the same time strive to clearly understand, simple interface function with portability and stability.Key words: Multiple Precision I
7、nteger,Comba,Montgomery,Binary search,Written calculationIII目录1 绪论 .11.1 题目的背景 .11.2 国内外研究状况 .11.3 本文研究内容 .22 大整数的结构 32.1 大整数的存取结构 .32.1.1 大整数 结构分析 32.1.2 大整数结构 42.2 预定义的变量 .52.3 大整数基本函数定义 .52.3.1 大整数初始化操作 52.3.2 大整数的销毁操作 62.3.3 大整数的扩展 62.3.4 大整数的输入和输出函数 62.4 大整数的移位函数 .72.4.1 字移位运算 72.4.2 比特移位运算 93
8、大整数加法和减法实现 133.1 符号相同的加法运算 .133.2 符号不相同的加法运算 .164 大整数乘法实现 194.1 笔算乘法 194.2 使用 COMBA 方法的快速乘法 .224.3 平方算法 .244.3.1 笔算平方算法 254.3.2 Comba 思想的平方算法 .27IV5 大整数模缩减实现 305.1 模 2 的幂 .305.2 BARRETT 缩减 .315.3 MONTGOMERY 缩减 336 大整数除法实现 376.1 使用减法替换除法运算 .376.2 模拟笔算除法 .387 大整数幂运算实现 437.1 单数位幂乘 .437.2 KRAY 幂乘 .457.3
9、 滑动窗口幂乘 .45结论 .47参考文献 .48致谢 .49附录 A .5011 绪论1.1 题目的背景科学技术和网络的发展使计算机深入到了各行业的方方面面，计算机在带来方便和提高了工作效率的同时却也带来了各种各样的新问题，其中信息安全问题最为突出，随着计算机信息安全要求的不断提高, 计算机保密系统已变得越来越重要。随着香农定理的发表，信息安全技术得到了迅猛的发展。密码学应用不再是局限于军事、国防等有限领域，而是迅速的走进了千家万户。RSA、ElGamal、DSA、ECC 等公钥密码算法和数字签名等算法得到了快速发展。其实现都是建立在大整数运算的基础上，耗时的大整数乘法、除法、模乘、幂运算、
10、模幂乘等运算更是被上述算法大量使用。而计算机微机的字长限制对信息安全中大整数的操作，带来了巨大的困难。它们的运算速度对这些算法的高效实现起着重要的作用，如何快速实现上述几种运算是公钥密码领域普遍关注的热点问题。1.2 国内外研究状况长期以来，各方面的工作者对大数基本运算的快速实现问题进行了大量研究，主要围绕模幂算法设计与优化、模乘算法设计与优化、专用芯片设计等，并且已经取得不少研究成果。模幂通常都转化为一系列模乘和模平方运算，目前较好的算法已经能够将 1 次 n 比特数的模幂运算转化为约 5n/4 次 n 比特的模乘运算，再减少模乘的次数的难度很大，因此，提高模乘的速度是模幂快速实现的关键 1
11、。目前，模乘主要有估商型、加法型和 Montgomery 型 3 类方法。1960 年，Pope 和 Stein 就提出了采用估商方法将模乘转化为一系列乘法和除法进行计算的思想；1981 年，Knuth 具体给出了一种转换为乘法和除法的估商型模乘算法 2；1987 年，Barrett 提出了一种转换为乘法和加法而没有除法的估商型模乘算法 3。1983 年，Blakley 提出了一种加法型模乘算法，其设计思想是将模乘转换为一系列加法。为减少加法次数，人们提出了窗口模乘算法和滑动窗口模乘算法，并相继提2出不少改进方法，获得较理想的结果。1985 年 Montgomery 提出了一种计算模乘的有效算
12、法，其设计思想是将普通模乘转换为易于计算的特殊模乘 4。此后，人们提出了不少基于 Montgomery 思想的改进算法，并得到了广泛的实际应用。大多数情况下，一种算法的理论描述和实际实现之间有不小的差距，是两个完全的不同的概念，因此，众多学者为这些优秀算法的具体的软、硬件实现、优化付出了辛勤的汗水，在软件实现方面这些算法多数被包含在某些算法库中，其中也有不少是开放源代码的算法函数库，最著名的就是 GNU 的号称地球上最快的多精度大数库GMP，还有 Miracl、openssl、cryptopp 、cryptlib、 flint 等优秀的算法库，而我国的郭先强先生的 HugeCalc.dll 库
13、也同样出名，虽然它不是开放源码的，但它的速度赶得上GMP 甚至在某些方面超越了 GMP。然而，正如 Tom St Denis 所说，不存在一个易懂的大数库！从目前收集到的信息看，凡是效率高的算法实现，要么结构复杂、层次庞大，要么编码风格奇特；所有速度快的库都使用了汇编，同时大部分都使用了 MMX、SSE2 、SIMD 系列指令作加速，也对多核心 CPU 进行了优化，使用了多线程等，甚至运行时监测 CPU 类型而使用相关的特殊优化，最大限度地榨取 CPU 的性能。然而，这些很好的优化技术却大大地降低了代码的可读性，极大地提高了理解、学习的门槛。在学术专著方面，大数算术也备受欢迎，Donald E
14、.Knuth 也用了一整本书 The Art of Computer Programming Volume 2来从理论的角度讲述了多精度算术，并讲解了高效的算法背后关键的数学概念；Handbook of Applied Cryptography 6也讲述了应用密码学相关的大数算术算法的有效实现；而Kryptographie in C und C+ 7和BigNum Math: Implementing Cryptographic Multiple Precision Arithmetic等则从编码学的角度对大数算术进行了多方面的讨论。1.3 本文研究内容本文基于 32 位的系统，首先采用模块化
15、的思想建立大整数运算库的基础框架，在实现一些辅助函数后在此运算库框架上讨论并实现多精度大整数的基本加法、减法、乘法、平方算法、缩减算法、模乘、模幂乘等算法。本文讨论的所用程序均采用3C/C+语言编写，所采用的优化也均建立在 C/C+语言这一层面上，在保证算法有足够高的效率的同时力求代码清晰易懂，函数接口简单明了，具有可移植性和稳定性。2 大整数的结构大整数运算是一些相当复杂的运算，本文要实现的是大整数基本运算，采用模块化思想按照层次结构来设计及实现一些其它的辅助函数，而不是把它们内嵌在算法函数内，这样既能够避免算法函数的程序代码的过分冗长，使代码清晰易懂、突出算法过程又能够使程序易于测试、排错
16、、更新和复用。由于本文是介绍大整数的基本算法，在本文开始之前只介绍一些关键的辅助函数，其它辅助函数是在相关算法实现的时候再简略介绍它的功能。2.1 大整数的存取结构2.1.1 大整数结构分析大数的存储方式主要是有两种链式存储方式和顺序存储方式。一是采用链表作为存储结构，这种方式可以适应不定长度的大数，这种方式的存储空间包括整数表示部分和链表指针部分，其空间利用率不高，而且其存储空间是离散的，所以随机访问效率也不高，而且频繁的堆操作和引用操作势必大量增加开销，不利于编译器优化速度；另外，根据内存管理方式的不同，顺序存储方式可以再分为静态存储方式和动态存储方式。静态存储方式数组的大小是事先已经知道
17、的，适合知道大小的整数运算。而因其数组长度是固定不变的，在运算的时候容易造成溢出。所以其不适合不定长度的整数运算。而动态方式，其可以动态分配内存空间。可以根据整数位数的变化调整大小。其是最常用的顺序存储方式，而且顺序存储方式是连续分配空间，所以其可以实现随机访问，提高速度，因为存储空间就是整数本身，没有其他额为开销，所以空间利用率也较高。由于受到计算机字长的限制制约着大整数的运算，所以必须要解决大整数表示问题，通常是采用 B 进制表示，如 ，其表示为：012(.)bxn4（ ） ，而且系数 （ ）必1210*.*nnnnnxbxbx00xin须是计算机可以表示的常数。在基选择上，最容易想到也是
18、最直接易懂的是用整数数组来保存大数，数组的每个元素保存大数的一个十进制数字，这种方式操作比较简单，但这种方式需要较多的额外运算，而且效率明显不高，也需要比较多的存储空间；效率比较高的，被采用的比较多的方式是用 B 进制数组存储大数，即把大数转化 B 进制表示，数组的每个元素存储这个 B 进制数的一个 B 进制的数位，这样既方便计算机处理又提高了内存的利用率，同时缩短了大数的实际表示长度，减少了循环的次数，有利于提高效率。为了更好的适应各种算法的需要及避免过度浪费存储空间，本文采用多精度的方式。1985 年，西安交通大学的王永祥副教授曾经采用过万进制的方法表示大数，并实现了自己的大数算法 11。
19、根据万进制的原理，本文决定将整数进制 B 取为 2 的某次幂。本文是基于 32 位系统，VC+有_int64 这样的 64 位整数类型，但 32 位硬件上毕竟不能直接处理 64 位整数，那势必需要附加内部操作来完成。为了匹配硬件操作，取 B 进制为半个 CPU 字长的 unsigned short int 型，即 216 进制. 由于 m 位的数乘以 n 位的数最多将产生 m+n 位的数,所以两个 216 进制数位的乘积可以用一个 232 数来保存而不超出CPU 的字长。2.1.2 大整数结构为了模块化编程,程序结构清晰易懂。本文在开始之前先对大整数做一个结构封装，使用一个结构体来表示大整数。
20、标示符 MBigInt 表示一个多精度的大整数结构：typedef struct /* 定义大整数的结构 */long int alloclen; /* 记录数组已经分配总的空间大小 */long int length; /* 记录大整数的长度，即系数的个数 */int sign; /* 记录大整数的符号 */un_short *pBigInt;/* 指向大整数的数据的指针，即指向系数的地址 */ MBigInt;上面大整数封装的结构中，分别记录大整数分配到的空间 alloclen 大小和已经在使5用到的空间 length 大小。使用这种结构能够有效地管理内存，减少堆操作的次数，减少相关操作带
21、来的性能损失；通过一个指针来指向保存大整数的数据的数组，这样有利于内存的动态调整，可以不移动数组的任何元素而交换两个大整数，其只需要交换三个内置的整数值和一个指针就可以了。本文所讨论的大整数的数位是按照低位在前的方式存放，则按从低位到高位的顺序把大整数的数位按下标从小到大顺序存放到数组中去，也就是十进制表示方式相反方向，这样既有利于扩展大数的长度也有利于缩减。2.2 预定义的变量因为在编程的时候很多的变量字符比较长，很容易略写某个字符，造成编程错误，所以本文首先对一些变量进行预定义。typedef unsigned short int un_short; /*16 位数的声明符号 */yped
22、ef unsigned long int un_long; /*32 位数的声明符号 */#define BIT_PRE_WORD 16UL /* 每个单精度数字含有的 bit 数 */* 定义信号标识 */RETURN_OK_BINT 1 /* 正常返回 */RETURN_FAILE_BINT 0 /*错误返回*/INITIAL_BINT 49 /*默认分配的大小*/STEP_BIINT 16 /*起跳值*/对于大整数的符号，本文只区分正数与负数，若大整数的 sign 分量为 1 则表示该大整数是正数，这要求其它函数在运算到使 sign 分量为-1 的保证设置该大整数为负。用 un_shor
23、t 表示大整数的一位数字，一下简称单字或字，两倍于 un_short 大小的则称为双字，三倍于 un_short 大小的则称为三字，由于双字用 un_long 来表示，已经达到了32 位 CPU 的字长。所以在需要三字的地方本程序通过模拟的方式达到相关效果，详细见后面介绍。2.3 大整数基本函数定义因为在大整数的基本运算中，很多的操作都是类似和重复的，所以本文在开始介绍6基本运算之前，先对这些基础函数进行说明。2.3.1 大整数初始化操作函数 initMBInt(MbigInt *bint, long int size = 49)的目的是初始化一个大整数，默认长度是 49 字节。因为采用了动态
24、分配内存的方式，所以大整数需要一些函数处理堆上的操作。为了在整数基本运算中减少多次进行堆操作，本函数分配的内存大小有个初始值 INITIAL_BINT，如果大整数的输入大小 size 小于该值，则都会分配该值大小的数位。否则在起跳值的基础上以步长 STEP_BIINT 为增量进行递加直到大于 size，然后分配该大小的数位。成功分配后所有数据位都被置为 0，符号标记为非负。在本实现中，初始值 INITIAL_BINT 在头文件中被定义为 48+1 即每个大整数的最小长度为48*16+1*16bit，而步长 STEP_BINT 在头文件中被定义为 16 即 256bit，因为 512bit 的公
25、钥算法已经不安全，所以本程序从 768bit 起跳，要多加 16bit 是因为很多公钥算法的长度都是 512 的倍数，如果大整数的长度刚好等于公钥算法的长度则很多时候会引起不必要的扩充内存的操作，所以本程序加了 16bit 的零头。2.3.2 大整数的销毁操作函数 deleteMBInt(MbigInt *bint)用于释放大整数所得的堆内存，并将符号标记为非负，要再次使用该数则必须先重新初始化；在较复杂的函数中，若某一步（例如调用子函数）执行失败，则需要调用 deleteMBInt 函数释放该一步之前初始化的所有的大整数以免做成内存泄漏。2.3.3 大整数的扩展扩展是指在运算操作的时候，结果
26、超出了现在大整数的大小就追加空间，使得大整数能够得到足够空间来存储数据。extendMBInt(MBigInt *bi,long int size)以STEP_BINT 为增量在原来的大小上递加直到大于 size，然后分配该大小的数位保存数据。72.3.4 大整数的输入和输出函数因为我们最为熟悉的数字和使用的最多的也是十进制表示形式，所以本文对大整数的输入和输出都是使用十进制形式。函数 read_radix(MBigInt *bi,char *str)是将十进制表示的字符串读入并转换为大整数结构表示；函数 write_radix(MBigInt *bi,char *str)将大整数转换为十进制
27、的字符串表示。当然可以通过对以上两个函数的修改，形成不同的 B 进制输入和输出方式。2.4 大整数的移位函数2.4.1 字移位运算字移位是指对大整数左移或右移动 16 个比特位，即大整数数组的一个元素的大小。从大整数的结构可以知道， 将大整数转换成了 B 进制数数组，此时可以将大整数看成是 B 的多项式，则大整数 a 可以很方便地看成 ，那么10nnaa或者 的运算就可以通过将数组的元素向左或向右移动 b 个 B 进制数位的*ba/b方式快速地完成，算法复杂度仅为 O(n)，速度大大得提高了。函数Left_shift_word(bigIntMP * dst,long int n)用于将大整数
28、dst 左移 n 个字即乘以 Bn；而函数 Right_shift_word(bigIntMP * dst,long int n)则用于将 dst 右移 n 个字，这相当于除以Bn。1、字左移运算字左移运算就是对整数进行一次乘法运算，乘数是基 B 的 n 次方。可以通过简单的节点左移得到。/左移 n 个字int Left_shift_word(MBigInt *dst, long int n) /为目标数组分配足够多空间if(dst-length + n dst-alloclen) int result = 0;if(result = extendMBInt(dst,dst-length+n)
29、 != RETURN_OK_BINT)8return result;for(int i = dst-length-1; i = 0; i-) /左移 n 个字dst-pBigInt i+ n = dst-pBigInt i ;for(i = 0;i pBigInt i = 0;dst-length += n; /大整数长度的设置return RETURN_OK_BINT;运算结果如下：图 2.1 左移字运算结果2、字右移运算，其实就是对整数进行一次除法运算。除数是基 B 的 n 次方。可以通过简单的节点左移得到。但其仅是计算商，而不保存余数。/右移 n 个字int Right_shift_wo
30、rd(MBigInt *dst, long int n)int len = dst-length; /大整数的长度for(int i = 0;i length; i+)/右移 n 个字9dst-pBigInt i = dst-pBigInt i+n ;for(i; i pBigInt i = 0;dst-length -= n; /重设置大整数长度return RETURN_OK_BINT;运算结果如下：图 2.2 右移字运算结果2.4.2 比特移位运算比特移位是指对大整数进行一个比特的左移或右移。根据大整数的结构定义，大整数的基是 2 的幂。所以大整数乘以或除以 2 都可以通过简单的左移一位
31、或右移一位得到结果。左移一个比特就是对大整数进行乘以 2 操作，右移的时候就是对大整数进行除以 2 操作，其也仅是得到除以 2 的商，而没有保存余数。1、大整数左移一个比特位/左移一位运算10int Left_shift_bit(MBigInt *dst, MBigInt *src)long int len = 0; /用来保存目的操作数的程度long int i = 0;unsigned int carry = 0;un_short *pDst = dst-pBigInt; /指向目的大整数存取数组的指针un_short *pSrc = src-pBigInt; /指向源大整数存取数组的指针
32、if(dst-alloclen length+1 ) / 确保目的大整数能够容纳移位后的结果int result = 0;if(result = extendMBInt(dst,src-length + 1) != RETURN_OK_BINT) /*扩大大整数空间，保证能够存储*/return result;pDst = dst-pBigInt;/记录目的操作数原来的已用空间，方便后面处理len = dst-length; dst-length = src-length;pSrc = src-pBigInt;/源大整数数组元素分别左移一位，并赋给目的大整数for(i = 0; i lengt
33、h; +i)pDsti = (un_short)(carry = (pSrci) 16);if(carry 16) != 0) /*若最高数位左移后溢出，则将溢出的比特存到下一个字*/11pDstsrc-length = (un_short)(carry 16);+(dst-length);for(i = dst-length; i sign = src-sign; /符号更改return RETURN_OK_BINT;运算结果如下：图 2.3 左移一比特运算结果1、大整数右移一个比特位/右移一位运算int Right_shift_bit(MBigInt *dst,MBigInt *src)l
34、ong int oldlen = 0; /保存目的大整数的原来长度long int i = 0;unsigned int carry = 0;un_short temp = 0; un_short *pDst = dst-pBigInt; /目的大整数数组指针un_short *pSrc = src-pBigInt; /源大整数数组指针if(dst-alloclen length) /确保目的大整数能够容纳移位后的结果12int result = 0;/扩展目的大整数长度if(result = extendMBInt(dst,src-length) != RETURN_OK_BINT) ret
35、urn result; pDst = dst-pBigInt;oldlen = dst-length;dst-length = src-length;pSrc = src-pBigInt;for(i = src-length -1; i = 0; -i) /* 分别右移一位 */temp = (un_short)(pSrci 1) | (un_short)(carry length; i sign = src-sign; /符号赋值oldlen=dst-alloclen; pDst=dst-pBigInt+dst-alloclen-1;while(0=*pDst-) /重新计算出目的大整数的长
36、度值oldlen-;dst-length=oldlen;return RETURN_OK_BINT;运算结果如下：13图 2.4 右移一比特运算结果3 大整数加法和减法实现由于大整数有正负之分，所以两个大整数相加有四种情况：a+b，a+（-b） ， （-a ）+b 和-a+ （-b ） 。对于 a+b 和-a+(-b)，可以通过对数组对应位元素相加即可，主要是考虑进位问题。而 a+（-b ）和（-a）+b ，其实就是两个大整数的减法，其主要考虑的问题是向高位借位的问题。所以减法可以通过加法的运算得到结果。3.1 符号相同的加法运算符号相同的加法运算的实现过程是基于这样的一个原因：因为两个整数的
37、符号相同，可以直接对两个整数数组对应位元素相加，并考虑进位问题。因为两个整数的数位存在以下三种情况（不妨设两个大整数分别为 src1 和 src2）：src1-length = src2-length；src1-length src2-length； src1-length length。所以在处理加法的时候，对数位较大的整数一分为二，第一部分是跟另一大整数的长度相等。可以先把相同长度部分先计算，然后再对多出的长度部分直接赋值，这样可以不用考虑谁大谁小问题，而且可以加快运算速度。如下所示：表 3.1 分段运算过程1 2 3 4 5 6 7 8 9 4 5 6 7 8 914 9 1 3 5 7
38、 81 2 3 9 1 3 5 7 8 /符号相同的加法运算实现int addMBInt1(MBigInt* dst, MBigInt* src1, MBigInt* src2)long int len = 0; /两个大整数的长度最大值long int dstlen = 0; /目标数组分配的最大值un_short mark = 0;/进位标志unsigned int result;/数组对应元素相加结果long sign = src1-sign;/整数的符号int i = 0;len = (src1-length = src2-length)?src1-length:src2-length
39、;dstlen = len;if(-1 = extendMBInt(dst,len)/扩展足够内存return 0; /较小数位的整数的长度len = (src1-lengthsrc2-length)? src2-length: src1-length;for(i = 0; i pBigInti + src2-pBigInti + mark;dst-pBigInti = (resultmark = result 16;/对第一个大整数数位多出部分进行计算if(src1-length src2-length) while(i length)15 result = src1-pBigInti +
40、mark;dst-pBigInti = result mark = result 16;i+;if(mark != 0)dst-pBigInti = mark;dstlen = src1-length + 1;/对第二个大整数数位多出部分进行计算else if(src1-length length)while(i length)result = src2-pBigInti + mark;dst-pBigInti = result mark = result 16;i+;if(mark != 0)dst-pBigInti = mark;dstlen = src2-length + 1;dst-s
41、ign = sign; /符号的更改16dst-length = dstlen; / 长度设置return 1;运行结果如下：图 3.1 符号相同的加法运算结果3.2 符号不相同的加法运算符号不同的加法运算的过程，也是跟符号相同的加法过程类似。也是对数位较大的整数进行分段，不过，它也有不同的地方。因为符号不相同，这个加法就是减法。因为不知道是哪个大整数较大，当它们做减法运算的时候，就有可能会在最高数位出现负数问题，而大整数各位元素都为正数，所以要多处理一步，即对负数转换成正整。而如果把无符号整数的较大者作为被减数，就可以省略这一步。所以此运算过程是：先对两个大整数进行无符号比较，最大者作为被减
42、数，然后进行减法操作。/符号不相同的加法运算（减法运算）int addMBInt2(MBigInt* dst, MBigInt* src1, MBigInt* src2)long int len = 0; /两个大整数的长度最大值un_short len1 = 0;long int dstlen = dst-alloclen; /目标数组分配的最大值un_short mark = 0;/进位标志17int result = 0;/数组对应元素相加结果int sign = src1-sign;/整数的符号int I = 0;len = (src1-length src2-length)? src
43、1-length: src2-length;len1 = (src1-length src2-length)? src2-length: src1-length;un_short *psrc1=src1 - pBigInt,*psrc2 = src2-pBigInt,*pdst = NULL;if(-1 = extendMBInt(dst,len)return 0;/扩展足够内存int re = compareMBInt(src1,src2); /对两个大整数进行无符号比较if(re = 0) /两个大整数大小相等，进行加法运算后位 0dst - length = 0;dst - sign =
44、 0;return 1;if(re = -1)/无符号比较，大整数 src2 比大整数 src1 大sign = src2-sign;psrc1 = src2-pBigInt;psrc2 = src1-pBigInt;len1 = src1-length;for(i = 0;ipBigInti = result;while(i = 0)dst-pBigInti + = result;mark = 0;else /表示向高一位借了 1dst-pBigInt i + = result + 65536; dst-length = len; /长度设置dst-sign = sign; /符号位设置tr
45、imMBInt(dst); /对结果去掉前面的 0return 1;运算结果如下：19图 3.2 减法运算结果4 大整数乘法实现现在流行的公钥算法很多都是以模幂为基础的，而计算模幂的过程中大量地使用了乘法，所以乘法非常重要。然而，通用的乘法都需要 O( )次的基本运算，直到2n1962 年发现了 Karatsuba 乘法人们才有了历史的突破，同时该技术也促使人们发现了更多的有效算法，如基于傅立叶变换的方法。4.1 笔算乘法对于 m 位和 n 位的输入，传统的乘法需要 m*n 次基本的乘法，也即算法复杂度为O( )。我们用纸和笔做乘法运算时，用乘数的每一位乘以被乘数的每一位并加上上一2n列的进位
46、而产生一行适当移位的中间结果，然后再将各行中间结果相加即得到乘法的最终结果，例如 10 进制下计算 456*32 的过程如表 4.1：表 4.1 笔算乘法的运算过程4 5 6* 3 29 1 2201 3 6 81 4 5 9 2本算法按照上述过程进行计算，但在计算机上最好是把内部的乘法和加法并行的执行，即计算每一行中间结果的同时将该行加到最终结果上，这样既可以省去不少步骤，也避免了中间结果的空间开销和内存管理操作，这个简单的优化能够一定程度上提高效率。本算法的伪代码如表 4.2：表 4.2 笔算乘法算法输入：分别含有 n 和 t 个 B 进制数位的大整数 x、y输出：B 进制的乘积 10nt
47、w1. i 从 0 到 n+t-1,置 Wi = 02. i 从 0 到 t2.1 c 02.2 j 从 0 到 n计算 ijjiuvwxyc置 ,ijc2.3 1in3. 返回 10tw由于两个 16 位的乘积将超过 16 位所能表示的最大数，我们知道一个 32 位的数能够完整地保存两个 16 位数的乘积，上面的伪码中用一个 32 位来保存两个 16 位的数乘积再加上两个单字后的结果，下面讨论一下它的安全性（是否会产生溢出）：上面表达式计算的最大的结果为 ，化简后为 ，正好1616161622321是一个 32 位数的最大值，所以一个 32 位数能够保存该表达式的结果而不产生溢出，程序不会损
48、失精度。/笔算乘法实现代码int mulBasicMBInt(MBigInt *product, MBigInt *biA, MBigInt *biB)MBigInt biT; /计算结果先保存在临时变量中21unsigned int carry = 0; /这个双字用于保存中间结果register un_short *pWordA = biA-pBigInt;register un_short *pWordB = biB-pBigInt;register un_short *pPrd = NULL;int result = 0;long int i = 0;long int j = 0;/初始化临时大整数变量if(result = initMBInt(biT.length = biA-length + biB-length;biT.sign = (biA-sign = biB-sign) ? 1: -1;/设置符号pPrd = biT.pBigInt;for(i = 0; i length; +i) / 按传统的纸笔乘法的方式计算乘积carry = 0;for(j = 0; j length
