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

liangbch · 发表于 2017-7-25 11:55:19

关于求平方根算法。我说一点儿。求平方根的算法，大体用到2个公式，

1.公式1，\(x_1=0.5*(x_0+a/x_0)\)
2.公式2，\(e=a*x_0*x_0-1, x_1=x_0-0.5*x_0*e， x0,x1\)为平方根的倒数。

Paul Zimmermann 的《Karatsuba Square Root》文中的方法，使用递归式的分治法，直接求平方根，在计算过程中用到除法，其基本原理源于公式1。

第二种方法，使用迭代公式2，逐步求精，求出平方根的倒数到指定精度，然后再做一次满精度的大数乘法。这个方法不需要使用除法。

Paul Zimmermann 的《Karatsuba Square Root》一文引用 BRENT, R. P的结论。
提到当前渐进复杂度最低的算法是使用以FFT乘法为基础的牛顿迭代法，其复杂度为5倍的M(n),M(n)是两个n位数乘法的复杂度。

再看Richard P Brent的那篇曾经有10个版本之多的那篇著名的论文《MULTIPLE-PRECISION ZERO-FINDING METHODS AND THE
COMPLEXITY OF ELEMENTARY FUNCTION EVALUATION》，文中说道，求平方根的倒数的复杂度为4.5M(n),求平方根还要另外一次乘法，这样总的复杂度是5.5M(n),稍高于Paul Zimmermann提到的5倍M(n).

这篇文章还提到，使用第一种迭代公式的算法，计算平方根需要8倍或者6倍的M(P), 比第二种迭代公式更慢。
我没有看懂其复杂度推导的过程。但据我自己的分析，在实践中，复杂度应该没有这么高。楼主的结果有力的支持了我的猜测。楼主可否给出更详细的数据，在不同规模下，两种开方算法的性能与同等规模乘法的性能比较。即。求2n位数(低n位为0）的平方根
和求两个n位数的乘法的速度之比。

补充一点，测试结果应该区分2种情况，第一种，只求平方根，不需要余数。第二种，需要平方根和余数。

Ickiverar · 发表于 2017-9-1 07:50:47

本帖最后由 Ickiverar 于 2017-9-1 08:20 编辑

你所说的，在不同规模下，两种开方算法求2n位数(低n位为0）的平方根与同等规模乘法比较，我可以先给你一个理论上的答案。实测如果我有时间再写代码吧。
记 Karatsuba 算法求 2n 位数平方根的用时为 K(n) ，牛顿法用时为 T(n)，乘法为 M(n)=O(n) ( 假设在 FFT 规模 )，2n 位数除以 n 位数获得商和余数用时为 D(n)，平方为 S(n)。

1. 平方复杂度分析。
因为乘法的算法是FFT，平方只需要做一次 FFT 和一次 IFFT 共两次；乘法需要三次，用时比 2/3。FFT 后的逐点系数相乘用时比应该也是一样的，因为平方的递归还用平方，乘法递归还用乘法。
故 S(n) = 2/3 M(n)。

2. 除法复杂度分析。
A/B 除法主要使用牛顿迭代法求逆 V=B^-1，然后 V*A 获得近似商 Q，然后 B*Q 获得余数 R=A-BQ。
求逆是牛顿迭代 V=V+V(1-BV) 计算的。假设我们已经有了 n/2 位精度的 V，要求 V 的 n 位精度。迭代需要首先计算 n 位的 B 与 n/2 位的 V 的积。结果虽有 3n/2 位，但高 n/2 位一定很接近 1，我们就可以模 Base^n±1 做这个乘法来获得其低 n 位。乘法的 SSA 算法本来就是用来计算模 Base^n±1 的乘法的。这个乘法用时相当于两个 n/2 位数相乘（SSA用时与结果长度成正比），即 M(n/2)。之后，V*(1-BV) 只需要计算 n/2 位精度的结果，用时 M(n/2)。
这样，求 n 位逆的复杂度 P(n) = P(n/2) + 2*M(n/2)， P(n)=2M(n).
我们一般只计算 n/2 位的 V，那么 P(n/2)=M(n). 然后计算 V*A 获得近似商 Q 的高 n/2 位（M(n/2)），然后用 Q 的高半部分乘以 n 位的 B ，结果与 A 做差获得部分余数（由于积的高 n/2 位和 A 一定一样，所以仍然使用模 Base^n±1 的乘法获得低 n 位即可，M(n/2))。对部分余数 R' 重复 Q=V*R' 及 R=R'-BQ，这样 D(n)=P(n/2)+(M(n/2)+M(n/2))*2=3M(n).
至于为什么只计算 n/2 位的V……如果我们计算 n 位的 V ，复杂度就要增加到 3.5M(n) ，证明留做习题。

3.Karatsuba算法复杂度分析
假设我们要计算 2n 位的 A 的 n 位平方根，则首先将 A 分成四份，a3*b^3+a2*b^2+a1*b+a0，b=Base^(n/2)，a0，a1，a2，a3 各有 n/2 位。
首先迭代计算 n 位数（a3*b+a2）的 n/2 位根 s1 和余数 r1，用时 K(n/2)。
然后计算（r1*b + a1）除以（2*s1）的商 q 和余数 u，除法用时 3*M(n/2)。
A的根大约就是 s1*b+q 。
A开根的余数 r = u*b + a0 - q^2，平方耗时 2/3 M(n/2)=M(n)/3。
故 K(n)=K(n/2)+3M(n)/2+M(n)/3=3.67M(n).
如果不需要余数，最后一次平方的 M(n)/3 有可能可以省去，复杂度就是 3.33 M(n)。

4.Newton开根法复杂度分析
牛顿法我们首先算 n 位的逆 V 。假设我们已经算好了 n/2 位。
迭代是 V=V+(1-AV^2)V/2，首先计算V^2，S(n/2)=M(n)/3；然后计算AV^2，仍然使用模乘法，复杂度是 M(3n/4)；然后和1做差，结果乘 V，复杂度 M(n/2)。
综上求逆复杂度 U(n)=U(n/2)+M(n)/3+3M(n)/4+M(n)/2=3.33M(n).
然后乘A得结果，M(n)，那就 T(n)=U(n)+M(n)=4.33M(n).
如果还要算余数，那更糟糕，得把结果平方做差，5M(n).
综上，Newton开根法确实远远不如Karatsuba算法。
但是！但是！
如果 A 的精度基本就是个 0 ，比如求 2 的高精度开根时， A=2，只有 1 位精度。
牛顿开根不算余数的复杂度就是（证明留做习题）
T(n)=T(n/2)+M(n)/3+M(n)/2=1.67M(n)
远远好于Karatsuba！

所以我才要实现两种算法，根据要求的精度选择具体使用哪个。

综上：

平方 S(n)= 0.67 M(n)
除法 D(n)= 3.00 M(n)

Karatsuba平方根：
要余数 K(n)= 3.67 M(n)
不要余数 K(n)= 3.33 M(n)

Newton迭代平方根：
要余数 T(n)= 5.00 M(n)
不要余数 T(n)= 4.33 M(n)
被开根数极短且不要余数 T(n)=1.67M(n)

实际测量中，对于除法和开根，复杂度会更低。
因为事实上，M(n/2) 略微比 M(n)/2 要小一些，因为 M(n)=O(n ln n) 而不是简单的 O(n)，所以我们的估计总是在不断放大系数。

liangbch · 发表于 2017-9-1 13:23:09

对于平方根运算，这里也顺便分享一下我的成果。

1.平方根依赖大数乘法，大数平方，和大数除法，所以这些基础运算的性能对平方根运算的性能都有影响。
2.大数乘法。我的大数乘法实现了basecase乘(使用了大量的汇编代码），和9种Toom-cook乘法，没有实现更高级的FFT，SSA，NTT（数论变换)算法。
3.大数平方。我没有对平方运算做特别的实现，平方运算仍然调用乘法函数。
4.大数除法。我实现了basecase除(主要用汇编实现)，以及通过使用牛顿法求除数倒数的方法。除数倒数的长度并非总是除数长度的一半，而是根据实际评测结果，根据除数长度的不同，除数倒数的长度取 dn/2+1,dn/3+1等，这里dn为除数的长度。我没有实现分治法的大数除法。
5. 和GMP不同，我实现了求平方根的函数的两种接口。第一种接口的函数包含“sqrtrem”字样，这类函数的功能和GMP函数mpn_sqrtrem相同，可以同时求出平方根和余数。另一类函数仅包含“sqrt”字样，这个接口的函数仅用来求平方根，不求余数，其结果的误差在绝大多数情况下不超过+-0.5，在极少数情况下，误差可超过+-0.5，但不超过+-1.0，对于每类函数，我都实现了两个算法。第一个算法，首先使用牛顿法求出平方根的倒数，再做一次乘法，求出平方根，这个版本的性能比我想象中的要差很多。第二个算法使用分治法（karatsuba）求平方根，这个版本给我带来很大的惊喜，在许多情况下，都比GMP快。

下面四个函数中，以_kara结尾的函数表示使用karatsuba算法，以__NR结尾的函数表示使用牛顿迭代法。
“_bin30_sqrt_kara” ，“_bin30_sqrt_NR”，“_bin30_sqrtrem_kara”，“_bin30_sqrtrem_NR” 。

性能测试结果。
下面的表格中，len表示被开方数的长度，单位为limb,一个limb对应于30比特（而楼主的大数库中，一个limb对应于64比特）。在和GMP的mpn_sqrtrem性能比较中，互有胜负，对于_bin30_sqrt_kara函数,当len为150-39822时，_bin30_sqrt_kara函数胜出，当len等于1500左右时，_bin30_sqrt_kara函数的性能可达到mpn_sqrtrem的180%左右。
对于_bin30_sqrtrem_kara函数，在len在300到9848之间时，_bin30_sqrtrem_kara胜出，在len为2200到3400之间时，_bin30_sqrtrem_kara函数的性能可达mpn_sqrtrem的130% 左右。

测试数据说明。
len列表示被开方数长度，一个单位对应于30个比特。
gmp_mpn_sqrtrem列表示gmp mpn_sqrtrem函数运行时间和我的大数乘法函数mpn_b30_mul的运行时间之比。数字越大，表示越慢，这里以mpn_b30_mul函数为基准。
bcmp_sqrt_kara列表示_bin30_sqrt_kara函数的运行时间和mpn_b30_mul函数的运行时间之比。数字越大，表示越慢。
其余的3列分别表示其他的三个函数的运行时间，仍以mpn_b30_mul函数为基准.

liangbch · 发表于 2017-9-1 13:30:42

len gmp_mpn_sqrtrem bcmp_sqrt_kara bcmp_sqrt_NR bcmp_sqrtrem_kara bcmp_sqrtrem_NR
2 5.711 3.850 3.843 5.617 5.611
4 8.241 13.325 21.785 12.009 25.906
6 12.244 18.145 18.573 16.318 23.133
8 7.907 11.868 19.205 10.937 22.523
10 10.017 11.904 13.980 11.612 16.649
12 8.820 9.687 12.779 9.445 15.481
14 8.604 8.672 13.216 8.459 15.662
16 7.519 9.129 12.817 9.376 15.056
18 10.006 8.999 11.933 8.893 13.643
20 9.573 8.609 11.162 9.939 13.268
22 8.738 9.005 10.460 9.164 12.531
24 7.601 7.886 9.897 8.133 11.841
26 8.210 7.353 9.354 7.582 11.186
28 7.429 7.107 9.891 7.688 11.734
30 6.513 7.160 8.969 7.462 10.739
32 5.959 6.048 9.216 6.200 10.906
34 4.918 6.377 8.528 6.699 10.160
36 6.584 5.925 8.199 6.537 9.801
38 5.611 5.413 7.540 5.674 9.126
40 5.590 5.673 7.387 6.186 8.883
42 5.107 5.015 6.984 5.331 8.553
44 5.088 5.270 6.856 5.577 8.271
46 4.742 5.058 6.558 5.350 8.793
48 4.560 4.900 6.441 5.182 7.810
50 4.174 4.606 6.247 4.794 7.611
52 4.601 4.170 6.054 4.518 7.412
54 4.178 4.592 6.079 4.726 7.395
56 4.118 4.401 6.119 4.639 7.433
58 3.856 4.149 5.926 4.415 7.221
60 3.906 4.104 5.891 4.511 7.178
62 3.713 4.303 5.765 4.580 7.038
64 3.363 3.671 5.757 3.860 7.056
66 3.373 4.113 5.585 4.423 6.825
68 3.213 4.209 5.300 4.416 6.512
70 3.749 3.870 5.361 4.116 6.527
72 3.410 3.576 5.330 3.809 6.540
74 3.437 3.630 5.145 3.971 6.275
76 3.248 3.603 5.138 3.960 6.285
78 3.243 3.454 5.022 3.861 6.149
80 3.110 3.187 5.039 3.509 6.179
82 3.148 3.418 4.871 3.645 5.957
84 3.106 3.160 4.845 3.530 5.965
86 3.086 3.108 4.687 3.476 5.782
88 3.043 2.995 4.722 3.319 5.835
90 2.945 2.823 4.516 3.206 5.812
92 2.874 2.916 4.556 3.187 5.644
94 2.967 3.067 4.655 3.203 5.780
96 2.714 2.902 4.542 3.204 5.662
98 2.644 2.738 4.373 3.078 5.512
100 2.702 2.601 4.582 2.981 5.510
102 2.596 2.798 4.535 3.042 5.692
104 2.756 2.838 4.603 3.165 5.700
106 2.525 2.608 4.415 2.909 5.515
108 2.625 2.558 4.476 2.851 5.598
110 2.559 2.631 4.389 2.880 5.512
112 2.571 2.529 4.522 2.800 5.646
114 2.442 2.448 4.394 2.731 5.487
116 2.448 2.476 4.303 2.771 5.388
118 2.366 2.321 4.227 2.604 5.310
120 2.450 2.356 4.275 2.615 5.392
122 2.374 2.237 4.183 2.506 5.281
124 2.420 2.445 4.239 2.715 5.320
126 2.337 2.230 4.149 2.510 5.232
128 2.216 2.467 4.306 2.710 5.413
130 2.195 2.446 4.034 2.686 5.106
132 2.194 2.324 4.054 2.574 5.107
134 2.062 2.247 3.882 2.607 4.942
136 2.062 2.273 3.982 2.520 5.036
138 2.135 2.115 3.814 2.412 4.904
140 2.167 2.306 3.554 2.482 4.554
142 2.259 2.302 3.889 2.529 4.897
144 2.322 2.370 4.041 2.703 5.079
146 2.315 2.342 3.994 2.615 5.042
148 2.257 2.289 3.794 2.565 4.781
150 2.208 2.132 3.757 2.511 4.744
152 2.306 2.237 3.786 2.486 4.768
154 2.243 2.183 3.848 2.703 4.779
156 2.268 2.221 3.790 2.505 4.772
158 2.206 2.093 3.811 2.398 4.782
160 2.203 2.157 3.868 2.460 4.855
162 2.167 2.114 3.854 2.383 4.810
164 2.237 2.152 3.786 2.480 4.627
166 2.197 2.046 3.735 2.430 4.730
168 2.204 2.113 4.107 2.833 5.046
170 2.729 2.006 3.716 2.281 4.678
172 2.268 2.073 3.745 2.381 4.748
174 2.196 2.021 3.781 2.488 4.794
176 2.217 2.037 3.820 2.352 4.817
178 2.127 1.982 3.779 2.334 4.895
180 2.149 1.986 3.621 2.325 4.558
182 2.056 1.890 3.480 2.237 4.417
184 2.050 1.953 3.455 2.184 4.371
186 1.982 1.847 3.453 2.377 4.350
188 2.089 1.907 3.616 2.155 4.484
190 2.121 1.870 3.552 2.150 4.466
192 2.005 1.892 3.708 2.157 4.611
194 2.093 1.837 3.551 2.102 4.452
196 1.993 1.793 3.602 2.080 4.489
198 2.070 1.814 3.548 2.096 4.447
200 2.074 1.853 3.714 2.126 4.613
202 2.010 1.785 3.528 2.061 4.423
204 2.031 1.815 3.616 2.088 4.520
206 2.096 1.801 3.552 2.346 4.471
208 2.083 1.847 3.725 2.097 4.624
210 2.042 1.819 3.538 2.200 4.442
212 1.887 1.758 3.584 2.037 4.468
214 1.913 1.809 3.547 2.330 4.472
216 1.925 1.827 3.711 2.093 4.630
218 1.886 1.740 3.552 2.136 4.463
220 1.884 1.754 3.599 2.003 4.532
222 1.942 1.771 3.575 2.034 4.510
224 1.876 1.765 3.742 2.026 4.669
226 1.859 1.722 3.565 2.031 4.492
228 1.871 1.664 3.585 1.929 4.481
230 1.846 1.697 3.565 1.966 4.501
232 1.911 1.725 3.687 2.032 4.609
234 1.885 1.648 3.551 1.913 4.503
236 1.919 1.775 3.641 2.039 4.581
238 1.856 1.668 3.594 1.938 4.525
240 1.895 1.693 3.687 1.990 4.651
242 1.884 1.669 3.570 1.942 4.516
244 1.918 1.650 3.606 1.948 4.540
246 1.886 1.687 3.557 1.965 4.513
248 1.909 1.668 3.721 1.974 4.670
250 1.863 1.634 3.567 1.968 4.518
252 1.827 1.607 3.613 1.938 4.561
254 1.857 1.585 3.590 1.879 4.525
256 1.812 1.692 3.900 2.272 4.939
258 2.325 1.889 3.706 2.186 4.666
260 2.330 1.902 3.690 2.163 4.658
262 1.846 1.638 3.547 1.918 4.497
264 1.865 1.635 3.794 1.886 4.631
266 1.772 1.619 3.547 1.926 4.376
268 1.705 1.553 3.465 2.417 5.811
270 1.754 1.221 2.674 1.447 3.410
272 1.761 1.577 3.770 2.069 4.780
274 2.217 1.773 3.677 2.025 4.492
276 1.779 1.553 3.585 1.837 4.562
278 1.780 1.586 3.711 1.881 4.633
280 1.784 1.506 3.295 1.790 4.229
282 1.900 1.593 3.489 1.900 4.505
284 2.037 1.689 3.690 1.999 4.739
286 2.020 1.680 3.578 1.986 4.610
288 1.955 1.706 4.053 2.177 5.168
290 2.558 1.657 3.716 2.083 4.428
292 1.912 1.695 3.608 1.999 4.776
294 2.068 1.715 3.655 2.020 4.683
296 2.044 1.727 3.779 2.065 4.798
298 1.984 1.700 3.665 1.987 4.692
300 1.983 1.681 3.588 1.977 4.599
302 1.989 1.629 3.643 1.940 4.673
304 1.963 1.640 4.026 2.083 5.132
306 1.925 1.583 3.567 1.896 4.606
308 2.037 1.576 3.611 1.871 4.604
310 2.075 1.603 3.617 1.890 4.640
312 2.042 1.616 3.724 1.927 4.769
314 2.031 1.613 3.647 1.953 4.686
316 1.986 1.624 3.561 1.916 4.611
318 1.995 1.582 3.620 1.923 4.675
320 1.979 1.581 3.670 1.917 4.656
322 1.947 1.701 3.593 1.930 4.631
324 2.053 1.692 3.683 2.030 4.797
326 2.085 1.642 3.698 1.967 4.680
328 2.023 1.614 3.673 1.878 4.640
330 2.061 1.666 3.665 1.959 4.718
332 1.972 1.610 3.595 1.945 4.655
334 2.013 1.592 3.585 1.891 4.664
336 2.105 1.621 3.871 1.935 4.909
338 1.955 1.532 3.571 1.823 4.623
340 2.023 1.550 3.551 1.845 4.597
342 2.033 1.581 3.630 1.893 4.671
344 2.054 1.560 3.686 1.857 4.745
346 2.016 1.624 3.821 2.173 4.919
348 2.876 1.816 3.670 1.833 4.508
350 2.164 1.543 3.660 1.878 4.699
352 2.054 1.582 3.842 1.868 4.808
354 3.698 1.622 3.813 2.470 6.195
356 2.053 1.604 3.563 1.966 5.200
358 1.474 1.166 2.729 1.428 3.683
360 2.079 1.491 3.569 1.814 4.620
362 2.029 1.604 3.665 1.929 4.731
364 2.066 1.551 3.599 1.949 4.867
366 2.814 1.713 3.697 1.802 4.514
368 2.200 1.590 3.694 1.878 4.763
370 2.047 1.497 3.597 1.822 4.657
372 2.426 1.551 3.574 1.874 4.643
374 2.005 1.498 3.616 1.811 4.666
376 2.076 1.528 3.698 1.837 4.765
378 2.029 1.509 3.641 1.811 4.722
380 2.092 1.605 3.689 1.892 4.802
382 2.070 1.504 3.618 1.807 4.687
384 2.018 1.526 3.707 1.862 4.767
386 2.033 1.490 3.721 1.973 4.909
388 2.022 1.525 3.567 1.837 4.631
390 2.163 1.545 3.733 1.799 4.692
392 2.080 1.504 3.683 1.813 4.758
394 2.054 1.477 3.643 1.790 4.721
396 2.080 1.510 3.558 1.821 4.635
398 2.099 1.479 3.604 1.779 4.666
400 2.166 1.491 3.850 1.845 4.803
402 2.175 1.520 3.606 1.884 4.853
404 2.069 1.493 3.669 1.819 4.689
406 2.024 1.417 3.482 1.742 4.525
408 2.084 1.479 3.846 2.775 4.867
410 2.121 1.463 3.169 1.551 4.122
412 2.081 1.481 3.581 1.860 4.803
414 2.167 1.435 3.615 1.753 4.695
416 2.083 1.438 3.813 1.940 4.990
418 2.059 1.496 3.779 2.378 6.030
420 2.074 1.497 3.587 1.807 4.661
422 2.006 1.429 3.609 1.757 4.690
424 1.957 1.453 3.575 1.744 4.654
426 1.823 1.134 2.712 1.365 3.513
428 1.866 1.417 3.705 1.914 4.823
430 2.014 1.378 3.539 1.696 4.585
432 1.838 1.424 3.639 1.743 4.706
434 1.857 1.401 3.804 1.873 4.882
436 1.795 1.443 3.466 1.744 4.518
438 1.835 1.396 3.619 1.737 4.697
440 1.818 1.407 3.591 1.709 4.719
442 1.806 1.417 3.601 1.731 4.687
444 1.849 1.384 3.546 1.697 4.635
446 1.924 1.433 3.469 1.700 4.525
448 1.908 1.409 3.685 1.719 4.864
450 1.809 1.425 3.623 1.689 4.711
452 1.842 1.410 3.554 1.718 4.660
454 1.809 1.400 3.701 1.753 4.812
456 1.828 1.366 3.693 1.728 4.804
458 1.860 1.405 3.755 1.689 4.871
460 2.171 1.419 3.442 1.682 4.512
462 1.847 1.407 3.605 1.720 4.704
464 1.857 1.418 3.634 1.739 4.709
466 1.847 1.362 3.600 1.666 4.689
468 1.841 1.371 3.671 1.772 4.944
470 1.875 1.411 3.484 1.950 5.348
472 1.831 1.385 3.621 1.680 4.733
474 1.808 1.428 3.738 1.685 4.719
476 1.796 1.365 3.530 1.674 4.694
478 1.888 1.358 3.614 1.693 4.693
480 1.794 1.396 3.670 1.702 4.738
482 1.860 1.313 3.538 1.632 4.601
484 1.835 1.328 3.555 1.661 4.638
486 1.872 1.348 3.640 1.609 4.573
488 1.845 1.340 3.746 1.714 5.018
490 1.968 1.355 3.480 1.622 4.537
492 1.891 1.366 4.143 1.939 4.827
494 1.851 1.351 3.585 1.674 4.666
496 1.891 1.380 3.813 1.820 4.967
498 2.377 1.302 3.554 2.148 5.095
500 1.849 1.224 3.215 1.506 4.205
502 1.825 1.341 3.605 1.652 4.695
504 1.850 1.348 3.611 1.652 4.700
506 1.835 1.300 3.586 1.620 4.677
508 1.838 1.326 3.527 1.643 4.622
510 1.816 1.348 3.631 1.664 4.736
512 1.833 1.344 3.744 1.667 4.912
512 2.307 1.453 3.693 1.634 4.832
528 1.819 1.353 3.702 1.819 4.952
546 1.913 1.407 3.798 1.729 4.772
564 1.944 1.349 3.590 1.689 4.613
584 2.211 1.465 4.013 1.752 4.941
602 2.185 1.364 3.588 1.665 4.955
622 2.121 1.209 3.567 1.563 4.487
642 2.345 1.417 3.927 1.778 5.092
664 2.063 1.297 3.366 1.505 4.314
686 2.201 1.284 3.588 1.596 4.484
708 2.463 1.296 3.656 1.619 4.687
732 2.124 1.276 3.623 1.591 4.652
756 2.337 1.302 3.572 1.624 4.611
782 2.167 1.226 3.571 1.533 4.609
806 2.493 1.287 3.615 1.616 4.627
834 2.281 1.214 3.568 1.506 4.494
862 1.805 1.212 3.518 1.564 4.562
890 1.843 1.209 3.472 1.512 4.482
918 1.877 1.237 3.573 1.590 4.583
950 2.027 1.229 3.571 1.560 4.587
980 2.069 1.221 3.457 1.514 4.445
1012 2.014 1.257 3.527 1.648 4.872
1046 2.070 1.638 3.464 1.709 4.446
1080 2.081 1.415 3.621 1.674 4.638
1116 2.103 1.448 3.556 1.714 4.558
1154 2.326 1.432 3.591 1.843 4.725
1192 2.208 1.431 3.649 1.835 4.682
1232 2.261 1.230 3.135 1.573 4.033
1272 2.490 1.449 3.721 1.866 4.792
1314 2.602 1.384 3.555 1.795 4.566
1358 2.392 1.569 3.872 2.069 4.922
1402 2.555 1.504 3.750 1.937 4.797
1448 2.610 1.352 3.440 1.784 4.415
1496 2.694 1.475 3.764 1.919 4.799
1546 2.849 1.473 3.735 1.861 4.621
1596 2.866 1.627 3.789 2.122 4.861
1650 2.880 1.612 3.742 2.054 4.722
1704 2.263 1.477 3.663 1.907 4.703
1760 2.219 1.532 3.867 2.024 4.939
1818 2.236 1.578 3.725 1.973 4.653
1878 2.315 1.549 3.769 1.896 4.945
1940 2.227 1.479 3.694 1.909 4.725
2004 2.268 1.500 3.699 1.904 4.651
2070 2.215 1.497 3.744 2.000 4.820
2138 2.375 1.498 3.752 1.891 4.801
2210 2.461 1.708 4.185 2.147 4.898
2282 2.577 1.588 3.793 1.966 4.863
2358 2.410 1.531 3.665 1.964 4.671
2436 2.491 1.547 3.708 2.001 4.749
2516 2.460 1.568 3.918 1.945 4.761
2600 2.450 1.634 3.754 1.984 4.947
2684 2.576 1.527 3.604 1.899 4.607
2774 2.614 1.637 3.808 2.119 4.773
2866 2.527 1.575 3.673 2.138 4.714
2960 2.518 1.569 3.691 1.929 4.755
3058 2.582 1.607 3.822 2.006 4.896
3158 2.516 1.581 3.668 1.864 4.681
3262 2.627 1.637 3.750 1.988 4.796
3370 2.589 1.589 3.700 1.901 4.741
3482 2.151 1.587 3.749 1.975 4.789
3596 2.166 1.595 3.727 1.923 4.783
3716 2.170 1.554 3.712 1.869 4.794
3838 2.121 1.518 3.654 1.887 4.685
3966 2.231 1.544 3.727 1.857 4.828
4096 2.256 1.631 3.938 2.014 5.020
4232 2.186 1.502 3.620 1.902 4.706
4372 2.095 1.524 3.653 1.817 4.682
4516 2.441 1.691 3.842 2.036 4.997
4664 2.369 1.625 3.799 1.981 4.858
4818 2.590 1.764 3.940 2.185 5.004
4978 2.554 1.747 3.910 2.104 4.993
5142 2.219 1.511 3.550 1.963 4.482
5312 2.457 1.701 3.837 2.062 4.860
5488 2.513 1.769 3.935 2.128 4.998
5668 2.439 1.724 3.880 2.112 4.776
5856 2.382 1.723 3.875 2.156 6.110
6050 2.401 1.698 3.830 2.095 4.783
6248 2.512 1.744 3.814 2.116 4.808
6456 2.483 1.710 3.849 2.092 4.887
6668 2.532 1.716 3.775 2.099 4.772
6888 2.186 1.687 3.766 2.059 4.766
7116 2.139 1.692 3.894 2.105 4.883
7352 2.284 1.725 3.971 2.113 5.060
7594 2.196 1.665 3.797 2.030 4.861
7844 2.234 1.673 3.826 2.060 4.926
8104 2.225 1.644 3.738 2.008 4.754
8372 2.236 1.616 3.839 2.030 4.856
8648 2.196 1.729 3.953 2.118 5.048
8934 2.144 1.700 3.819 2.133 4.822
9228 2.128 1.738 3.842 2.114 4.881
9534 2.138 1.726 3.880 2.121 4.920
9848 2.148 1.741 3.856 2.146 4.891
10174 2.113 1.739 3.850 2.118 4.850
10510 2.131 1.750 3.852 2.191 4.866
10856 2.096 1.761 3.912 2.144 4.898
11214 1.554 1.281 2.846 1.558 3.640
11586 2.135 1.758 3.867 2.145 4.918
11968 2.049 1.743 3.801 2.109 4.841
12364 2.058 1.762 3.762 2.134 4.786
12772 2.095 1.782 3.728 2.138 4.763
13194 2.144 1.798 3.741 2.150 4.909
13628 2.012 1.754 3.792 2.126 4.860
14078 1.983 1.725 3.764 2.110 4.799
14544 2.170 1.896 4.193 2.132 4.865
15024 1.978 1.692 3.790 2.041 4.747
15520 2.006 1.690 3.764 2.070 4.813
16032 1.953 1.700 3.817 2.075 4.798
16562 1.993 1.739 3.858 2.088 4.936
17110 2.027 1.662 3.766 2.031 4.829
17674 1.954 1.663 3.746 2.002 4.783
18258 2.079 1.804 3.840 2.190 4.874
18862 2.104 1.778 3.840 2.153 4.856
19484 2.047 1.775 3.844 2.144 4.846
20128 2.061 1.804 3.889 2.185 4.968
20792 2.144 1.794 3.874 2.181 4.901
21478 2.065 1.775 3.841 2.160 4.846
22188 2.024 1.795 3.861 2.186 4.876
22920 1.954 1.788 3.864 2.160 4.873
23678 1.972 1.801 3.897 2.210 4.936
24460 1.938 1.825 3.844 2.225 4.901
25268 2.011 1.870 3.965 2.297 5.003
26102 1.961 1.869 3.894 2.240 4.931
26964 1.868 1.784 3.848 2.175 4.846
27854 1.889 1.798 3.947 2.240 4.966
28774 1.840 1.792 3.985 2.251 4.951
29724 1.859 1.792 3.844 2.137 4.833
30706 1.831 1.787 3.842 2.223 4.881
31720 1.776 1.756 3.887 2.145 4.955
32768 1.835 1.956 3.918 2.119 5.002
33850 1.855 1.780 3.998 2.184 5.057
34968 1.953 1.839 4.069 2.214 5.090
36122 1.912 1.826 3.929 2.225 5.026
37316 1.889 1.817 3.871 2.289 4.929
38548 1.816 1.784 3.881 2.160 4.963
39822 1.799 1.783 3.884 2.167 4.944
41136 1.717 1.829 3.914 2.213 5.020
42494 1.776 1.763 3.867 2.124 4.849
43898 1.794 1.815 3.808 2.206 4.925
45348 1.693 1.791 3.828 2.161 4.868
46846 1.687 1.794 3.815 2.181 4.854
48392 1.651 1.814 3.832 2.196 4.874
49990 1.572 1.813 3.825 2.187 4.869
51642 1.685 1.837 3.831 2.249 4.895
53348 1.668 1.836 3.878 2.206 4.905
55108 1.615 1.770 3.774 2.157 4.876
56928 1.561 1.786 3.773 2.142 4.794
58810 1.521 1.759 3.785 2.124 4.774
60752 1.449 1.614 3.496 1.969 4.414
62758 1.546 1.862 3.887 2.180 4.892
64830 1.527 1.821 3.969 2.198 4.953

Ickiverar · 发表于 2017-9-1 14:39:30

liangbch 发表于 2017-9-1 13:23
对于平方根运算，这里也顺便分享一下我的成果。

1.平方根依赖大数乘法，大数平方，和大数除法，所以这些 ...

你的库和 GMP 库各是 64 位还是 32 位的？
让我首先假设你传给 GMP 库 mpn_sqrtrem 的参数长度经过了适当减少，以使得你的库和它计算的事实上是相同 bit 长度的平方根。
那，180%的效率说明你的乘法汇编（可能还有ToomCook）做得非常非常优秀。我阅读了 GMP x64 的 basecase 乘法代码（四百多行汇编），感觉它在不考虑 MMX，SSE，AVX 什么黑魔法的情况下已经是最优了（它的 basecase 比我旧版的 Toom-Cook33 都快）（我不想使用高级汇编指令加速，头大）。

对你所关心的位数，我实现了除法的 basecase，分治，牛顿法（事实上我也不总是计算 n/2 位的逆，但对 2n 除以 n 的除法一定是这样）。
我在这个帖子第一页发的图里，那个所有算法和乘法用时比较的对数图，里面的开根就是求余数的Karatsuba算法。可以看到在100~2000位间，它大概总是乘法的 1.2~1.5 倍，这个倍数比你的库似乎稍微低一些。
这意味着我的平方根算法比你的稍好（也许归因于除法的分治法），而更基础的比如说乘法比你的烂。
（2000位后我的比例猛增到2，是因为进入了 FFT 范围，乘法复杂度增速骤减，从幂函数向线性转化。）

牛顿法是真的慢，而且是死慢死慢，没救的那种。所以我基本没有想过给牛顿法做 √2n = n 的测速。没多大意义。
但我测过了 √n = 4n ~ 20n 这样的速，最终选择了大概这么一个分支： nf>=10*na+const，（na是输入精度，输出具有 na/2 的整数精度加上 nf 的小数精度），也就是对 n 位精度数要计算它的 10n 位到 11n 位以上精度根的时候，应该选择牛顿法。（const是为了避免对过小的数用牛顿法，因为牛顿法只对大数很有效。我把它随便设成 50 位）

Ickiverar · 发表于 2017-9-1 14:51:30

我的所有平方根用同一个接口。
void ilmp_sqrt_(mp_ptr dsts,mp_ptr dstr,mp_srcptr numa,mp_size_t na,mp_size_t nf);
计算 na 位整数 numa 的平方根直到 nf 位小数，结果（长度 nf+na/2+1 位）放到 dsts；
如果 dstr 不是 NULL，计算余数（长度 nf+na/2+1 位）存入 dstr，且保证根必然是floor舍入的。
如果 dstr 是 NULL，不输出余数，此时根可能采用 floor 或 round 舍入，保证最后一位误差不超过1（且保证 0xFF....FF 的根一定有形式 0xFF...FF，而不会进位成 0x100...00）。

具体内部用什么算法是自动选的，call 一些不留接口的内部过程。选择条件楼上也说了。

liangbch · 发表于 2017-9-1 21:28:58

你的库和 GMP 库各是 64 位还是 32 位的？
让我首先假设你传给 GMP 库 mpn_sqrtrem 的参数长度经过了适当减少，以使得你的库和它计算的事实上是相同 bit 长度的平方根

我的库和GMP都是32位的。在测试速度的时候，总是在同一个程序调用我的库中的函数和GMP的函数。大数使用字符串赋值，所以我的库中的函数和GMP函数使用的是同样的数据，但是GMP表示的大数比我的库的大数要短。

wayne · 发表于 2020-1-26 17:42:14

GMP 6.2.0 版本出来了, 差不多四年没更过了.

xiaozhang · 发表于 2022-3-6 14:27:51

想看楼主的ilmp库的源码

Ickiverar · 发表于 2023-3-8 20:58:17

xiaozhang 发表于 2022-3-6 14:27
想看楼主的ilmp库的源码

很久没动力继续搞了，我直接开源吧……

https://github.com/himisawww/ilmp

账号		自动登录	找回密码
密码			欢迎注册

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

评分

点评