数学研发网设为首页收藏本站

数学研发论坛

 找回密码
 欢迎注册
查看: 985|回复: 22

[原创] 基于x64的二进制乘法库-超越最新的GMP-6.1.2

[复制链接]
发表于 2017-4-26 16:38:55 | 显示全部楼层 |阅读模式

马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。

您需要 登录 才可以下载或查看,没有帐号?欢迎注册

x
本帖最后由 Ickiverar 于 2017-4-26 16:44 编辑

我此时只想仰天长笑。
经过我不知道多少日夜的努力,现在我的64位乘法库已经超越了GMP的最新版6.1.2在windows-x64(我的机子)下编译的结果。
结果:

这幅图给出了直到700万字规模的几个算法的计时。
蓝线是我原来的成果, http://bbs.emath.ac.cn/thread-8680-1-1.html 这个帖子里有详情。
紫线是Mathematica 8.0.4,也就是MMA,专业的数学软件大家都知道。
黄线是 GMP6.1.2,是我从GMP官网上下的最新版,在我的机子上编译的64位版本。
绿线是我现在的成果。
一字就是64个二进制位,相当于19.3十进制位。700万字意味着这是一个 1.35亿位乘1.35亿位得到2.7亿位结果的乘法。
从大约70万字开始,我的算法快于GMP,平均快 9% 左右。
当然在小字长上我不如GMP,最多能慢一半(用时多一半)。因为GMP实现了一长串十来个Toom-Cook算法,而我为了偷懒只写了 Toom-Cook22, 32, 33和42,这四个。
说实话我最底层的汇编代码借鉴了GMP的源码……但,GMP的FFT(SSA算法)实现得非常糟糕,它的源码我基本看不下去,我重写了SSA算法。事实证明在FFT的字长范围我实现得还不错。
给个度盘的链接吧,包括我的库,GMP的库以及一个测试例程。不知道为什么我发不了附件
https://pan.baidu.com/s/1bpnb3pP

点评

发帖(或编辑)时,去掉右上角的“纯文本”勾选项,切换到所见即所得模式试试。。。  发表于 2017-4-26 16:50
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-4-26 17:12:42 | 显示全部楼层
恭喜楼主,贺喜楼主。
实不相瞒,我也在实现一个同时支持\(10^9\)进制和\(2^{30}\)进制的大数库,目前只支持32位,将来计划支持64位。我的库主要目标是追求极致的性能。为了达到高性能,不惜花费大量的时间,并使用了超多的汇编代码(超过10万行)。目前在在小字长(1000*30比特以下)的情况下,完胜GMP6.12.  即使在10000*30比特下,也只比GMP慢一点点。在大数乘法方面,目前已经实现了9种Toom-Cook算法,暂时不打算实现NTT和SSA算法,而将主要精力放在任意精度对数函数的实现上。

点评

我是非常佩服能狠下心来实现Toom-Cook一系列算法的人的。像我,基本上就把Toom-Cook后面那几对数字写在纸上(22,32,33,42,43,52,44,53,62,54,63,72,6h,8h...),画一画最优算法图,就不想实现了……又多又难写…(哭  发表于 2017-4-26 17:32
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-4-26 17:38:50 | 显示全部楼层
liangbch 发表于 2017-4-26 17:12
恭喜楼主,贺喜楼主。
实不相瞒,我也在实现一个同时支持\(10^9\)进制和\(2^{30}\)进制的大数库,目前只支 ...

对我而言,写这些程序同样是比较困难的任务。我实现了Toom22,toom32,toom33,toom42,toom43,toom44,toom52,toom53,toom62 共9种Toom-Cook算法,完全重新实现,中间结果的内存布局和GMP不同,比GMP占用的内存空间更省。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-4-27 15:34:45 | 显示全部楼层
在我的电脑上运行你的程序,测试结果表明,你的实现并无绝对优势,和GMP相比,互有胜负。
  1. D:\download\测试程序\测试程序>test.exe
  2. seed=1493277974
  3. k>1 means ilmp is faster, k<1 means gmp is faster.
  4. n=         1    gmp=          12ns      ilmp=           8ns     k=  1.621015
  5. n=         6    gmp=          38ns      ilmp=          36ns     k=  1.047551
  6. n=        11    gmp=         101ns      ilmp=         107ns     k=  0.943279
  7. n=        17    gmp=         224ns      ilmp=         246ns     k=  0.909822
  8. n=        23    gmp=         384ns      ilmp=         410ns     k=  0.938326
  9. n=        30    gmp=         624ns      ilmp=         636ns     k=  0.980198
  10. n=        38    gmp=         894ns      ilmp=         974ns     k=  0.918196
  11. n=        46    gmp=        1267ns      ilmp=        1462ns     k=  0.867156
  12. n=        55    gmp=        1892ns      ilmp=        1843ns     k=  1.026708
  13. n=        65    gmp=        2246ns      ilmp=        2617ns     k=  0.858353
  14. n=        76    gmp=        3014ns      ilmp=        3361ns     k=  0.896703
  15. n=        88    gmp=        3902ns      ilmp=        4162ns     k=  0.937578
  16. n=       101    gmp=        4560ns      ilmp=        5174ns     k=  0.881345
  17. n=       116    gmp=        5922ns      ilmp=        6322ns     k=  0.936774
  18. n=       132    gmp=        6920ns      ilmp=        7561ns     k=  0.915310
  19. n=       150    gmp=        8459ns      ilmp=        9436ns     k=  0.896426
  20. n=       170    gmp=       10084ns      ilmp=       11739ns     k=  0.859010
  21. n=       192    gmp=       12927ns      ilmp=       14504ns     k=  0.891302
  22. n=       216    gmp=       14802ns      ilmp=       17875ns     k=  0.828108
  23. n=       242    gmp=       20788ns      ilmp=       20724ns     k=  1.003088
  24. n=       271    gmp=       21229ns      ilmp=       24055ns     k=  0.882547
  25. n=       303    gmp=       24990ns      ilmp=       28498ns     k=  0.876910
  26. n=       338    gmp=       28940ns      ilmp=       33688ns     k=  0.859062
  27. n=       376    gmp=       33555ns      ilmp=       41815ns     k=  0.802470
  28. n=       418    gmp=       40058ns      ilmp=       45008ns     k=  0.890023
  29. n=       464    gmp=       47097ns      ilmp=       55326ns     k=  0.851271
  30. n=       515    gmp=       53844ns      ilmp=       61917ns     k=  0.869620
  31. n=       571    gmp=       62324ns      ilmp=       85177ns     k=  0.731697
  32. n=       633    gmp=       77967ns      ilmp=       90826ns     k=  0.858419
  33. n=       701    gmp=       84184ns      ilmp=      111939ns     k=  0.752049
  34. n=       776    gmp=       93546ns      ilmp=      119015ns     k=  0.785999
  35. n=       858    gmp=      115892ns      ilmp=      138998ns     k=  0.833766
  36. n=       948    gmp=      122903ns      ilmp=      162649ns     k=  0.755631
  37. n=      1047    gmp=      143072ns      ilmp=      186824ns     k=  0.765812
  38. n=      1156    gmp=      165828ns      ilmp=      214333ns     k=  0.773692
  39. n=      1276    gmp=      196053ns      ilmp=      249764ns     k=  0.784954
  40. n=      1408    gmp=      210600ns      ilmp=      284253ns     k=  0.740891
  41. n=      1553    gmp=      264641ns      ilmp=      327403ns     k=  0.808305
  42. n=      1713    gmp=      286842ns      ilmp=      387051ns     k=  0.741097
  43. n=      1889    gmp=      329077ns      ilmp=      446205ns     k=  0.737502
  44. n=      2082    gmp=      467912ns      ilmp=      546575ns     k=  0.856080
  45. n=      2295    gmp=      544486ns      ilmp=      573932ns     k=  0.948695
  46. n=      2529    gmp=      575779ns      ilmp=      589931ns     k=  0.976011
  47. n=      2786    gmp=      688270ns      ilmp=      679206ns     k=  1.013345
  48. n=      3069    gmp=      720444ns      ilmp=      806878ns     k=  0.892878
  49. n=      3380    gmp=      811015ns      ilmp=      892991ns     k=  0.908200
  50. n=      3723    gmp=      884408ns      ilmp=      990613ns     k=  0.892789
  51. n=      4100    gmp=     1018786ns      ilmp=     1231631ns     k=  0.827184
  52. n=      4515    gmp=     1293999ns      ilmp=     1201320ns     k=  1.077147
  53. n=      4971    gmp=     1276848ns      ilmp=     1366523ns     k=  0.934378
  54. n=      5473    gmp=     1387762ns      ilmp=     1443015ns     k=  0.961710
  55. n=      6025    gmp=     1565099ns      ilmp=     1525972ns     k=  1.025641
  56. n=      6632    gmp=     1694864ns      ilmp=     1740041ns     k=  0.974037
  57. n=      7300    gmp=     1862608ns      ilmp=     2049918ns     k=  0.908626
  58. n=      8035    gmp=     2145797ns      ilmp=     2064504ns     k=  1.039376
  59. n=      8843    gmp=     2521734ns      ilmp=     2608956ns     k=  0.966568
  60. n=      9732    gmp=     2729006ns      ilmp=     2653740ns     k=  1.028362
  61. n=     10710    gmp=     2905987ns      ilmp=     3319500ns     k=  0.875429
  62. n=     11786    gmp=     3496002ns      ilmp=     3256494ns     k=  1.073548
  63. n=     12969    gmp=     3592527ns      ilmp=     3781459ns     k=  0.950037
  64. n=     14270    gmp=     4112878ns      ilmp=     3939986ns     k=  1.043881
  65. n=     15702    gmp=     4228194ns      ilmp=     4480131ns     k=  0.943766
  66. n=     17277    gmp=     4793725ns      ilmp=     5763703ns     k=  0.831709
  67. n=     19009    gmp=     5622742ns      ilmp=     6321174ns     k=  0.889509
  68. n=     20914    gmp=     6267117ns      ilmp=     7257836ns     k=  0.863497
  69. n=     23010    gmp=     7047754ns      ilmp=     6851196ns     k=  1.028690
  70. n=     25316    gmp=     7686810ns      ilmp=     8646967ns     k=  0.888960
  71. n=     27852    gmp=     8822022ns      ilmp=     9560754ns     k=  0.922733
  72. n=     30642    gmp=     9730001ns      ilmp=     9866100ns     k=  0.986205
  73. n=     33711    gmp=    11806869ns      ilmp=    13034998ns     k=  0.905782
  74. n=     37087    gmp=    14785851ns      ilmp=    14274669ns     k=  1.035810
  75. n=     40800    gmp=    14711108ns      ilmp=    14564730ns     k=  1.010050
  76. n=     44885    gmp=    14675430ns      ilmp=    16245389ns     k=  0.903360
  77. n=     49378    gmp=    17524302ns      ilmp=    19957500ns     k=  0.878081
  78. n=     54320    gmp=    18171355ns      ilmp=    20565260ns     k=  0.883595
  79. n=     59757    gmp=    21623690ns      ilmp=    21640439ns     k=  0.999226
  80. n=     65737    gmp=    26945191ns      ilmp=    27814515ns     k=  0.968746
  81. n=     72315    gmp=    28145044ns      ilmp=    28892253ns     k=  0.974138
  82. n=     79551    gmp=    29616903ns      ilmp=    32303177ns     k=  0.916842
  83. n=     87511    gmp=    32707198ns      ilmp=    36009928ns     k=  0.908283
  84. n=     96267    gmp=    38551550ns      ilmp=    35829139ns     k=  1.075983
  85. n=    105898    gmp=    47034065ns      ilmp=    46249187ns     k=  1.016971
  86. n=    116492    gmp=    49244010ns      ilmp=    48172893ns     k=  1.022235
  87. n=    128146    gmp=    52930197ns      ilmp=    46322172ns     k=  1.142654
  88. n=    140965    gmp=    61661354ns      ilmp=    64796445ns     k=  0.951616
  89. n=    155066    gmp=    69173000ns      ilmp=    68268951ns     k=  1.013242
  90. n=    170577    gmp=    74071286ns      ilmp=    83825307ns     k=  0.883639
  91. n=    187639    gmp=    82018278ns      ilmp=    90033694ns     k=  0.910973
  92. n=    206407    gmp=   100268819ns      ilmp=   112791650ns     k=  0.888974
  93. n=    227052    gmp=   112261152ns      ilmp=   107102429ns     k=  1.048166
  94. n=    249762    gmp=   122833906ns      ilmp=   105580855ns     k=  1.163411
  95. n=    274743    gmp=   139121327ns      ilmp=   150434341ns     k=  0.924798
  96. n=    302222    gmp=   159464139ns      ilmp=   163827009ns     k=  0.973369
  97. n=    332449    gmp=   178752797ns      ilmp=   192130925ns     k=  0.930370
  98. n=    365698    gmp=   205912753ns      ilmp=   187744962ns     k=  1.096768
  99. n=    402272    gmp=   226229478ns      ilmp=   217542227ns     k=  1.039934
  100. n=    442504    gmp=   262347815ns      ilmp=   217557622ns     k=  1.205877
  101. n=    486759    gmp=   253564770ns      ilmp=   216938815ns     k=  1.168831
  102. n=    535439    gmp=   290248031ns      ilmp=   339676073ns     k=  0.854485
  103. n=    588987    gmp=   330180993ns      ilmp=   340253398ns     k=  0.970397
  104. n=    647890    gmp=   351209325ns      ilmp=   392852881ns     k=  0.893997
  105. n=    712684    gmp=   377366446ns      ilmp=   425663785ns     k=  0.886536
  106. n=    783957    gmp=   444466634ns      ilmp=   458319452ns     k=  0.969775
  107. n=    862357    gmp=   509146759ns      ilmp=   505154661ns     k=  1.007903
  108. n=    948597    gmp=   532963788ns      ilmp=   540649914ns     k=  0.985784
  109. n=   1043461    gmp=   580735974ns      ilmp=   560308916ns     k=  1.036457
  110. n=   1147812    gmp=   657908140ns      ilmp=   710752237ns     k=  0.925650
  111. n=   1262598    gmp=   730722995ns      ilmp=   739181882ns     k=  0.988556
  112. n=   1388862    gmp=   835570006ns      ilmp=   847162274ns     k=  0.986316
  113. n=   1527753    gmp=   936996256ns      ilmp=   900178714ns     k=  1.040900
  114. n=   1680533    gmp=  1020859410ns      ilmp=  1025445513ns     k=  0.995528
  115. n=   1848591    gmp=  1172358169ns      ilmp=  1143482914ns     k=  1.025252
  116. n=   2033455    gmp=  1263941234ns      ilmp=  1109889414ns     k=  1.138799
  117. n=   2236805    gmp=  1434605915ns      ilmp=  1415366864ns     k=  1.013593
  118. n=   2460490    gmp=  1470467235ns      ilmp=  1512510215ns     k=  0.972203
  119. n=   2706544    gmp=  1694172305ns      ilmp=  1686090604ns     k=  1.004793
  120. n=   2977203    gmp=  1829641942ns      ilmp=  1894271605ns     k=  0.965882
  121. n=   3274928    gmp=  2109515221ns      ilmp=  2132744233ns     k=  0.989108
  122. n=   3602425    gmp=  2422466257ns      ilmp=  2280946249ns     k=  1.062044
  123. n=   3962672    gmp=  2801679717ns      ilmp=  2493125764ns     k=  1.123762
  124. n=   4358944    gmp=  2848544434ns      ilmp=  2937701021ns     k=  0.969651
  125. n=   4794843    gmp=  3131557508ns      ilmp=  3264304307ns     k=  0.959334
  126. n=   5274332    gmp=  3501409743ns      ilmp=  3621818027ns     k=  0.966755
  127. n=   5801770    gmp=  3962427230ns      ilmp=  3998511783ns     k=  0.990976
  128. n=   6381952    gmp=  4412224075ns      ilmp=  4396991236ns     k=  1.003464
复制代码

点评

那个测试用的exe,是gmp64和我的实现比较的。我没有试过gmp32,原则上同样规模的乘法,32应该慢一倍左右。我不太明白你的话,难道是说gmp64和32完全一样?那为何还要开发64位版本?(我没有编译gmp32)  发表于 2017-5-3 14:05
我弄错了,64位GMP的乘法的速度正好等于乘数为2倍长度的32位GMP的乘法的速度。  发表于 2017-4-28 09:09
从测试结果来看,貌似64位GMP和32位GMP一样快。楼主可否测试一下你的GMP64和GMP 32的的性能。  发表于 2017-4-27 15:56
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-4-28 17:10:06 来自手机 | 显示全部楼层
这个结果很正常,不同的机器配置不同,所以在不同机器上性能结果会有区别是正常的。这充分说明了测试的重要性
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-4-29 10:35:51 | 显示全部楼层
对。我也时常注意到这种情况。两个不同的实现,在测试软硬件环境A,实现1快于实现2,在另一个环境则相反。甚至在同样的硬件环境,只是测试上下文环境稍有不同,就会得到不同的结果。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-7-8 09:41:04 | 显示全部楼层
这几天把除法写完了,主要是3/2试商法,分治法和牛顿迭代法。
最终速度图像如下,我的CPU是 intel i5 3320M@2.6GHz 的,L1=128k。内存 4G DDR3 @1.6GHz.
在大于百万的字长,我的库稳定快于GMP,但不多。按楼上的测试,基本上换个平台这个优势就不稳定了。


在几十万字长的范围内,因为FFT类算法的通病,只支持特定字长乘法,所以一些字长范围都要向上ceiling到最小可行字长。
所以可以看到一系列的平台。
同时我的除法和GMP的除法在很大一段范围内互相交错……这是因为平台范围刚好错开。


在几万的范围,我的除法稳定慢于GMP。乘法基本没有差别。


几千。
同样是除法明显慢,乘法没有区别。


几百,Toom-Cook范围。
完全慢于GMP。
我认为这是可以接受的,考虑到GMP实现了13种Toom-Cook算法来加速这一段,而我非常懒只写了4个(而且未来也不会再写更多,我就是这么懒)。


几十,汇编范围。
我的汇编技术基本上是向GMP学习的,比如那个3/2试商法,还有各种汇编过程中的四路循环展开以及多通道跳入等等,所以速度差别不大。
它的代码充满了各种诡异迭代而难以理解的宏,而且还是用AT&T语法写的,简直是一场灾难。我花了很多个晚上才理顺了那些玩意,然后用intel语法重写。
(还发现了它乘法 mul_basecase 以及平方 sqr_basecase 汇编里的一个小bug,读内存越界。因为没有越界写,所以如果没有当场访问冲突崩溃,就不会造成什么后果。)


总之,除了乘法的SSA/FFT算法是完全由我自己写的,其他代码都或多或少借鉴了GMP。
一个原则是,一切算法必须经过我亲手证明其正确性,才能把它写进ilmp。(所以借鉴绝不是复制粘贴代码
(有些算法真的很难证明,我觉得最难的就是那个3/2试商。不知道为什么没人回复,太长懒得看吗
有些实在是没法证明,那只能我自己想一个等效的算法补上去。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
发表于 2017-7-8 15:06:05 | 显示全部楼层
楼主辛苦了,性能真的很不错。
毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-7-12 13:31:38 | 显示全部楼层
最后就剩下 Base Conversion。
内部使用$2^64$进制,使用十进制输入输出就必须进行转进制。
输出基本就是把一个高精度数不断地除以10的某个次方,余数和商就是输出的低位部分和高位部分。
输入就是把高位部分乘以10的某个次方的积加上低位部分。
由于除法总是慢于乘法,所以输出总是慢于输入。
我对输出做了一点点优化:除法就是乘以逆(倒数),所以我提前算了所有10的高次幂的逆,这样就把除法转换成了乘法,此时理论输出耗时是输入的两倍(需要两次乘法,一次把被除数乘以逆得到商,一次把商乘以除数(然后和被除数做差)来得到余数)。
所以可以看到在我的输入和GMP差不多快的情况下,大字长的输出比GMP快了10%~15%,而且基本就是两倍的输入耗时。

我还做了n=100万以上的测试,对这样大的n,我的输出快于GMP20%~30%。
因为过于耗时,点没取得很密,画不出来图,列个表:

n        ilmp_in   gmp_in  ilmp_out  gmp_out   输出快于GMP
1000000     1995     1850     3864     4605         19.17%
2000000     4513     4338     8767    10697         22.01%
4000000    10219     9969    19908    25096         26.06%
8000000    22964    22558    44955    58039         29.10%


假设除法的耗时是乘法的2.6倍(实测GMP在大字长范围大致如此),那么理论上我就会快30%。

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
 楼主| 发表于 2017-7-24 11:31:31 | 显示全部楼层
这应该是最后四张速度测试图。
ilmp 的核心到这里就结束了。
也许以后会写 gcd,divexact 什么的新函数,但现在我可以休息了。

乘,除,平方,开根:

新写了开根的函数。测速过程是求一个 2n-limb 整数 x 的平方根 s 以及余数 r,满足 x = s*s + r 。
开根使用的是 Karatsuba 算法。

对数比例图:

对 n=1 到 700万,这是各种运算耗时与一次乘法的耗时比。
可以看到在 n 很大时,一次平方大概耗时相当于 0.67 次乘法,除法是 2.8,平方根是 3.

进制转换输入输出:

该说的楼上都说过了,花了很长时间测到了700万。

求2的平方根(不输出):

求2的近似平方根时,舍入可能是 floor 也可能是 round,保证最后一位与精确值不会差 1 以上即可。当然此时也不用求余数了。
gmp 使用的算法仍然是 Karatsuba,将 2 后面补上 2n 个 0,然后看作一个大整数去求它的根。
我使用的算法是 牛顿迭代法求逆平方根(1/√x),再乘原数得到平方根 (简称牛顿法)。
显然,如果原数的精度很低(比如 2 只有 1 位精度),而我们要求很高精度的根,那牛顿法是远远好于Karatsuba方法的,实测加速达到了 100% 以上(如图)。
进一步测速表明,当根的精度大约是原数精度的十倍以上时,牛顿法要更好(否则 Karatsuba 更好)。
在 ilmp 里,Karatsuba 和 牛顿法 是由 sqrt 函数自动选择的,调用者并不需要担心到底使用何种算法。

毋因群疑而阻独见  毋任己意而废人言
毋私小惠而伤大体  毋借公论以快私情
您需要登录后才可以回帖 登录 | 欢迎注册

本版积分规则

小黑屋|手机版|Archiver|数学研发网 ( 苏ICP备07505100号  

GMT+8, 2017-9-20 08:27 , Processed in 0.340190 second(s), 22 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.

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