3. 浮点的圆整是依赖于基数的,而且也是不明显的错误的一个根源。
4. 一个事实是整数和浮点数据类型的尺寸由计算机架构决定,这使语言设计者难以保持它们的语言能够真正可移植。
5. 数字只能以固定的依赖机器的精度来表示。在许多应用程序中这是严重的缺陷,例如在新兴的正在成长的公钥加密算法领域。
6. 基数依赖的现象是不易考虑的。例如很难从传统的整形中访问十进制数的某一位。
这里描述了一组用于直接操作多倍精度有理数的C例程,多倍精度整数是其一个兼容子集。同时也可以进行近似的实数运算。
两个新的数据类型是big和flash。前者用于存储多倍精度整数,后者以\不固定斜杠)的形式用分子和分母来存储小数。两者的格式都是一个固定长度的数字数组,外加一个包含符号和长度信息的独立32位整数。
需要用一种内置类型来储存数字,这个类型定义为mr_small,它可以是int、long甚至double。 struct bigtype {
mr_unsign32 L; mr_small *d; };
typedef struct bigtype *big; typedef struct bigtype *flash; 定义变量:
big x, y[10], z[10][10];
要注意一个big变量仅仅是一个指针,big(或flash)变量实体所需要的内存是从堆(也可以从栈)中分配的。因此使用前必须进行初始化。
分母长度为0意味着实际上分母是1,分子长度为0意味着实际上分母是1。零本身定义为一个首元素是0的数字。
flash中的斜杠并不在一个固定位置,它依赖于分子和分母的相对尺寸。
操作flash时,会将它拆分为分别等于分子和分母的两个big值。big的操作由提取并操作其中的每个mr_small元素组合而成。为避免可能的溢出,每个元素中的值通常限定在一个比完整的变量所能表示的最大值小的范围,例如在16位机上是0-32767(2的15次方 - 1),但基于2的16次方进行编程也是可以的,C语言不会报告整数溢出错误。
系统初始化时用户必须指定一个固定的尺寸(以字或字节为单位)和基数,这两个值将应用于所有big和flash变量。基数的上限与机器字长有关,可以使用不超过上限的任意值作为基数。
当使用一个偏小的基数b时,系统实际上将使用b*n作为基数(这能够更高的效率),n是一个满足使b*n能被单个机器字表示的最大整数。
当使用一个完整宽度的基数时,程序通常运行得最快。若使用这个基数,请在调用mirsys时指定基数为0。需要注意的是对某些时间关键的例程,这种模式在部分常用的编译器/处理器组合中是由大量的汇编代码支持的。例如,对于使用Borland/Turbo C和80x86处理器,可以查看mrarth1.c中的源程序。
将符号和分子、分母的大小信息编码到单个字中是可行的,C语言包含位操作的完整结构。
第五章 实现
big类型的算术运算的实现并没有太大的创造性,使用的算法是Knuth[Knuth81]。但是做了一些速度优化的工作。
在耗时的乘除法内部,典型的实现是对操作数的每一位相乘,加上上一次操作的“进位”,将结果分为一个数字和对下一个操作的“进位”。以10进制乘法为为例,考虑8723536221 * 9:
为了正确地处理数字5所在的列,需要计算5 * 9 = 45,加上前一列的进位3,得到48,将其中的8作为
结果,4作为对下一列的进位。
基本的原始操作是计算(a * b + c) / m的商和余数。在上例中: a = 5, b = 9, c = 3, m = 10
这个操作具有惊人的广泛用途,并且它处算法的最里层循环中,因此它的高效实现是非常必要的。 对这个MAD(Multiply, Add and Divide)操作,高级语言在实现时通常具有三个难点: 1. 慢。
2. 商和余数不能同时作为一次除法的结果,因此计算过程执行两次,一次用来获得商,一次用来获得余数。 3. 即使操作结果是两个单值,但中间结果ab + c可能是双倍长度的。很庆幸地C语言拥有long整形可以用来临时保存这个乘积。
因为这些原因,最好用汇编语言来实现这些关键操作。即使C语言的long类型不是双倍长度(在32位机器上,大部分C编译器的int和long都是32位),在机器语言层也通常能处理一个临时的双倍长度的结果。有关详细资料可以查看mrmuldv.any文件。
MIRACL对big和flash类型使用固定长度的数组,这是为了避免可变长度方式下需要进行内存分配和碎片收集的困难度和时间消耗。这也意味着在对big值进行计算时,所有的结果和中间数者必须小于等于在mirsys中指定的固定尺寸。
在实践中,稳态计算中涉及的数字的尺寸大多相近。除非两个数相乘导致产生一个双倍长度的中间结果,但这也通常会接下来立即在一次除法操作中缩小。这种情况的一个典型的例子是Pollard-Brentl因式分解。
请注意上面提到的问题的另一个微观表现。只是为了拷贝这些偶然出现的中间值,而为每个变量分配所需大小的两倍的存储是很可惜的。因此,MIRACL库提供了一个特殊的Multiply,Add,Divide例程(名为mad)。它已经被证明在内存受限的机器上实现大型程序(像Pomerance-Silverman-Montgomery因式分解)时非常有用。
除了基本的算法运算,以下例程也有提供:
1. 生成和判断质数,使用概率质数判断方法[Knuth81]。
2. 生成随机big和flash,基于subtract-with-borrow生成器[Marsaglia]。但是请注意内部实现的的基本
随机数生成器并不具有秘密安全性(cryptographicly secure,保密性?)。在真实的加密应用程序中,它是不够可靠的。mrstrong.c中提供了一个强保密生成器。
3. 计算幂和根。
4. 实现了普通的和扩展的欧几里得最大公约数(Euclidean GCD)算法。
5. 实现了中国剩余定理[Knuth81],以及计算雅可比符号(Jacobi Symbol)[Reisel]。 6. 超大数的乘法,使用快速傅里叶变换(Fast Fourier Transform)[Pollard71]。
在进行大规模模运算(modular arithmetic)时,一个时间上的关键操作是\模乘\,就是对两个数进行乘法之后再除以一个固定的n,取其余数。一个明显的解决方案是用之前提到的mad例程来完成,但Montgomery提出了一种替代方案,这要求要求先将两个数转换为特殊的n-residue样式,在n-residue样式下可以通过一个完全没有用到除法的特殊例程更快的完成模乘操作,然后再将结果转换为普通格式。请注意n-residue的模加与模减操作与普通的类似,使用相同的例程。 鉴于要求对变量进行n-residue格式转换,Montgomery算法只应考虑用于对同一个模数进行大量模运算的计算场合。这在C++环境中使用更为方便,因为屏蔽了这些细节。
MIRACL库中的许多要进行大量模运算的例程在内部使用了蒙哥马利算法,如高度优化的模幂(modular exponentiation)函数powmod,以及实现GF(P)椭圆曲线运算的函数,可以在参考手册中找到相关信息。
为了达到最快的模运算,同时还必须使用汇编,针对特定的模数或特定尺寸的模数进行优化等方法。有相当数量的技术可供使用:
首先是Comba和KCM方法,其实现分别在mrcomba.c和mrkcm.c文件中。这些文件通过向模板文件mrcomba.tpl和mrkcm.tpl中插入.mcs文件中的宏来创建。这是通过宏展开程序mex自动完成的。在目标机器上编译并运行config.c以自动创建合适的mirdef.h和获取继续操作的建议。同时也请阅读kcmcomba.txt。为使嵌入的程序获得尽可能快的性能,当没有与您的编译器/处理器对应的x.mcs文件时,建议你自己创建一个。
另外有两种在一定程序上属于实验性的技术,分别实现在mr87v.c和mr87f.c文件中,仅适合于80x86系列处理器和Borland C++编译器。
在满足条件时,相应的代码将被自动调用。
请务必注意上述四种技术要求编译器支持内联汇编,而且后两者仅用Borland C++ V4.5编译器在80x86
系列处理器上进行过测试。
第一种方法的思路是对乘法和归约过程(multiplication and reduction process)中隐含的循环进行完全的分解和重组,由[Comba]首先提出,[Scott96]对其进行了修改,见mrcomba.tpl。这种方法必须使用固定宽度的模数并在编译时确定(在mirdef.h中定义MR_COMBA宏为模数尺寸,以字为单位)。对中小尺寸的模数这种方法工作得很好,尤其是在GF(p)椭圆曲线加密算法中。若要获得比这还要快的速度,归约过程可对某些具有特别简单的样式(particularly simple form)的模数进行优化。这可以通过手工向mrcomba.tpl插入相应的代码来完成。针对模数p = 2 ^ 192 - 2 ^ 64 - 1的示例代码在例程comba_redc中给出,要运行此特殊的代码,必须在mirdef.h中定义MR_SPECIAL。
这项技术可以与Karatsuba的快速乘法[Knuth81]联合用于提高对大模数的模乘运算速度。可在mirdef.h中定义MR_KCM以调用此Karatsuba-Comba-Montgomery (KCM)方法。以机器字为单位的模数尺寸限定为必须等于MR_KCM*2n,n可为(合理范围内的)任意正整数,这是使用Karatsuba算法的一个后果。例如若在32位机器上定义MR_KCM为8,则可用的模数尺寸为512,1024,2048位等。
另一个选择是利用浮点协处理器(floating point co-processor),如果存在的话。它的乘法指令通常要比整数单元的要快[Rubin]。Intel Pentium处理器内嵌的协处理器可以用三个周期完成一次乘法,而一次整数乘法需要10个周期。不过对Pentium Pro,II或III情况可能不同。同时协处理器有8个额外的寄存器,可直接操作64位数字。这个特性给程序员提供了一些额外的可以加以利用的灵活性。一些实验性的代码在mr87f.c和mr87v.c中提供,可通过在mirdef.h中定义MR_PENTIUM加以利用。使用config.c来生成mirdef.h——这时底层类型必须选择double。模块mr87v.c实现了简洁的循环代码(compact looping code),可工作于小于某个最大值的任意模数。模块mr87f.c展开了这些循环以获得更快的速度,但也因此更为宠大并且要求固定尺寸的模数。请注意这些运行方式与完全宽度基数(full-width base)是不兼容的,它们通常在基数为2^28或2^29(config.c会为你计算出这个值)的数字上工作得最好。同时也请注意即使这些方法可以加速Pentium上的模幂运算,但它也可能在大部分的其它80x86处理器上更慢,所以请小心使用。在一次测试中,对一个2048位的数字求2048位的幂再对一个2048位的模数取模的操作在60MHz的Pentium上耗时2.4秒。