大整数乘法的实现与分析.doc

大整数乘法的实现与分析.doc

《大整数乘法的实现与分析.doc》由会员分享，可在线阅读，更多相关《大整数乘法的实现与分析.doc（64页珍藏版）》请在冰点文库上搜索。

本科毕业设计（论文）

大整数乘法的实现与分析

课程名称

题目名称

学生学院

专业班别

学号

学生姓名

指导教师

2008年6月

摘要

随着计算机信息安全要求的不断提高,密码学被大量应用到生活中。

RSA、ElGamal、DSA、ECC等公钥密码算法和数字签名算法都建立在大整数运算的基础上，比较耗时的大整数乘法、模乘、幂运算、模幂乘运算等却被上述算法大量使用，它们的运算速度对这些算法的高效实现起着重要的作用，如何快速实现上述几种运算是公钥密码领域普遍关注的热点问题。

本文基于32位的系统，首先采用模块化的思想建立大整数运算库的基础框架，在实现一些辅助函数后在此框架上讨论并实现多精度大整数的基本乘法、Comba乘法、Karatsuba乘法、各种平方算法、Barrett缩减、Mentgomery缩减、模乘、Mentgomery模幂乘等算法及相关的算法。

本文讨论的所用程序均采用C语言编写，所采用的优化也均建立在C语言这一层面上，在保证算法有足够高的效率的同时力求代码清晰易懂，函数接口简单明了，具有可移植性和稳定性。

关键词：

多精度大整数，缩减，模幂乘，滑动窗口

注：

本设计（论文）题目来源于实际应用。

Abstract

Nowadays,ascomputerinformationsecurityrequirementsimprovecontinuously,thecryptologyiswidelyappliedtolife.PublickeycryptographicalgorithmsanddigitalsignaturealgorithmssuchasRSA,ElGamal,DSA,ECCareallbaseonmultipleprecisionarithmetic.Multipleprecisionmultiplication,modularmultiplication,exponentiation,modularexponentiationwhichneedmoreworkingtimeisusedbypublickeycryptographicalgorithmswidely,theirspeedisveryimportanttotheimplementationsofthosealgorithms.Howtofastimplementthosearithmeticaboveisthehottopicinthepublickeycryptographicfield.

Thispaperisbasedonthe32bitsystem.First,wefoundthemodularfoundationofmultipleprecisionarithmeticlibrary;Aftersomeauxiliaryfunctionisformed,wediscussandimplementthemultipleprecisionintegerbasicmultiplication,Combamultiplication,Karatsubamultiplication,kindsofsquarealgorithms,Barrettreduction,Montgomeryreduction,MontgomeryModularExponentiationalgorithmandsomerelationalfunction.AllthealgorithmsdiscussinthispaperisimplemententirelyinportableISOCandtheoptimizationofthosealgorithmsimplementationsisbuiltontheclanguagelevel.Clearcode,simpleapplicationprogramminginterfaceisasimportantastheefficiency,therobustnessandtheportability.

Keywords:

MultiplePrecisionInteger，Reduction，ModularExponentiation，SlidingWindow

目录

1绪论 1

1.1题目背景 1

1.2国内外研究状况 1

1.3本文构成及研究内容 2

2基础设置 3

2.1大整数的表示 3

2.2部分预定义的量 4

2.3底层函数 5

2.3.1初始化大整数 5

2.3.2清除大整数 6

2.3.3扩展大整数的精度 6

2.3.4把一个大整数赋值给另一大整数 6

2.3.5格式化的表示 6

3特殊优化 7

3.1字移位 7

3.2乘以2或除以2 7

3.3乘以或除以2的幂 8

3.4模2的幂 10

4乘法 12

4.1传统的乘法 12

4.2使用Comba思想的乘法 14

4.3只计算低半部或高半部的乘法 17

4.4Karatsuba乘法 17

4.5平方 22

4.5.1基本的平方算法 22

4.5.2使用Comba思想的平方 24

4.5.3Karatsuba平方 26

5模缩减 27

5.1Barrett缩减 27

5.2Montgomery缩减 30

6幂乘 33

6.1幂乘概述 33

6.1.1字幂乘 34

6.2k-ray幂乘 35

6.2.1滑动窗口幂乘 36

6.3模幂乘 36

结论 43

参考文献 44

致谢 46

附录 47

1绪论

1.1题目背景

科技的发展特别是网络的发展使计算机深入到了各行各业的方方面面，计算机在带来方便和提高了工作效率的同时却也带来了各种各样的新问题，其中信息安全问题最为突出，随着计算机信息安全要求的不断提高,计算机保密系统已变得越来越重要，密码学应用不再是局限于军事、国防等有限领域，而是迅速的走进了千家万户，如CA认证、电子政务、电子商务、数字签名、身份认证、密钥分发等。

RSA、ElGamal、DSA、ECC等公钥密码算法和数字签名算法都建立在大整数运算的基础上，耗时的大整数乘法、模乘、幂运算、模幂乘运算等更是被上述算法大量使用，它们的运算速度对这些算法的高效实现起着重要的作用，如何快速实现上述几种运算是公钥密码领域普遍关注的热点问题。

1.2国内外研究状况

长期以来，各方面的工作者对大数模幂的快速实现问题进行了大量研究，主要围绕模幂算法设计与优化、模乘算法设计与优化、专用芯片设计等，并且已经取得不少研究成果。

模幂通常都转化为一系列模乘和模平方运算，目前较好的算法已经能够将1次n比特数的模幂运算转化为约5n/4次n比特的模乘运算，再减少模乘的次数的难度很大，因此，提高模乘的速度是模幂快速实现的关键[1]。

目前，模乘主要有估商型、加法型和Montgomery型3类方法。

1960年，Pope和Stein就提出了采用估商方法将模乘转化为一系列乘法和除法进行计算的思想；1981年，Knuth具体给出了一种转换为乘法和除法的估商型模乘算法[2]；1987年，Barrett提出了一种转换为乘法和加法而没有除法的估商型模乘算法[3]。

1983年，Blakley提出了一种加法型模乘算法，其设计思想是将模乘转换为一系列加法[4]。

为减少加法次数，人们提出了窗口模乘算法和滑动窗口模乘算法，并相继提出了不少改进方法，获得较理想的结果。

1985年Montgomery提出了一种计算模乘的有效算法，其设计思想是将普通模乘转换为易于计算的特殊模乘[5]。

此后，人们提出了不少基于Montgomery思想的改进算法，并得到了广泛的实际应用。

大多数情况下，一种算法的理论描述和实际实现之间有不小的差距，是两个完全的不同的概念，因此，众多学者为这些优秀算法的具体的软、硬件实现、优化付出了辛勤的汗水，在软件实现方面这些算法多数被包含在某些算法库中，其中也有不少是开放源代码的算法函数库，最著名的就是GNU的号称地球上最快的多精度大数库GMP，还有Miracl、openssl、cryptopp、cryptlib、flint等优秀的算法库，而我国的郭先强先生的HugeCalc.dll库也同样出名，虽然它不是开放源码的，但它的速度赶得上GMP甚至在某些方面超越了GMP[6]。

然而，正如TomStDenis所说，不存在一个易懂的大数库[7]！

从目前收集到的信息看，凡是效率高的算法实现，要么结构复杂、层次庞大，要么编码风格奇特；所有速度快的库都使用了汇编，同时大部分都使用了MMX、SSE2、SIMD系列指令作加速，也对多核心CPU进行了优化，使用了多线程等，甚至运行时监测CPU类型而使用相关的特殊优化，最大限度地榨取CPU的性能。

然而，这些很好的优化技术却大大地降低了代码的可读性，极大地提高了理解、学习的门槛。

在学术专著方面，大数算术也备受欢迎，DonaldE.Knuth也用了一整本书——《TheArtofComputerProgrammingVolume2》来从理论的角度讲述了多精度算术，并讲解了高效的算法背后关键的数学概念；《HandbookofAppliedCryptography》[8]也讲述了应用密码学相关的大数算术算法的有效实现；而《KryptographieinCundC++》[9]和《BigNumMath:

ImplementingCryptographicMultiplePrecisionArithmetic》等则从编码学的角度对大数算术进行了多方面的讨论。

1.3本文构成及研究内容

本文基于32位的系统，首先采用模块化的思想建立大整数运算库的基础框架，在实现一些辅助函数后在此框架上讨论并实现多精度大整数的基本乘法、Comba乘法、Karatsuba乘法、各种平方算法、Barrett缩减、Mentgomery缩减、模乘、Mentgomery模幂乘等算法及相关的算法。

本文所讨论的所有算法的实现都不使用汇编、多媒体加速指令、多线程等技术，本文讨论的所用程序均采用C语言编写，所采用的优化也均建立在C语言这一层面上，在保证算法有足够高的效率的同时力求代码清晰易懂，函数接口简单明了，具有可移植性和稳定性。

2基础设置

大整数运算是一些相当复杂的运算，虽然本文要实现的只是乘法类运算，但还是很有必要采用模块化思想按照层次结构来设计及实现一些其它的辅助函数，而不是把它们内嵌在算法函数内，这样既能够避免算法函数的程序代码的过分冗长、使代码清晰易懂、突出算法过程又能够使程序易于测试、排错、更新和复用。

由于本文重点在乘法类算法，下面只介绍一些关键的辅助函数，其它辅助函数要到相关算法使用到时再简略介绍。

2.1大整数的表示

大数的表示方式有很多种，最容易想到也是最直接易懂的跟手工计算方式最接近的方式是通过字符数组来保存大数，数组的每个元素保存大数的一个十进制数字，这种方式操作比较简单，但这种方式需要较多的额外运算，而且效率明显不高，也需要比较多的存储空间；另一种方式是采用链表作为存储结构，这种方式可以适应不同长度的大数，但这种方式的存储效率不高，对本身就需要不少内存空间的大数运算来说是个明显的负担，而且频繁的堆操作和解引用操作势必大量增加开销，而链表的不连续的存储方式也大大的降低了CPU的高速缓存cache的命中率,不利于编译器优化速度；效率比较高的，被采用的比较多的方式是用B进制数组存储大数，即把大数转化B进制表示，数组的每个元素存储这个B进制数的一个B进制的数位，这样既方便计算机处理又提高了内存的利用率，同时缩短了大数的实际表示长度，减少了循环的次数，有利于提高效率。

根据内存管理方式的不同，B进制数组的方式还可以细分为两种方式——多精度和固定精度。

在多精度算法中，我们通过分配所需要的新的内存来表示运算的结果，试图容纳任意大小的整数，增加了对数字大小的兼容性，这对于处理未知大小的数字很有用，但是，它有在计算过程中可能需要执行堆操作，这会带来一定的开销；固定精度算法只有一个有限（固定的）空间来存储一个大数的表示，所以计算过程不需要堆操作，但它需要注意溢出问题对运算结果的准确性所带来的影响。

固定精度非常适合用于输入的大小是事先已知的任务中[10]。

为了更好的适应各种算法的需要及避免过度浪费存储空间，本文采用多精度的方式。

1985年，西安交通大学的王永祥副教授曾经采用过万进制的方法表示大数，并实现了自己的大数算法[11]。

不过本文决定将进制B取为2的某次幂，以方便优化及提供空间利用率。

本文讨论32位系统，由于m位的数乘以n位的数最多将产生m+n位的数，虽然C99标准增加了longlong这一64位的整数类型，而VC、BCB等IDE也有__int64这样的64位整数类型，但32位硬件上毕竟不能直接处理64位整数，那势必需要附加内部操作来完成，为了可以方便硬件操作，取B为半个CPU字长的unsignedshortint，即216进制，这样两个B进制数位的乘积可以用一个unsignedlongint来保存而不超出CPU的字长，而且这样处理对于提取乘积的结果及进位都比较方便。

为了方便操作，用一些额外的信息来描述大整数是一个不错的做法，本文的大整数结构主要参考LibTomMath大数库[12]和MPI大数库[13]。

结构体bigIntMP表示一个多精度的大整数：

typedefstruct /*定义多精度大整数的结构*/

{

longintused; /*记录已经使用的元素个数*/

longintallocated; /*记录已经分配的空间*/

intsign; /*记录大整数的符号*/

wordBi*pWordBis; /*指向大整数的数据的指针*/

}bigIntMP;

分别记录大整数分配到的空间和已经分配到的空间能够有效地管理内存，减少堆操作的次数，减少相关操作带来的性能损失；通过一个指针来指向保存大整数的数据的数组，这样有利于内存的动态调整，这样的存储方式的另一个好处是可以不移动数组的任何元素而交换两个大整数，只需要交换三个内置的整数值和一个指针就可以了。

大整数的数位按低位在前的方式存放，则按从低位到高位的顺序把大整数的数位按下标从小到大顺序存放到数组中去，也就是跟我们手写的方式方向相反，这样既有利于扩展大数的长度也有利于缩减。

2.2部分预定义的量

BASE_BI 0x10000UL /*数的进制的基*/

BIT_PRE_WORD 16UL /*每个单精度数字含有的bit数*/

KARATSUBA_MUL_CUTOFF 128 /*达到这一长度则改用Karatsuba方法做乘法运算*/

KARATSUBA_SQR_CUTOFF 192 /*达到这一长度则改用Karatsuba方式做平方运算*/

/*定义信号标识*/

NON_NEG_BI 0 /*非负数*/

NEG_BI 1 /*负数*/

R_OK_BI 0 /*正常返回*/

E_BAD_ALLOCATE_BI -1 /*内存分配失败*/

E_BAD_ARGUMENT_BI -2 /*参数不正确*/

GREATER_BI 1 /*大于*/

EQUAL_BI 0 /*等于*/

LESS_BI -1 /*小于*/

/*定义类型*/

typedefunsignedshortwordBi; /*单精度的数字*/

typedefunsignedlongdbWordBi; /*双精度的数字*/

对于大整数的符号，本文只区分负数与非负数，若大整数的used分量为0则表示该大整数为0，这要求其它函数在运算到使used分量为0的保证设置该大整数为非负。

用wordBi表示大整数的一位数字，一下简称单字或字，两倍于wordBi大小的则称为双字，三倍于wordBi大小的则称为三字，由于双字用unsignedlong来表示，已经达到了32位系统的字长，所以不能够用一个变量方便地表示一个三字，在需要三字的地方本程序通过模拟的方式达到相关效果，详细见后面介绍。

本文中函数的参数传递方式模仿C语言的表达式赋值的方式，即目的操作数放在左边，例如mulMP（c,a,b）表示c=a*b。

本程序所使用到的数量大部分都在头文件中定义为常量，这样可以仅修改头文件而使程序在不同的系统上获得最好的性能，例如可以通过修改wordBi、dbWordBi等几个定义就可以使程序充分地利用64位系统的运算能力。

2.3底层函数

2.3.1初始化大整数

因为采用了动态分配内存的方式，所以大整数需要一些函数处理堆上的操作。

函数initMP（bigIntMP*bi,longintsize）的目的是初始化一个大整数，为避免多次进行堆操作，本函数分配的内存大小有个起跳值INITIAL_BI，如果size不大于该值，则都会分配该值大小的数位，否则在起跳值的基础上以INCREMENT_BI为增量进行递加直到大于size，然后分配该大小的数位。

成功分配后所有数据位都被置为0，符号标记为非负。

在本实现中，INITIAL_BI在头文件中被定义为48+1即每个大整数的最小长度为48*16+1*16bit，而INCREMENT_BI在头文件中被定义为16即256bit，因为512bit的公钥算法已经不安全，所以本程序从768bit起跳，要多加16bit是因为很多公钥算法的长度都是512的倍数，如果大整数的长度刚好等于公钥算法的长度则很多时候会引起不必要的扩充内存的操作，所以本程序加了16bit的零头。

2.3.2清除大整数

函数clearMP（bigIntMP*bi）用于释放大整数所得的堆内存，并将符号标记为非负，要再次使用该数则必须先重新初始化；在较复杂的函数中，若某一步（例如调用子函数）执行失败，则需要调用clearMP函数释放该一步之前初始化的所有的大整数以免做成内存泄漏。

2.3.3扩展大整数的精度

函数extendMP（bigIntMP*bi,longintsize）以INCREMENT_BI为增量在原来的大小上递加直到大于size，然后分配该大小的数位

2.3.4把一个大整数赋值给另一大整数

函数assignMP（bigIntMP*dst,bigIntMP*src）先判断大整数dst能否容纳大整数src，如果不能则先将大整数dst扩大到能够容纳src再复制相应的数据位；若出现自赋值的情况，则该函数能够检测出来，从而省掉复制数据的开销。

2.3.5格式化的表示

函数readRadixMP（bigIntMP*bi,char*str,intradix）将radix进制的表示的字符串读入并转换为大整数；函数toRadixMP（bigIntMP*bi,char*str,intradix）将大整数转换为radix进制的字符串表示。

3特殊优化

3.1字移位

由于将大整数转换成了B进制数数组，则此时可以将大整数看成是B的多项式，则大整数a可以很方便地看成，那么或者的运算就可以通过将数组的元素向左或向右移动b个B进制数位的方式快速地完成，算法复杂度仅为O（n），速度大大得提高了。

函数lShiftWordsMP（bigIntMP*bi,longintpw）用于将大整数bi左移pw个字即乘以Bpw，而函数rShiftWordsMP（bigIntMP*bi,longintpw）则用于将bi右移pw个字，这相当于除以Bpw。

3.2乘以2或除以2

在所有的基为2的幂的二进制系统中，乘以（除以）2都可以通过左（右）移一个bit来快速实现。

/*位左移函数，左移一个二进制位，即乘以2*/

intlShiftBitMP（bigIntMP*dst,bigIntMP*src）

{

longintoldUsed=0;

longinti=0;

dbWordBicarry=0;

wordBi*pDst=dst->pWordBis;

wordBi*pSrc=src->pWordBis;

if（dst->allocatedused+1） /*确保能够容纳移位后的结果*/

{

intresult=0;

if（（result=extendMP（dst,src->used+1）!

=R_OK_BI））/*扩大大整数*/

returnresult;

pDst=dst->pWordBis;

}

oldUsed=dst->used;/*记录目的操作数原来的已用空间，方便后面处理*/

dst->used=src->used;

pSrc=src->pWordBis;

for（i=0;iused;++i） /*分别左移*/

{

pDst[i]=（wordBi）（carry=（（（dbWordBi）pSrc[i]）<<1）

|（carry>>BIT_PRE_WORD））;

}

if（（carry>>BIT_PRE_WORD）!

=0）

{/*若最高数位左移后溢出，则将溢出的比特存到下一个字*/

pDst[src->used]=（wordBi）（carry>>BIT_PRE_WORD）;

++（dst->used）;

}

for（i=dst->used;i

pDst[i]=0;

dst->sign=src->sign;

returnR_OK_BI;

}

注意carry是一个双字，它能够安全地容纳一个字左移位后的结果，carry右移BIT_PRE_WORD个位就得到上一个字被移出的最高位，当前字左移一位后再与该值按位或就完成了前一个字的最高比特的向前传递。

该算法需要仅需要2n次移位和n次或运算（n为大整数的字数）。

除以2与乘以2的算法基本相同，只是移位方向不一样，这里就不列出相关代码了。

左移函数的运行结果如图3.1：

图3.1左移函数运行结果

3.3乘以或除以2的幂

因为本文将大整数的基B定为2