阶乘和圆周率

无心人

刚反复的测试，觉得结果推进到3.1415926535897，需要的精度是1.388E-14
所以，每次误差不能超过5E-15
我无论怎么测试，gcc每次浮点误差都大于1E-17，所以斯特拉公式调整频率得降低到100次1调整才行

无心人

无心人 发表于 2019-11-3 10:15
刚反复的测试，觉得结果推进到3.1415926535897，需要的精度是1.388E-14
所以，每次误差不能超过5E-15
我 ...

找到原因了，把double改long double, log10改log10l, 相应的输入输出格式化字符串%lf改%Lf
编译参数加-m128bit-long-double
后，测试精度高于5E-19，已经满足1000次修正后，计算到3.1415926535897精度大于1E-14的要求

现在重新复核上次计算到22.8万亿的过程

Ickiverar

如果我们写一个十进制定点数的乘法算法，只保留十进制前几十位数，使得从1乘到很大的N都不会产生严重误差，是不是可以避免浮点斯特令公式？

Ickiverar

我写了一个算法，可以单线程在18秒内搜索10亿的范围（测试起点 0xdeadbeefbadcab1e ， 测试长度 10^9 ），可以任意并行化，最高可以搜索到 2^64 - 2 的阶乘。
但我的电脑正在全线程计算海绵球的体积，所以这个项目得过几天开算。

Ickiverar

我先把我的代码贴出来吧。
使用到了mpfrc++，可以在 http://www.holoborodko.com/pavel/mpfr/ 下载。
用64位无符号定点数存储阶乘的十进制的开头，具体算法见 einit。
然后用emul逐项乘下一个数进行递推（可以每递推百万次用einit校正一次）。
搜索比较接近 pi/10*2^64 的值就行了。

1. #include<mpreal.h>
   
2. #pragma comment(lib,"mpfr-4.0.1.lib")
   
3. #include<iostream>
   
4. #include<intrin.h>
   
5. 
   
6. typedef mpfr::mpreal real;
   
7. 
   
8. //10^frac(log10(factorial(n)))*2^64/10
   
9. uint64_t einit(uint64_t n){
   
10.     static const real log10e=1/mpfr::log(10);
    
11.     static const real p264=mpfr::pow(2,64);
    
12.     return (uint64_t)(p264*mpfr::pow(10,mpfr::frac(mpfr::lngamma(n+1)*log10e)-1));
    
13. }
    
14. 
    
15. //2^64/10 ~< e*n/10^k < 2^64
    
16. uint64_t emul(uint64_t e,uint64_t n){
    
17.     //10^pos
    
18.     static uint64_t pow10[]={
    
19.         1,10,100,1000,10000,
    
20.         100000,1000000,10000000,100000000,1000000000,
    
21.         10000000000,100000000000,1000000000000,10000000000000,100000000000000,
    
22.         1000000000000000,10000000000000000,100000000000000000,1000000000000000000,10000000000000000000
    
23.     };
    
24.     static int shrn[]={3,6,9,13,16,19,23,26,29,33,36,39,43,46,49,53,56,59,63,66};
    
25.     //round(2^(64+shrn[pos])/10^(pos+1))
    
26.     static uint64_t inv10[]={
    
27.         0xcccccccccccccccd,0xa3d70a3d70a3d70a,0x83126e978d4fdf3b,0xd1b71758e219652c,
    
28.         0xa7c5ac471b478423,0x8637bd05af6c69b6,0xd6bf94d5e57a42bc,0xabcc77118461cefd,
    
29.         0x89705f4136b4a597,0xdbe6fecebdedd5bf,0xafebff0bcb24aaff,0x8cbccc096f5088cc,
    
30.         0xe12e13424bb40e13,0xb424dc35095cd80f,0x901d7cf73ab0acd9,0xe69594bec44de15b,
    
31.         0xb877aa3236a4b449,0x9392ee8e921d5d07,0xec1e4a7db69561a5,0xbce5086492111aeb
    
32.     };
    
33. 
    
34.     //res = e*n
    
35.     uint64_t res[2];
    
36.     res[0]=_umul128(e,n,res+1);
    
37. 
    
38.     if(res[1]==0)return res[0];
    
39. 
    
40.     //shrt = res >> shrn[pos]
    
41.     int pos=0;
    
42.     uint64_t shrt;
    
43.     if(res[1]>=pow10[19]){
    
44.         pos=19;
    
45.         shrt=res[1]>>2;
    
46.         res[0]=0;
    
47.     }
    
48.     else{
    
49.         while(res[1]>=pow10[pos+1])++pos;
    
50.         shrt=res[1]<<64-shrn[pos]|res[0]>>shrn[pos];
    
51.         res[0]=res[1]>>shrn[pos]?inv10[pos]:0;
    
52.     }
    
53.    
    
54.     _umul128(shrt,inv10[pos],res+1);
    
55.     return res[0]+res[1];
    
56. }
    
57. 
    
58. int main(){
    
59.     mpfr::mpreal::set_default_prec(192);
    
60.     std::cout.precision(39);
    
61. 
    
62.     uint64_t iinit=0xdeadbeefbadcab1e,isize=1000000000;
    
63.     uint64_t eval=einit(iinit);
    
64.     for(uint64_t i=iinit+1;i<=iinit+isize;++i){
    
65.         eval=emul(eval,i);
    
66.     }
    
67.     std::cout<<eval<<std::endl;
    
68.     std::cout<<einit(iinit+isize)<<std::endl;
    
69.     return 0;
    
70. }
    

复制代码

mathe

最关键是误差分析，乘百万次再调整精度应该不够了

liangbch

Mark, 以后有时间看

Ickiverar

2019.12.11 16:42，在我的老电脑上四线程并行开算。
16秒40亿，一天应该能算21万亿。
更新：
18:18算出前12个数，验证了20楼的结果：
d =  0, n = 9
d =  1, n = 62
d =  2, n = 62
d =  3, n = 10044
d =  4, n = 50583
d =  5, n = 1490717
d =  6, n = 5573998
d =  7, n = 65630447
d =  8, n = 688395641
d =  9, n = 5777940569
d = 10, n = 77773146302
d = 11, n = 1154318938997
d = 12, n = 1544607046599
done to n =        1552000000000, maxerr = 0.109788

liangbch

我尝试了把105#的代码修改为可在64位Linux下可用gcc编译通过的代码。
1. 在GCC下，不需要这句 “#pragma comment(lib,"mpfr-4.0.1.lib")”
2. 在GCC下，没有"intrin.h"这个文件，所以也就没有_umul128这个函数。
3. 在GCC下，10000000000000000000要改为 10000000000000000000ULL

修改思路，
首先可实现 c语言版的_umul128。修改后的部分代码如下。

1. #include <stdio.h>
   
2. #include <stdlib.h>
   
3. #include <stdint.h>
   
4. #include <iostream>
   
5. #include <mpreal.h>
   
6. 
   
7. typedef unsigned __int128  uint128_t;
   
8. 
   
9. uint64_t _umul128( uint64_t x,  uint64_t y,  uint64_t *HighPart)
   
10. {
    
11.     uint128_t __res = (uint128_t)x * y;
    
12.     *HighPart = (uint64_t) (__res >> 64);
    
13.      return (uint64_t) __res;
    
14. }
    

复制代码

代码通译通过后，可考虑改为汇编语言的实现。
假设上面的包含 _umul128的源代码编译为 test.o(目标文件，非可执行文件），执行下面的命令可看到这个函数汇编后的代码。

1. objdump -d test.o -M intel

复制代码

稍加修改，改为nasm格式的汇编文件，保存为umul128.asm

1. global _umul128
   
2. 
   
3. section .text
   
4. _umul128:
   
5.      mov    rax,rdi
   
6.      mov    rcx,rdx
   
7.      mul    rsi
   
8.      mov    QWORD [rcx],rdx
   
9.      ret
   

复制代码

3. 将_umul128的c语言的实现删除，改为声明

1. extern "C" uint64_t _umul128( uint64_t x,  uint64_t y,  uint64_t *HighPart);

复制代码

4. 现在，就可以使用汇编版的_umul128了，下面的编译链接命令

1. g++ -m64 -c test.cpp -O2 -o test.o
   
2. nasm -f elf64 -o umul128.o umul128.asm
   
3. g++ test.o umul128.o -lmpfr -o test
   

复制代码

liangbch

liangbch 发表于 2019-12-13 15:03
我尝试了把105#的代码修改为可在64位Linux下可用gcc编译通过的代码。
1. 在GCC下，不需要这句 “#pragma c ...

从性能上将，c语言的版本并不比汇编版本差。不过，将此函数替换为nasm格式的汇编代码，使得代码在跨平台(windows/Liunx)方面，更有优势，因为nasm本身就是跨平台的。