阶乘和圆周率

liangbch

回复29楼，那是很久以前的事了，代码不好找了。

无心人

#include #include #include #define DefaultBitLength 128 int main(void) { mpfr_t a, b, low, hi, lf, t; double dLow, dHi, dLf, dA; int find = 0; long i = 0; mpfr_init2(a, DefaultBitLength); mpfr_set_d(a, 2.0, GMP_RNDD); mpfr_init2(b, DefaultBitLength); mpfr_set_d(b, 1.0, GMP_RNDD); mpfr_init2(low, DefaultBitLength); mpfr_init2(hi, DefaultBitLength); mpfr_init2(lf, DefaultBitLength); mpfr_set_d(lf, 0.0, GMP_RNDD); mpfr_init2(t, DefaultBitLength); printf("请输入上界: "); scanf("%lf", &dHi); printf("请输入下界: "); scanf("%lf", &dLow); // printf("当前界限[%f, %f]\n", dLow, dHi); mpfr_set_d(low, dLow, GMP_RNDD); mpfr_set_d(hi, dHi, GMP_RNDD); mpfr_log10(low, low, GMP_RNDD); mpfr_log10(hi, hi, GMP_RNDD); while (1) { mpfr_log10(t, a, GMP_RNDD); mpfr_add(lf, lf, t, GMP_RNDD); if (mpfr_cmp_d(lf, 1.0) >= 0) mpfr_sub_ui(lf, lf, 1, GMP_RNDD); if ((mpfr_cmp(lf, low) >= 0) && (mpfr_cmp(lf, hi) < 0)) { find = 1; break; } mpfr_add(a, a, b, GMP_RNDD); if (mpfr_cmp_d(a, 10.0) >= 0) { i++; mpfr_div_ui(b, b, 10, GMP_RNDD); mpfr_add_ui(a, b, 1, GMP_RNDD); } } if (find) { dLf = mpfr_get_d1(lf); while (i > 0) { mpfr_mul_ui(a, a, 10, GMP_RNDD); i--; } dA = mpfr_get_d1(a); printf("当前界限下的阶乘%.0f!, 对应的对数%.15f", dA, dLf); } mpfr_clear(a); mpfr_clear(b); mpfr_clear(low); mpfr_clear(hi); mpfr_clear(lf); mpfr_clear(t); return 1; } //////////////////////////////////////////////////////// 利用mpfr可写出适合超过100000000的阶乘的一个寻找程序 但似乎速度慢很多哦

无心人

使用这个库之后 速度巨慢了 上次的程序，十几分钟就能计算出 314159265的结果

无心人

经过一夜的运算找到了 起始数字314159265的结果 在20# 用PARI计算的0.497149872468489（程序输出的，转换为Double, 真实的是96Bit浮点）对应的数字 希望有人进行复检 另外，下个数字运算程序可能就要写定时状态保存了，而且要用128位浮点了 唯一的好处是这个程序基本是线性时间

mathe

没有问题，我用striling公式带入windows计算器验算了一下，结果为: 3.14159265196... 而Windows计算机浮点精度挺高的。

无心人

32#程序是在MinGW+ GCC4.21下编译通过的 预先编译了GMP 4.2.4和MPFR2.3.2，均是静态库形式 gcc factPi.c -std=c99 -O2 -lmpfr -lgmp -o factPi.exe 编译 可脱离MinGW运行 有兴趣的童鞋 我可以上传一个包 覆盖在MINGW的 include和lib就能使用

无心人

.h放mingw\include 其他放 mingw\lib 安装msys进入模拟unix环境gcc编译 GMP-MPFR-MinGW-static.rar (433.45 KB, 下载次数: 0)

2008-12-9 08:56 上传

点击文件名下载附件

无心人

在把浮点数保存在文件里时遇到困难 不知道如何写 MPFR这方面的资料太少了 mpfr_out_str不能用 mpfr_get_str怎么不能保存小数点？

无心人

今天想到一个新思路 计算了 log(3.141592653590) = 0.4971498726941624371131490186099 log(3.141592653589) = 0.49714987268172081667353310595166 二者的差是1.3824022710488533020286478680411e-13 约等于2^(-42.717887740820254037211382936331 ) 考虑double的精度是2^(-53) 显然只要最终数字的精度大于2^(-43)即可 那么如果我们用MPFR每隔1024隔数字 用strling公式 计算一个精度在2^(-53)的准确的阶乘的log值 而之间的数字则使用double类型的log进行连加 则中间状态的精度肯定大于2^(-43) =========================== 为了防止边界溢出现象 我们可以考虑使用long double类型 则此时的机器精度是2^(-63)

无心人

在long double情况下 有 Prelude> let a = 3.14159265358979 Prelude> let b = 3.14159265358980 Prelude> log(b) / log(10) - log(a) / log(10) 1.4432899320127035e-15 Prelude> log(it) / log(2) -49.29956028185891 Prelude> it + 10 -59.29956028185891 即保守估计可以计算到314159265358979