重大消息-一种新的任意精度对数算法研制成功

wayne

liangbch 发表于 2018-1-1 18:53
关于扣除到空循环这种写法，这不是我的发明，我借鉴了 mpfr的测试代码， 请参阅 http://www.mpfr.org/mpfr- ...

1） 嗯，关于空循环， 我只是想表达一个意思，就是针对Mathematica这种函数式语言来说(相比于指令式编程)，恐怕是多余的考虑，我们只需尽管套用AbsoluteTiming函数就足够了，无需减去循环带来的上下文开销。（反而 当用在较高的测试次数的场合下，比如10^8次，其影响是突出的，喧宾夺主了。）
2）

我猜想mathematica把log(97)当做一个常数，预存了它，调用log(97)时，直接取出来给你

我不认为存在某种语言[包括Mathematica]，有能力将 log(97) 的任意精度的值能存起来，我们输入一个位数，它就能直接返回给我们预期的精度，这个是难以想象的。
（我当然知道，很多语言会把常用的函数值预先计算出来，但这恐怕只是标准的浮点精度，跟咱们这里讨论的任意精度[or高精度]恐怕不是一回事）
3） 关于你选择的 $\sqrt{3}-1$, 没有本质差别的。 因为你最终经过了N运算，即求值运算，变成了带有特定精度【可以查看Mathematica的NumberMarks的文档】的小数(暂且这么叫)。意思就是显式的要求Mathematica进行高精度浮点运算。但是我们其实大可圆滑一点，将$\sqrt{3}-1$ 乘以$10^10000$，再取其整数部分。然而咱们这里的测试应该不用涉及到这么大的整数吧。
一句话，直接用N[Log[97],10000]就足够了，这样Mathematica内部的算法自动选择机制就会有可能采用高精度整数运算的算法。

liangbch

Mathematica 经常做符号运算，而不是高精度数值计算，这我相信。

liangbch

1. 我不认为存在某种语言[包括Mathematica]，有能力将 log(97) 的任意精度的值能存起来，我们输入一个位数，它就能直接返回给我们预期的精度，这个是难以想象的。

复制代码

任意精度是不可能，一定精度是可能的，比如缓冲100个1万位精度的常数，开销并非不可接受。如果，用户调用某个函数，当精度低于特定位，直接给出，否则重新计算。

liangbch

这个应该不是问题的聚焦点吧．　我甚至觉得Mathematica这种符号计算系统，没有理由在缓存高精度常数做大量的文章

对于常用的常数，Mathematica应该会缓冲，比如圆周率,log(2), 至于会存储多少位，我没有调查过。之所以用圆周率和log(2)做例子，是因为如果使用AGM算法计算对数的话，会用到这两个常数。另外，如果mathematica会存储刚刚计算过的某个函数的值的话？那我写的那个测试程序就没有办法应对了，因为，这个程序总是计算重复许多次的总的时间花费。

wayne

liangbch 发表于 2018-1-2 11:17
对于常用的常数，Mathematica应该会缓冲，比如圆周率,log(2), 至于会存储多少位，我没有调查过。之所以用 ...

原则上，每计算一次，应该进行ClearSystemCache[] 操作．

liangbch

ClearSystemCache也是需要开销的。
我也不至于这么不相信Mathematica.
另外，缓冲几个常数是可以理解的，这也是通行的做法，我的算法也会缓冲几个常数。

liangbch

在我这里，用Timing得到时间值确实是离散的，它总是15.6ms的总倍数。AbsoluteTiming不是
下面是我的代码，你可测试下。

1. Print["Mathematica:",$VersionNumber];
   
2. 
   
3. digArr={20,30,40,50,60,70,80,90,100,125,158,199,251,316,398,501,630,794,
   
4.   1000};
   
5. 
   
6. For [i=1,i <= Length[digArr],i++,
   
7.   prec=digArr[[i]];
   
8.   n=IntegerPart[N[2^24/(prec^1.5),10]];
   
9.   n=Max[n,3];
   
10.   (*Print["n=",n];*);
    
11.   x=N[Sqrt[3]-1,prec];
    
12.   y=N[Sqrt[5],prec];
    
13.   
    
14.   While [True,
    
15.         (*t=(AbsoluteTiming[Do[Log[x],{n}]][[1]]-AbsoluteTiming[Do[Null,{n}]][[1]]);*)
    
16.     t=Timing[Do[Log[x],{n}]][[1]]；
    
17.     (* convert to ms *)
    
18.     t *= 1000;
    
19.     (*Print["t=",t];    *)
    
20.     If [t>=10.0, Break[]];
    
21.     n*=2;
    
22.   ];
    
23.   Print[prec,":took ",t/n, " ms (", n, " eval in ", t, "ms)" ];
    
24. ]
    
25. 
    

复制代码

下面是测试结果

1. Mathematica:10.3
   
2. 
   
3. 20:took 0.0025782 ms (187574 eval in 483.603ms)
   
4. 30:took 0.00275021 ms (102102 eval in 280.802ms)
   
5. 40:took 0.00305806 ms (66317 eval in 202.801ms)
   
6. 50:took 0.00328748 ms (47453 eval in 156.001ms)
   
7. 60:took 0.00388944 ms (36098 eval in 140.401ms)
   
8. 70:took 0.00435666 ms (28646 eval in 124.801ms)
   
9. 80:took 0.0053229 ms (23446 eval in 124.801ms)
   
10. 90:took 0.00555757 ms (19649 eval in 109.201ms)
    
11. 100:took 0.0055791 ms (16777 eval in 93.6006ms)
    
12. 125:took 0.00779745 ms (12004 eval in 93.6006ms)
    
13. 158:took 0.00923411 ms (8447 eval in 78.0005ms)
    
14. 199:took 0.0130523 ms (5976 eval in 78.0005ms)
    
15. 251:took 0.0147938 ms (4218 eval in 62.4004ms)
    
16. 316:took 0.0208977 ms (2986 eval in 62.4004ms)
    
17. 398:took 0.0295456 ms (2112 eval in 62.4004ms)
    
18. 501:took 0.0521394 ms (1496 eval in 78.0005ms)
    
19. 630:took 0.0735854 ms (1060 eval in 78.0005ms)
    
20. 794:took 0.0833116 ms (749 eval in 62.4004ms)
    
21. 1000:took 0.117737 ms (530 eval in 62.4004ms)
    
22. 
    

复制代码

各次运行时间除以15.6， 在6位精度范围内，他确实是15.6的正倍数。

1. 总时间  总时间/15.6
   
2. 483.603        31.00019231
   
3. 280.802        18.00012821
   
4. 202.801        13.0000641
   
5. 156.001        10.0000641
   
6. 140.401        9.000064103
   
7. 124.801        8.000064103
   
8. 124.801        8.000064103
   
9. 109.201        7.000064103
   
10. 93.6006        6.000038462
    
11. 93.6006        6.000038462
    
12. 78.0005        5.000032051
    
13. 78.0005        5.000032051
    
14. 62.4004        4.000025641
    
15. 62.4004        4.000025641
    
16. 78.0005        5.000032051
    
17. 78.0005        5.000032051
    
18. 62.4004        4.000025641
    
19. 62.4004        4.000025641
    

复制代码

如果时间的分辨率是毫秒，该如何解释 124.801，109.201，93.6006这样的数值呢？ 如果正好碰到一个15.6的倍数不奇怪，但100%的15.6的倍数就难以解释了。

wayne

Mathematica:11.1.0 for Microsoft Windows (64-bit) (March 13, 2017), time unit = 1/1000

20:took 0.00358192 ms (187574 eval in t/15.6 = 43.0689ms)

30:took 0.00428493 ms (102102 eval in t/15.6 = 28.0449ms)

40:took 0.00424099 ms (66317 eval in t/15.6 = 18.0288ms)

50:took 0.00428055 ms (47453 eval in t/15.6 = 13.0208ms)

60:took 0.00476134 ms (36098 eval in t/15.6 = 11.0176ms)

70:took 0.00545451 ms (28646 eval in t/15.6 = 10.016ms)

80:took 0.0053314 ms (23446 eval in t/15.6 = 8.01282ms)

90:took 0.00636165 ms (19649 eval in t/15.6 = 8.01282ms)

100:took 0.00558801 ms (16777 eval in t/15.6 = 6.00962ms)

125:took 0.00650825 ms (12004 eval in t/15.6 = 5.00801ms)

158:took 0.00739908 ms (8447 eval in t/15.6 = 4.00641ms)

199:took 0.00784388 ms (5976 eval in t/15.6 = 3.00481ms)

251:took 0.0111131 ms (4218 eval in t/15.6 = 3.00481ms)

316:took 0.0156983 ms (2986 eval in t/15.6 = 3.00481ms)

398:took 0.0147964 ms (2112 eval in t/15.6 = 2.00321ms)

501:took 0.0313336 ms (1496 eval in t/15.6 = 3.00481ms)

630:took 0.0442217 ms (1060 eval in t/15.6 = 3.00481ms)

794:took 0.0417223 ms (749 eval in t/15.6 = 2.00321ms)

1000:took 0.0884434 ms (530 eval in t/15.6 = 3.00481ms)

liangbch

https://stackoverflow.com/questi ... to-15-ms-resolution
上文提到

The default timer resolution on Windows 7 is 15.6 milliseconds (ms). Some applications reduce this to 1 ms, which reduces the battery run time on mobile systems by as much as 25 percent.

liangbch

38楼数据比我的环境稍微大了那么一点点。 所有是时间都是15.625 毫秒的整倍数，15.625 正好等于1/64秒，即39楼链接提到的64HZ.