技术宅的结界

 找回密码
 立即注册→加入我们

QQ登录

只需一步,快速开始

搜索
热搜: 下载 VB C 实现 编写
查看: 838|回复: 2
收起左侧

【C】大数计算GILIB库

[复制链接]

85

主题

263

帖子

3415

积分

用户组: 管理员

No. 418

UID
418
精华
13
威望
52 点
宅币
1969 个
贡献
1014 次
宅之契约
0 份
在线时间
252 小时
注册时间
2014-8-9
发表于 2016-8-17 17:52:01 | 显示全部楼层 |阅读模式

欢迎访问技术宅的结界,请注册或者登录吧。

您需要 登录 才可以下载或查看,没有帐号?立即注册→加入我们

x
gilib_logo.bmp

232636roy52q92559zy1f5.png

        大整数运算在现代密码学与高精度计算领域有着非常重要的作用。
        在RSA算法中,算出密钥必须对一个大整数进行素因数分解。在高精度计算方面,有时候我们会要求精确到小数点后几百位。这些时候,现有编译系统所能支持的内置最大数据类型long long都不足以满足要求。而我们恰恰希望计算一个几百位,甚至上千位的大整数时候效率还有保证。所以,一个大整数计算的工具库是必须的。
市面上又有多少我们可以用的库呢?
        付费使用的数学软件几乎都有此功能,像是Mathematica,Maple,Matlab。一些免费开源的数学软件:SAGA。还有一些专门的大数计算库:GNU dc、GNU gmp,GAP4。
        GILIB库开发的目的在简单易懂,体积小,容易扩展,纯C写成(有C++接口)可跨平台。
        我主要在此帖介绍一下GILIB整个库的设计与开发GILIB时遇到的问题与总结出的经验。
        GILIB的核心写在 GILIB_KRNL.C 中。核心文件实现了大整数结构(GINT结构)的管理与加减乘除四则基本运算。
        每个大整数分别存储在名为 SN_GINT 的结构体中:
typedef struct SN_GINT
{
        unsigned int   incl;
        unsigned int   icap;
                 int   flag;
        unsigned int * data;
}
GINT, * PGINT;
        其中flag代表了一个GINT所占的实际“节”数。(一节等于一个int的大小,那么被32位编译器编译后一个int占32位,x64上编译,一节占64bit)flag在32位平台上最大可以表示 2147483647 节,flag的最高位用于标示GINT的符号。data是一个指向无符号整数的指针,其地址存储了GINT实际的值。那么如果一个GINT的flag等于-1,data所指向的内存地址中存放着0x0000000A,那么这个GINT表示了-10这个十进制整数。icap是一个GINT实际所占用的节的数量。为了在计算时加速,我们不能每次进位的时候使用realloc重新声请内存空间。通常我们位一个GINT一次声请足够的空间,这个空间大小就是icap所示的大小。所以flag的绝对值表示了一个GINT大整数的有效位数,这些有效位数存储在data所指向的内存区域后。icap表示了该GINT最多可以存储多少节而不需要重新声请内存空间。最后一个incl表示了每一次我需要声请data的空间时,我可以再多声请incl个节。如果一个GINT原来的icap为10,也就是说在32位平台上可以存放10*32bit而不用重新开地址。此时我需要重新relloc为此GINT再加上16个节,同时我将incl设置为4。那么如果声请成功此GINT新的内存空间icap变成了10(原来的)+16+4(增长量incl)=30.我们可以手动控制incl实现灵活的内存空间增长控制。
        整数0在GINT中这样表示:flag=1 且 *data=0.如果一个GINT的flag为0,那么此整数不能被计算(NaN)。
        四则运算是所有高级别运算的基础,因此,级别越低的运算被调用的次数越多,低级别运算的函数实现时也就更要注意效率。
        加法运算基本可以理解为笔算算法。这里需要注意的是为了GINT进位时扩展内存的方便,我们将整数“倒置”存储于GINT结构体data指针中。
        计算456+789的原理如下:
654
987+
5421
        实际执行过程中,一节最大可以容纳 4294967295,我们也可以理解为:在32位平台上GINT是一个 “4294967295进制” 的整数。每逢4294967295,GINT就会进位1位。
这是加法函数giadd在处理两个正整数相加时的情况。如果加数或者被加数内存在负数,giadd函数先判断两数绝对值的大小,而后利用giadd自身或者减法函数gisub进行计算。比如(-7)+3,giadd会计算7-3=4,然后将结果取相反数。
        减法函数gisub原理与giadd类似,故不多做介绍。
        对于乘法运算,我们可以用加法运算累加。
9*5=5+5+5+5+5+5+5+5+5
但是这样的做法效率奇低。如果我们换一种思维来看乘法运算。9*5=5*9=101b*1001b将位数较短的整数作为比较位,5只有三位二进制整数,所以此时5作为比较位。从5的低位依次向5的高位扫描。
5的第一位是1,所以两数之积等于9左移0位(左移0位就是不左移,将乘数9原数加到位0的积中)
5的第二位是0,因此两数之积等于0移动1位。(数字0移位等于不做任何操作,因此这一步判断过后可以很快抛弃)
5的第三位是1,将乘数9左移2位后加到积中。因为积已经等于9,9左移2位等于9*2^2=36,所以积成为9+36=45
计算结束,所以9*5的积为45.
        实际运算过程中,因为0xFFFFFFFFFFFFFFFF这样每一位都是1的被乘数较少,而当有一个二进制位为0,我们就不需要做移位与加法两项操作,因此相比较线性的累加乘法操作大大减少了运算时间。
        解释除法运算的实现之前,我们先说一说比较运算gicompare. 我们不需要逐位比较两个GINT,这样做效率太低了。加速比较算法有三步步:第一步,正数一定大于负数。第二步,谁的有效位多,谁的绝对值肯定大。如果GINT a的flag为-1,而GINT b的flag为-3.首先b的绝对值大于a,其次因为a和b同为负数,则a<b. 第三步我们才需要逐节比较两个整数,但是我们从一个整数的最高位开始比较,如果当前位数能得出大小则立即退出比较循环输出比较结果。111159999999 肯定小于111168888888 此时我们只需要比较5个节(循环5次)。在最坏的情况下 +2222222222222222 与 +2222222222222223 相比,才需要16次循环。而且要记住一节的大小是32或者64位整数,循环一次能比较1节,一个GINT最大不能超过 2147483647 或者 在64位平台上的 9223372036854775807 个节。gicompare的效率因此能够保证。
        除法运算,其实就是笔算方法,而且试商用的是二分法,原理十分简单。
        光有加减乘除四则运算是不行的。我们要实现将GINT转为字符串,和将合法字符串转为GINT的函数。GILIB_KRNL.C 文件中,gitosz函数用来将内存中存储的GINT转为零结尾ASCII字符串。至于sztogi函数的作用正好相反。对于八进制与十六进制的转换相对简单。其中八进制要稍微复杂一点,因为需要将节数按照3的倍数取余下的二进制位。而转为十进制是非常难的算法。其难点在于保证效率的同时,不能增加过大的内存负载。事实上,在GILIB库中,转换十进制的算法是效率最低的。2000的阶乘一共5000多位十进制整数,计算出来花费了63毫秒左右。但是光是从GINT在内存中的二进制形式转换为10进制,就要耗费将近600多毫秒。为了实现时间与空间之间的平衡,我们需要先将GINT转为BCD code(binary code decimal)再将BCD做进一步转码,生成字符串。这个算法,在犹他州大学校园的网站上有公开参考资料:http://www.eng.utah.edu/~nmcdona ... /BCDConversion.html
那么二进制GINT到BCD可以使用这个算法,至于BCD到GINT将这个算法倒过来即可。
        在核心文件 GILIB_KRNL.C 中,剩下的一些函数,我介绍一下我们的“神父”。因为一些功能比较零散所以我将它们放在了“神父”gipriest函数中。作为“大祭司”而不是上帝gipriest无权初始化或者销毁一个GINT,但是他可以使一个GINT变的富有(增加icap)或者贫穷(缩减icap)。还可以判断一个GINT成为“字符串”后所占的字符位数。gipriest的具体使用方法在C文件的函数实现出有详细说明。
        剩下的函数是一些常规整型数据类型转为GINT的函数。还有接生婆giinit与死神gidestory.
        因为种种原因,我没有来得急将GILIB对接CUDA或者opencl实现GPU平行加速计算。但是我在 GILIB_MULTHD.C 文件中实现了一个多线程加速乘法运算的例子。乘法多线程加速运算的原理如下:
        如果我们要计算9*5=45,那么我们还是交换乘数与被乘数9*5=5*9。5*9=2*9+2*9+1*9.这时候我们开两个线程分别计算2*9.每个线程运算完毕后将结果加起来,最后再加上剩下的1个9.多线程并行计算,大大减少了单位时间内单次乘法运算所比较的次数。但是不要迷信多线程。开每个线程所耗费的空间,时间也是非常可观的。多线程并行计算的数据同步问题也是非常难以解决的。我们在编写每个多线程算法之前,还是不要嫌麻烦,画好相容矩阵,做好逻辑分析,尽量避免数据不一致引发的各种惨案。实际测试中,在积为4096节以内的两GINT之间做乘法运算,一般不需要多线程加速。flag大于4096节而小于5000节时多线程乘法耗费的总时间仅仅比单线程乘法运算快100多毫秒。为了实现跨操作系统平台的多线程,我们需要为类UNIX与NT分别编写多线程函数的实现方式。NT系统下,我们可以参考MSDN。类UNIX系统下,可以使用POSIX pthread接口。关于多线程并行计算,以及POSIX pthread的资料,很多地方均有参考,这里我推荐美国劳伦斯利弗莫尔国家实验室的公开资料:https://computing.llnl.gov/tutorials/pthreads/
        一个最基本的大数计算库,实现至此已经尚可,但是,为了GILIB的实用价值,我们需要为其编写数学库。
在 GILIB_MATH 文件中包含了若干函数,其中有素数检测(Miller-Rabin方法)、随机数发生器、求平方根、最大公因数 a的b次幂与c的模数函数(用在密码学方面)、阶乘、幂。关于素数检测,Miller-Rabin强伪素数并不是真正的素数,尤其是整数位数一大,检测结果有少许疏漏。建议自己完成Lehmer检测,进而实现Baillie-PSW检测以增大检测结果正确概率。求幂函数可以改进为二进方法。最大公因数我使用了最简单的Euclid方法,大家可以改进为二进方法,或者Lehmer加速方法。数学库的实现、优化可不止一日两日能完成,需要时间的积累。不管怎样,维基百科上就可以找到许多通俗易懂的指导。这里推荐一个网站:http://mathworld.wolfram.com
        在 GILIB_EXTMEM.C 文件中,实现了可以将内存中的GINT整数保存至外部存储器上的功能。gemSaveGI 函数将GINT保存为 xxx.gid 文件。GINT data 文件格式在 SN_FGINT 结构体中有声明。注意:32位平台与64位平台的GINT data文件不可以通用。一个GINT data文件只能保存一个GINT整数。这里如果大家使用文本编辑器打开一个GINT data文件,并且文本编辑器编码格式选择为UTF-8,就会看到我留的“彩蛋”之一。
        GILIB库内其他一些文件,因为时间关系,这些文件并没有经过很好的测试,bug很多。这些文件有:GILIB_INTERFACE.CXX 是 GILIB 库的C++接口。GILIB_UTILITIES.C 是高精度浮点数运算工具。
        GILIB库提供了三个demo。第一个是一个控制台程序,理论上可以在任何平台上编译。第二个是一个Win32计算器,只能编译运行在NT平台上。第三个是一个动态链接库编译的例子。动态链接库(DLL)文件只能编译运行于NT平台。
以下是GILIB的下载地址。(注意:GILIB遵循LGPLv3协议,这意味着他的使用权限比GPL,“L” [良い!])
https://sourceforge.net/projects/gilib/files/
In the beginning I was not the best.
And the world was also not the best.
But I still know that I am who I am.
Because I think that it is good.
I have been working hard.
I have been keeping growth with the world.
And it was so.

0

主题

22

帖子

57

积分

用户组: 小·技术宅

UID
1919
精华
0
威望
0 点
宅币
35 个
贡献
0 次
宅之契约
0 份
在线时间
2 小时
注册时间
2016-8-24
发表于 2016-8-24 21:06:58 | 显示全部楼层
新人学习一下

0

主题

22

帖子

57

积分

用户组: 小·技术宅

UID
1919
精华
0
威望
0 点
宅币
35 个
贡献
0 次
宅之契约
0 份
在线时间
2 小时
注册时间
2016-8-24
发表于 2016-8-24 21:07:24 | 显示全部楼层
新人学习一下

本版积分规则

QQ|申请友链|Archiver|手机版|小黑屋|技术宅的结界 ( 滇ICP备16008837号|网站地图

GMT+8, 2018-4-26 00:13 , Processed in 0.096814 second(s), 21 queries , Gzip On, Memcache On.

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表