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

Ickiverar

我此时只想仰天长笑。
经过我不知道多少日夜的努力，现在我的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

补充内容 (2020-3-29 14:11):
emmm百度盘和下面的测速图全挂了，补个链接，包括ilmp的库，圆周率计算器（见https://bbs.emath.ac.cn/thread-8760-1-1.html）的例程，和测速图。
https://pan.baidu.com/s/1e2CB23nK0s2i3Fn7MBMimA
提取码py8f

liangbch

恭喜楼主，贺喜楼主。
实不相瞒，我也在实现一个同时支持\(10^9\)进制和\(2^{30}\)进制的大数库，目前只支持32位，将来计划支持64位。我的库主要目标是追求极致的性能。为了达到高性能，不惜花费大量的时间，并使用了超多的汇编代码（超过10万行）。目前在在小字长（1000*30比特以下）的情况下，完胜GMP6.12.  即使在10000*30比特下，也只比GMP慢一点点。在大数乘法方面，目前已经实现了9种Toom-Cook算法，暂时不打算实现NTT和SSA算法，而将主要精力放在任意精度对数函数的实现上。

liangbch

liangbch 发表于 2017-4-26 17:12
恭喜楼主，贺喜楼主。
实不相瞒，我也在实现一个同时支持\(10^9\)进制和\(2^{30}\)进制的大数库，目前只支 ...

对我而言，写这些程序同样是比较困难的任务。我实现了Toom22,toom32,toom33,toom42,toom43,toom44,toom52,toom53,toom62 共9种Toom-Cook算法，完全重新实现，中间结果的内存布局和GMP不同，比GMP占用的内存空间更省。

liangbch

在我的电脑上运行你的程序，测试结果表明，你的实现并无绝对优势，和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
     

复制代码

mathe

这个结果很正常，不同的机器配置不同，所以在不同机器上性能结果会有区别是正常的。这充分说明了测试的重要性

liangbch

对。我也时常注意到这种情况。两个不同的实现，在测试软硬件环境A，实现1快于实现2，在另一个环境则相反。甚至在同样的硬件环境，只是测试上下文环境稍有不同，就会得到不同的结果。

Ickiverar

这几天把除法写完了，主要是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试商。不知道为什么没人回复，太长懒得看吗）
有些实在是没法证明，那只能我自己想一个等效的算法补上去。

liangbch

楼主辛苦了，性能真的很不错。

Ickiverar

最后就剩下 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%。

Ickiverar

这应该是最后四张速度测试图。
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 函数自动选择的，调用者并不需要担心到底使用何种算法。