回复29楼，那是很久以前的事了，代码不好找了。
#include <gmp.h>
#include <mpfr.h>
#include <stdio.h>
#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位浮点了
唯一的好处是这个程序基本是线性时间
没有问题，我用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编译
在把浮点数保存在文件里时遇到困难
不知道如何写
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