基于x64的二进制乘法库-超越最新的GMP-6.1.2
本帖最后由 Ickiverar 于 2017-4-26 16:44 编辑我此时只想仰天长笑。
经过我不知道多少日夜的努力，现在我的64位乘法库已经超越了GMP的最新版6.1.2在windows-x64（我的机子）下编译的结果。
结果：
http://chuantu.biz/t5/74/1493194311x2890149494.png
这幅图给出了直到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 恭喜楼主，贺喜楼主。
实不相瞒，我也在实现一个同时支持\(10^9\)进制和\(2^{30}\)进制的大数库，目前只支持32位，将来计划支持64位。我的库主要目标是追求极致的性能。为了达到高性能，不惜花费大量的时间，并使用了超多的汇编代码（超过10万行）。目前在在小字长（1000*30比特以下）的情况下，完胜GMP6.12.即使在10000*30比特下，也只比GMP慢一点点。在大数乘法方面，目前已经实现了9种Toom-Cook算法，暂时不打算实现NTT和SSA算法，而将主要精力放在任意精度对数函数的实现上。 liangbch 发表于 2017-4-26 17:12
恭喜楼主，贺喜楼主。
实不相瞒，我也在实现一个同时支持\(10^9\)进制和\(2^{30}\)进制的大数库，目前只支 ...
对我而言，写这些程序同样是比较困难的任务。我实现了Toom22,toom32,toom33,toom42,toom43,toom44,toom52,toom53,toom62 共9种Toom-Cook算法，完全重新实现，中间结果的内存布局和GMP不同，比GMP占用的内存空间更省。 在我的电脑上运行你的程序，测试结果表明，你的实现并无绝对优势，和GMP相比，互有胜负。
D:\download\测试程序\测试程序>test.exe
seed=1493277974
k>1 means ilmp is faster, k<1 means gmp is faster.
n= 1 gmp= 12ns ilmp= 8ns k=1.621015
n= 6 gmp= 38ns ilmp= 36ns k=1.047551
n= 11 gmp= 101ns ilmp= 107ns k=0.943279
n= 17 gmp= 224ns ilmp= 246ns k=0.909822
n= 23 gmp= 384ns ilmp= 410ns k=0.938326
n= 30 gmp= 624ns ilmp= 636ns k=0.980198
n= 38 gmp= 894ns ilmp= 974ns k=0.918196
n= 46 gmp= 1267ns ilmp= 1462ns k=0.867156
n= 55 gmp= 1892ns ilmp= 1843ns k=1.026708
n= 65 gmp= 2246ns ilmp= 2617ns k=0.858353
n= 76 gmp= 3014ns ilmp= 3361ns k=0.896703
n= 88 gmp= 3902ns ilmp= 4162ns k=0.937578
n= 101 gmp= 4560ns ilmp= 5174ns k=0.881345
n= 116 gmp= 5922ns ilmp= 6322ns k=0.936774
n= 132 gmp= 6920ns ilmp= 7561ns k=0.915310
n= 150 gmp= 8459ns ilmp= 9436ns k=0.896426
n= 170 gmp= 10084ns ilmp= 11739ns k=0.859010
n= 192 gmp= 12927ns ilmp= 14504ns k=0.891302
n= 216 gmp= 14802ns ilmp= 17875ns k=0.828108
n= 242 gmp= 20788ns ilmp= 20724ns k=1.003088
n= 271 gmp= 21229ns ilmp= 24055ns k=0.882547
n= 303 gmp= 24990ns ilmp= 28498ns k=0.876910
n= 338 gmp= 28940ns ilmp= 33688ns k=0.859062
n= 376 gmp= 33555ns ilmp= 41815ns k=0.802470
n= 418 gmp= 40058ns ilmp= 45008ns k=0.890023
n= 464 gmp= 47097ns ilmp= 55326ns k=0.851271
n= 515 gmp= 53844ns ilmp= 61917ns k=0.869620
n= 571 gmp= 62324ns ilmp= 85177ns k=0.731697
n= 633 gmp= 77967ns ilmp= 90826ns k=0.858419
n= 701 gmp= 84184ns ilmp= 111939ns k=0.752049
n= 776 gmp= 93546ns ilmp= 119015ns k=0.785999
n= 858 gmp= 115892ns ilmp= 138998ns k=0.833766
n= 948 gmp= 122903ns ilmp= 162649ns k=0.755631
n= 1047 gmp= 143072ns ilmp= 186824ns k=0.765812
n= 1156 gmp= 165828ns ilmp= 214333ns k=0.773692
n= 1276 gmp= 196053ns ilmp= 249764ns k=0.784954
n= 1408 gmp= 210600ns ilmp= 284253ns k=0.740891
n= 1553 gmp= 264641ns ilmp= 327403ns k=0.808305
n= 1713 gmp= 286842ns ilmp= 387051ns k=0.741097
n= 1889 gmp= 329077ns ilmp= 446205ns k=0.737502
n= 2082 gmp= 467912ns ilmp= 546575ns k=0.856080
n= 2295 gmp= 544486ns ilmp= 573932ns k=0.948695
n= 2529 gmp= 575779ns ilmp= 589931ns k=0.976011
n= 2786 gmp= 688270ns ilmp= 679206ns k=1.013345
n= 3069 gmp= 720444ns ilmp= 806878ns k=0.892878
n= 3380 gmp= 811015ns ilmp= 892991ns k=0.908200
n= 3723 gmp= 884408ns ilmp= 990613ns k=0.892789
n= 4100 gmp= 1018786ns ilmp= 1231631ns k=0.827184
n= 4515 gmp= 1293999ns ilmp= 1201320ns k=1.077147
n= 4971 gmp= 1276848ns ilmp= 1366523ns k=0.934378
n= 5473 gmp= 1387762ns ilmp= 1443015ns k=0.961710
n= 6025 gmp= 1565099ns ilmp= 1525972ns k=1.025641
n= 6632 gmp= 1694864ns ilmp= 1740041ns k=0.974037
n= 7300 gmp= 1862608ns ilmp= 2049918ns k=0.908626
n= 8035 gmp= 2145797ns ilmp= 2064504ns k=1.039376
n= 8843 gmp= 2521734ns ilmp= 2608956ns k=0.966568
n= 9732 gmp= 2729006ns ilmp= 2653740ns k=1.028362
n= 10710 gmp= 2905987ns ilmp= 3319500ns k=0.875429
n= 11786 gmp= 3496002ns ilmp= 3256494ns k=1.073548
n= 12969 gmp= 3592527ns ilmp= 3781459ns k=0.950037
n= 14270 gmp= 4112878ns ilmp= 3939986ns k=1.043881
n= 15702 gmp= 4228194ns ilmp= 4480131ns k=0.943766
n= 17277 gmp= 4793725ns ilmp= 5763703ns k=0.831709
n= 19009 gmp= 5622742ns ilmp= 6321174ns k=0.889509
n= 20914 gmp= 6267117ns ilmp= 7257836ns k=0.863497
n= 23010 gmp= 7047754ns ilmp= 6851196ns k=1.028690
n= 25316 gmp= 7686810ns ilmp= 8646967ns k=0.888960
n= 27852 gmp= 8822022ns ilmp= 9560754ns k=0.922733
n= 30642 gmp= 9730001ns ilmp= 9866100ns k=0.986205
n= 33711 gmp= 11806869ns ilmp= 13034998ns k=0.905782
n= 37087 gmp= 14785851ns ilmp= 14274669ns k=1.035810
n= 40800 gmp= 14711108ns ilmp= 14564730ns k=1.010050
n= 44885 gmp= 14675430ns ilmp= 16245389ns k=0.903360
n= 49378 gmp= 17524302ns ilmp= 19957500ns k=0.878081
n= 54320 gmp= 18171355ns ilmp= 20565260ns k=0.883595
n= 59757 gmp= 21623690ns ilmp= 21640439ns k=0.999226
n= 65737 gmp= 26945191ns ilmp= 27814515ns k=0.968746
n= 72315 gmp= 28145044ns ilmp= 28892253ns k=0.974138
n= 79551 gmp= 29616903ns ilmp= 32303177ns k=0.916842
n= 87511 gmp= 32707198ns ilmp= 36009928ns k=0.908283
n= 96267 gmp= 38551550ns ilmp= 35829139ns k=1.075983
n= 105898 gmp= 47034065ns ilmp= 46249187ns k=1.016971
n= 116492 gmp= 49244010ns ilmp= 48172893ns k=1.022235
n= 128146 gmp= 52930197ns ilmp= 46322172ns k=1.142654
n= 140965 gmp= 61661354ns ilmp= 64796445ns k=0.951616
n= 155066 gmp= 69173000ns ilmp= 68268951ns k=1.013242
n= 170577 gmp= 74071286ns ilmp= 83825307ns k=0.883639
n= 187639 gmp= 82018278ns ilmp= 90033694ns k=0.910973
n= 206407 gmp= 100268819ns ilmp= 112791650ns k=0.888974
n= 227052 gmp= 112261152ns ilmp= 107102429ns k=1.048166
n= 249762 gmp= 122833906ns ilmp= 105580855ns k=1.163411
n= 274743 gmp= 139121327ns ilmp= 150434341ns k=0.924798
n= 302222 gmp= 159464139ns ilmp= 163827009ns k=0.973369
n= 332449 gmp= 178752797ns ilmp= 192130925ns k=0.930370
n= 365698 gmp= 205912753ns ilmp= 187744962ns k=1.096768
n= 402272 gmp= 226229478ns ilmp= 217542227ns k=1.039934
n= 442504 gmp= 262347815ns ilmp= 217557622ns k=1.205877
n= 486759 gmp= 253564770ns ilmp= 216938815ns k=1.168831
n= 535439 gmp= 290248031ns ilmp= 339676073ns k=0.854485
n= 588987 gmp= 330180993ns ilmp= 340253398ns k=0.970397
n= 647890 gmp= 351209325ns ilmp= 392852881ns k=0.893997
n= 712684 gmp= 377366446ns ilmp= 425663785ns k=0.886536
n= 783957 gmp= 444466634ns ilmp= 458319452ns k=0.969775
n= 862357 gmp= 509146759ns ilmp= 505154661ns k=1.007903
n= 948597 gmp= 532963788ns ilmp= 540649914ns k=0.985784
n= 1043461 gmp= 580735974ns ilmp= 560308916ns k=1.036457
n= 1147812 gmp= 657908140ns ilmp= 710752237ns k=0.925650
n= 1262598 gmp= 730722995ns ilmp= 739181882ns k=0.988556
n= 1388862 gmp= 835570006ns ilmp= 847162274ns k=0.986316
n= 1527753 gmp= 936996256ns ilmp= 900178714ns k=1.040900
n= 1680533 gmp=1020859410ns ilmp=1025445513ns k=0.995528
n= 1848591 gmp=1172358169ns ilmp=1143482914ns k=1.025252
n= 2033455 gmp=1263941234ns ilmp=1109889414ns k=1.138799
n= 2236805 gmp=1434605915ns ilmp=1415366864ns k=1.013593
n= 2460490 gmp=1470467235ns ilmp=1512510215ns k=0.972203
n= 2706544 gmp=1694172305ns ilmp=1686090604ns k=1.004793
n= 2977203 gmp=1829641942ns ilmp=1894271605ns k=0.965882
n= 3274928 gmp=2109515221ns ilmp=2132744233ns k=0.989108
n= 3602425 gmp=2422466257ns ilmp=2280946249ns k=1.062044
n= 3962672 gmp=2801679717ns ilmp=2493125764ns k=1.123762
n= 4358944 gmp=2848544434ns ilmp=2937701021ns k=0.969651
n= 4794843 gmp=3131557508ns ilmp=3264304307ns k=0.959334
n= 5274332 gmp=3501409743ns ilmp=3621818027ns k=0.966755
n= 5801770 gmp=3962427230ns ilmp=3998511783ns k=0.990976
n= 6381952 gmp=4412224075ns ilmp=4396991236ns k=1.003464
这个结果很正常，不同的机器配置不同，所以在不同机器上性能结果会有区别是正常的。这充分说明了测试的重要性 对。我也时常注意到这种情况。两个不同的实现，在测试软硬件环境A，实现1快于实现2，在另一个环境则相反。甚至在同样的硬件环境，只是测试上下文环境稍有不同，就会得到不同的结果。 这几天把除法写完了，主要是3/2试商法，分治法和牛顿迭代法。
最终速度图像如下，我的CPU是 intel i5 3320M@2.6GHz 的，L1=128k。内存 4G DDR3 @1.6GHz.
在大于百万的字长，我的库稳定快于GMP，但不多。按楼上的测试，基本上换个平台这个优势就不稳定了。
http://chuantu.biz/t5/129/1499476069x2728282818.png
在几十万字长的范围内，因为FFT类算法的通病，只支持特定字长乘法，所以一些字长范围都要向上ceiling到最小可行字长。
所以可以看到一系列的平台。
同时我的除法和GMP的除法在很大一段范围内互相交错……这是因为平台范围刚好错开。
http://chuantu.biz/t5/129/1499476126x2165260516.png
在几万的范围，我的除法稳定慢于GMP。乘法基本没有差别。
http://chuantu.biz/t5/129/1499476176x2165260516.png
几千。
同样是除法明显慢，乘法没有区别。
http://chuantu.biz/t5/129/1499476197x2165260516.png
几百，Toom-Cook范围。
完全慢于GMP。
我认为这是可以接受的，考虑到GMP实现了13种:LToom-Cook算法来加速这一段，而我非常懒只写了4个（而且未来也不会再写更多，我就是这么懒:sleepy:）。
http://chuantu.biz/t5/129/1499476210x2165260516.png
几十，汇编范围。
我的汇编技术基本上是向GMP学习的，比如那个3/2试商法，还有各种汇编过程中的四路循环展开以及多通道跳入等等，所以速度差别不大。
它的代码充满了各种诡异迭代而难以理解的宏，而且还是用AT&T语法写的，简直是一场灾难。我花了很多个晚上才理顺了那些玩意，然后用intel语法重写。
（还发现了它乘法 mul_basecase 以及平方 sqr_basecase 汇编里的一个小bug，读内存越界。因为没有越界写，所以如果没有当场访问冲突崩溃，就不会造成什么后果。）
http://chuantu.biz/t5/129/1499476221x2165260516.png
总之，除了乘法的SSA/FFT算法是完全由我自己写的，其他代码都或多或少借鉴了GMP。
一个原则是，一切算法必须经过我亲手证明其正确性，才能把它写进ilmp。（所以借鉴绝不是复制粘贴代码:)）
（有些算法真的很难证明，我觉得最难的就是那个3/2试商。不知道为什么没人回复，太长懒得看吗:sleepy:）
有些实在是没法证明，那只能我自己想一个等效的算法补上去。
楼主辛苦了，性能真的很不错。 最后就剩下 Base Conversion。
内部使用$2^64$进制，使用十进制输入输出就必须进行转进制。
输出基本就是把一个高精度数不断地除以10的某个次方，余数和商就是输出的低位部分和高位部分。
输入就是把高位部分乘以10的某个次方的积加上低位部分。
由于除法总是慢于乘法，所以输出总是慢于输入。
我对输出做了一点点优化：除法就是乘以逆（倒数），所以我提前算了所有10的高次幂的逆，这样就把除法转换成了乘法，此时理论输出耗时是输入的两倍（需要两次乘法，一次把被除数乘以逆得到商，一次把商乘以除数（然后和被除数做差）来得到余数）。
所以可以看到在我的输入和GMP差不多快的情况下，大字长的输出比GMP快了10%~15%，而且基本就是两倍的输入耗时。
http://chuantu.biz/t5/132/1499835620x2728282818.png
我还做了n=100万以上的测试，对这样大的n，我的输出快于GMP20%~30%。
因为过于耗时，点没取得很密，画不出来图，列个表：
n ilmp_in gmp_inilmp_outgmp_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%。
这应该是最后四张速度测试图。
ilmp 的核心到这里就结束了。
也许以后会写 gcd，divexact 什么的新函数，但现在我可以休息了。
乘，除，平方，开根：
http://chuantu.biz/t5/147/1500865595x2794031223.png
新写了开根的函数。测速过程是求一个 2n-limb 整数 x 的平方根 s 以及余数 r，满足 x = s*s + r 。
开根使用的是 Karatsuba 算法。
对数比例图：
http://chuantu.biz/t5/147/1500865769x2794031223.png
对 n=1 到 700万，这是各种运算耗时与一次乘法的耗时比。
可以看到在 n 很大时，一次平方大概耗时相当于 0.67 次乘法，除法是 2.8，平方根是 3.
进制转换输入输出：
http://chuantu.biz/t5/147/1500865429x2728281592.png
该说的楼上都说过了，花了很长时间测到了700万。
求2的平方根（不输出）：
http://chuantu.biz/t5/147/1500865904x2794031223.png
求2的近似平方根时，舍入可能是 floor 也可能是 round，保证最后一位与精确值不会差 1 以上即可。当然此时也不用求余数了。
gmp 使用的算法仍然是 Karatsuba，将 2 后面补上 2n 个 0，然后看作一个大整数去求它的根。
我使用的算法是 牛顿迭代法求逆平方根（1/√x），再乘原数得到平方根 （简称牛顿法）。
显然，如果原数的精度很低（比如 2 只有 1 位精度），而我们要求很高精度的根，那牛顿法是远远好于Karatsuba方法的，实测加速达到了 100% 以上（如图）。
进一步测速表明，当根的精度大约是原数精度的十倍以上时，牛顿法要更好（否则 Karatsuba 更好）。
在 ilmp 里，Karatsuba 和 牛顿法 是由 sqrt 函数自动选择的，调用者并不需要担心到底使用何种算法。