阶乘和圆周率

liangbch

回复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位浮点了
唯一的好处是这个程序基本是线性时间

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