刚反复的测试，觉得结果推进到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万亿的过程
如果我们写一个十进制定点数的乘法算法，只保留十进制前几十位数，使得从1乘到很大的N都不会产生严重误差，是不是可以避免浮点斯特令公式？
我写了一个算法，可以单线程在18秒内搜索10亿的范围（测试起点 0xdeadbeefbadcab1e ， 测试长度 10^9 ），可以任意并行化，最高可以搜索到 2^64 - 2 的阶乘。
但我的电脑正在全线程计算海绵球的体积，所以这个项目得过几天开算。
我先把我的代码贴出来吧。
使用到了mpfrc++，可以在 http://www.holoborodko.com/pavel/mpfr/ 下载。
用64位无符号定点数存储阶乘的十进制的开头，具体算法见 einit。
然后用emul逐项乘下一个数进行递推（可以每递推百万次用einit校正一次）。
搜索比较接近 pi/10*2^64 的值就行了。
#include<mpreal.h>
#pragma comment(lib,"mpfr-4.0.1.lib")
#include<iostream>
#include<intrin.h>
typedef mpfr::mpreal real;
//10^frac(log10(factorial(n)))*2^64/10
uint64_t einit(uint64_t n){
static const real log10e=1/mpfr::log(10);
static const real p264=mpfr::pow(2,64);
return (uint64_t)(p264*mpfr::pow(10,mpfr::frac(mpfr::lngamma(n+1)*log10e)-1));
}
//2^64/10 ~< e*n/10^k < 2^64
uint64_t emul(uint64_t e,uint64_t n){
//10^pos
static uint64_t pow10[]={
1,10,100,1000,10000,
100000,1000000,10000000,100000000,1000000000,
10000000000,100000000000,1000000000000,10000000000000,100000000000000,
1000000000000000,10000000000000000,100000000000000000,1000000000000000000,10000000000000000000
};
static int shrn[]={3,6,9,13,16,19,23,26,29,33,36,39,43,46,49,53,56,59,63,66};
//round(2^(64+shrn)/10^(pos+1))
static uint64_t inv10[]={
0xcccccccccccccccd,0xa3d70a3d70a3d70a,0x83126e978d4fdf3b,0xd1b71758e219652c,
0xa7c5ac471b478423,0x8637bd05af6c69b6,0xd6bf94d5e57a42bc,0xabcc77118461cefd,
0x89705f4136b4a597,0xdbe6fecebdedd5bf,0xafebff0bcb24aaff,0x8cbccc096f5088cc,
0xe12e13424bb40e13,0xb424dc35095cd80f,0x901d7cf73ab0acd9,0xe69594bec44de15b,
0xb877aa3236a4b449,0x9392ee8e921d5d07,0xec1e4a7db69561a5,0xbce5086492111aeb
};
//res = e*n
uint64_t res;
res=_umul128(e,n,res+1);
if(res==0)return res;
//shrt = res >> shrn
int pos=0;
uint64_t shrt;
if(res>=pow10){
pos=19;
shrt=res>>2;
res=0;
}
else{
while(res>=pow10)++pos;
shrt=res<<64-shrn|res>>shrn;
res=res>>shrn?inv10:0;
}
_umul128(shrt,inv10,res+1);
return res+res;
}
int main(){
mpfr::mpreal::set_default_prec(192);
std::cout.precision(39);
uint64_t iinit=0xdeadbeefbadcab1e,isize=1000000000;
uint64_t eval=einit(iinit);
for(uint64_t i=iinit+1;i<=iinit+isize;++i){
eval=emul(eval,i);
}
std::cout<<eval<<std::endl;
std::cout<<einit(iinit+isize)<<std::endl;
return 0;
}
最关键是误差分析，乘百万次再调整精度应该不够了
Mark, 以后有时间看
本帖最后由 Ickiverar 于 2019-12-11 18:18 编辑
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
我尝试了把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。修改后的部分代码如下。
#include <stdio.h>
#include <stdlib.h>
#include <stdint.h>
#include <iostream>
#include <mpreal.h>
typedef unsigned __int128uint128_t;
uint64_t _umul128( uint64_t x,uint64_t y,uint64_t *HighPart)
{
uint128_t __res = (uint128_t)x * y;
*HighPart = (uint64_t) (__res >> 64);
return (uint64_t) __res;
}
代码通译通过后，可考虑改为汇编语言的实现。
假设上面的包含 _umul128的源代码编译为 test.o(目标文件，非可执行文件），执行下面的命令可看到这个函数汇编后的代码。
objdump -d test.o -M intel
稍加修改，改为nasm格式的汇编文件，保存为umul128.asm
global _umul128
section .text
_umul128:
mov rax,rdi
mov rcx,rdx
mul rsi
mov QWORD ,rdx
ret
3. 将_umul128的c语言的实现删除，改为声明
extern "C" uint64_t _umul128( uint64_t x,uint64_t y,uint64_t *HighPart);
4. 现在，就可以使用汇编版的_umul128了，下面的编译链接命令
g++ -m64 -c test.cpp -O2 -o test.o
nasm -f elf64 -o umul128.o umul128.asm
g++ test.o umul128.o -lmpfr -o test
liangbch 发表于 2019-12-13 15:03
我尝试了把105#的代码修改为可在64位Linux下可用gcc编译通过的代码。
1. 在GCC下，不需要这句 “#pragma c ...
从性能上将，c语言的版本并不比汇编版本差。不过，将此函数替换为nasm格式的汇编代码，使得代码在跨平台(windows/Liunx)方面，更有优势，因为nasm本身就是跨平台的。