【C】大数计算GILIB库

cyycoish

        大整数运算在现代密码学与高精度计算领域有着非常重要的作用。
        在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/

jason

新人学习一下

jason

新人学习一下