阶乘和圆周率

无心人

20#
为检索方便
记录下20#的结果的链接
有时间考虑下上面提出的每隔1024个数字，或者是1000个数字用striling公式计算精确结果
再进行累加
的方法

==========================================================
还有个前提，数字必须转化成小于等于10的小数，否则会影响精度

无心人

#include <gmp.h>
#include <mpfr.h>
#include <stdio.h>
#include <math.h>

#define DefaultBitLength 128

  mpfr_t ln10, A;
//=1/2ln(2π)+(n+1/2)ln(n)-n+1/12n
//n > 10^8 误差<1/360Z^3 = 10^(-26)
long double stirling(double n)
{
    mpfr_t a, b, c, s;
    long double r;
   
    mpfr_init2(a, DefaultBitLength);
    mpfr_init2(b, DefaultBitLength);
    mpfr_set_d(a, n, GMP_RNDD); //a=n
    mpfr_init2(c, DefaultBitLength);
    mpfr_set_d(b, 12.0, GMP_RNDD); //b=12
    mpfr_set_d(c, 1.0, GMP_RNDD); //c=1
    mpfr_init2(s, DefaultBitLength);

    mpfr_div(b, c, b, GMP_RNDD); // b = c / b = 1/12
   
    mpfr_div(b, b, a, GMP_RNDD); // b = b / a = 1/12n
   
    mpfr_sub(s, b, a, GMP_RNDD); //s = b - a = - n + 1/12n
   
    mpfr_log(c, a, GMP_RNDD); //c = ln(n)

    mpfr_set_d(b, 0.5, GMP_RNDD);
    mpfr_add(b, a, b, GMP_RNDD);
    mpfr_mul(b, b, c, GMP_RNDD); //b = (n + 1/2)ln(n)
   
    mpfr_add(s, s, b, GMP_RNDD); //s = (n + 1/2)ln(n) - n + 1/12n

    mpfr_add(s, s, A, GMP_RNDD); //最终结果s=1/2ln(2*pi) + (n + 1/2)ln(n) - n + 1/12n

    mpfr_div(s, s, ln10, GMP_RNDD); //转换为log10
    mpfr_floor(b, s);
    mpfr_sub(s, s, b, GMP_RNDD); //只要小数
   
    r = mpfr_get_ld(s, GMP_RNDD);
   
    mpfr_clear(a);
    mpfr_clear(b);
    mpfr_clear(c);
    mpfr_clear(s);
    return r;
}

int main(void)
{
  double a, b, low, hi, lf, t;
  double start;
  int find = 0;
  long i = 0, j = 0;
  
  mpfr_init2(ln10, DefaultBitLength);
  mpfr_init2(A, DefaultBitLength);
  mpfr_set_d(ln10, 10.0, GMP_RNDD);
  mpfr_log(ln10, ln10, GMP_RNDD);

  mpfr_const_pi(A, GMP_RNDD);
  mpfr_mul_ui(A, A, 2, GMP_RNDD);
  mpfr_log(A, A, GMP_RNDD);
  mpfr_div_ui(A, A, 2, GMP_RNDD);
  
  printf("请输入起始搜索数字：");
  scanf("%lf", &start);
  printf("请输入上界: ");
  scanf("%lf", &hi);
  printf("请输入下界: ");
  scanf("%lf", &low);
  
  hi = log10(hi);
  low = log10(low);
  
  printf("起始位置：%.0f\n当前界限[%.15f, %.15f]\n", start, low, hi);
  a = start;
  b = 1.0;
  while (a >= 10.0)
  {
    a /= 10.0;
    b /= 10.0;
  }

  lf = stirling(start);
  
  printf("a: %.15f\nb: %.15f\nlf: %.15f\n", a, b, lf);
  while (1)
  {
    a += b;
    start += 1.0;
    if (a >= 10.0)
    {
        a /= 10.0;
        b /= 10.0;
    }
   
    i ++;
    lf += log10(a);
        
    if (lf >= 1.0)
        lf -= 1.0;

//    printf("[%.0f, %.15f]  ", start, lf);
   
    if ((lf >= low) && (lf < hi))
     {
       find = 1;
       break;
     }
   
    if (i >= 100)
    {
      lf = stirling(start);
      i = 0;
      j ++;
      if (j >= 100000)
      {
          j = 0;
          printf(".");
      }
    }
  }
  
  if (find)
  {
    printf("\n当前界限下的阶乘%.0f!, 对应的对数%.15f", start, lf);
  }
  
  mpfr_clear(ln10);
  mpfr_clear(A);
  return 1;
}

===========================
终于搞定
顺便说下
16#公式是错误的，害死我了
因为偷懒没去找，照抄的
应该是(z+1/2)
===========================
大概每4秒搜索10^8个数据

mathe

原帖由 无心人 于 2008-12-17 09:56 发表
...
顺便说下
16#公式是错误的，害死我了
因为偷懒没去找，照抄的
应该是(z+1/2)

...

公式应该没有错，抄自wolframe的网站上，可能是$\Gamma(z)$和$\Gamma(z+1)$的区别，注意$n! =\Gamma(n+1)$

无心人

呵呵
知道了

得到了最新的结果
当前界限下的阶乘13655917103!, 对应的对数0.497149872668885
尾数3141592653407
=============================================
经检验，似乎不对？？
精度达不到么？
还是程序错误？

无心人

(11:38) gp > n = 13655917103
%14 = 13655917103
(11:39) gp > (n+1/2)*log(n) + 1/2*log(2*pi) - n + 1/(12*n) - 1/(360*n^3)
%15 = 305038211903.8966681480956355012460677100090297729522643543
(11:39) gp > %15 / log(10)
%16 = 132476412199.4971477766334402275974580068172942982614808834

这是PARI的检验结果

mathe

应该程序错误。128比特长度应该够用了，用Windows计算器都能够算出来结果不对

无心人

我想了，我那种预先处理到10以下的方法不对
那么做，增量必然是个很小的值
恐怕是在那里出问题了

应该每次都用除以10^n的方式，以保持精度

无心人

#include <gmp.h>
#include <mpfr.h>
#include <stdio.h>
#include <math.h>

#define DefaultBitLength 128

  mpfr_t ln10, A;
//=1/2ln(2π)+(n+1/2)ln(n)-n+1/12n
//n > 10^8 误差<1/360Z^3 = 10^(-26)
long double stirling(double n)
{
    mpfr_t a, b, c, s;
    long double r;
   
    mpfr_init2(a, DefaultBitLength);
    mpfr_init2(b, DefaultBitLength);
    mpfr_set_d(a, n, GMP_RNDD); //a=n
    mpfr_init2(c, DefaultBitLength);
    mpfr_set_d(b, 12.0, GMP_RNDD); //b=12
    mpfr_set_d(c, 1.0, GMP_RNDD); //c=1
    mpfr_init2(s, DefaultBitLength);

    mpfr_div(b, c, b, GMP_RNDD); // b = c / b = 1/12
   
    mpfr_div(b, b, a, GMP_RNDD); // b = b / a = 1/12n
   
    mpfr_sub(s, b, a, GMP_RNDD); //s = b - a = - n + 1/12n
   
    mpfr_log(c, a, GMP_RNDD); //c = ln(n)

    mpfr_set_d(b, 0.5, GMP_RNDD);
    mpfr_add(b, a, b, GMP_RNDD);
    mpfr_mul(b, b, c, GMP_RNDD); //b = (n + 1/2)ln(n)
   
    mpfr_add(s, s, b, GMP_RNDD); //s = (n + 1/2)ln(n) - n + 1/12n

    mpfr_add(s, s, A, GMP_RNDD); //最终结果s=1/2ln(2*pi) + (n + 1/2)ln(n) - n + 1/12n

    mpfr_div(s, s, ln10, GMP_RNDD); //转换为log10
    mpfr_floor(b, s);
    mpfr_sub(s, s, b, GMP_RNDD); //只要小数
   
    r = mpfr_get_ld(s, GMP_RNDD);
   
    mpfr_clear(a);
    mpfr_clear(b);
    mpfr_clear(c);
    mpfr_clear(s);
    return r;
}

int main(void)
{
  double a, b, low, hi, lf, t;
  double start;
  int find = 0;
  long i = 0, j = 0;
  
  mpfr_init2(ln10, DefaultBitLength);
  mpfr_init2(A, DefaultBitLength);
  mpfr_set_d(ln10, 10.0, GMP_RNDD);
  mpfr_log(ln10, ln10, GMP_RNDD);

  mpfr_const_pi(A, GMP_RNDD);
  mpfr_mul_ui(A, A, 2, GMP_RNDD);
  mpfr_log(A, A, GMP_RNDD);
  mpfr_div_ui(A, A, 2, GMP_RNDD);
  
  printf("请输入起始搜索数字：");
  scanf("%lf", &start);
  printf("请输入上界: ");
  scanf("%lf", &hi);
  printf("请输入下界: ");
  scanf("%lf", &low);
  
  hi = log10(hi);
  low = log10(low);
  
  printf("起始位置：%.0f\n当前界限[%.15f, %.15f]\n", start, low, hi);
  a = start;
  b = 1.0;
  while (a >= 10.0)
  {
    a /= 10.0;
    b *= 10.0;
  }

  lf = stirling(start);
  
  printf("a: %.15f\nb: %.0f\nlf: %.15f\n", a, b, lf);
  while (1)
  {
    start += 1.0;
    a = log10(start/b);
   
    if (a >= 1.0)
    {
      b *= 10.0;
      a -= 1.0;
    }
   
    i ++;
    lf += a;
        
    if (lf >= 1.0)
        lf -= 1.0;

//    printf("[%.0f, %.15f]  ", start, lf);
   
    if ((lf >= low) && (lf < hi))
     {
       find = 1;
       break;
     }
   
    if (i >= 100)
    {
      lf = stirling(start);
      i = 0;
      j ++;
      if (j >= 100000)
      {
          j = 0;
          printf(".");
      }
    }
  }
  
  if (find)
  {
    printf("\n当前界限下的阶乘%.0f!, 对应的对数%.15f", start, lf);
  }
  
  mpfr_clear(ln10);
  mpfr_clear(A);
  return 1;
}
===========================================
修改了程序
重新来过

无心人

另外一个问题，似乎目前的库，对long double都不支持
我修改程序用long double，取不到对数值

而目前的double精度恐怕支持不了多久了

无心人

终于正确了
当前界限下的阶乘5777940569!, 对应的对数0.497149872614068
起始数字314159265301