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个数据
原帖由 无心人 于 2008-12-17 09:56 发表 http://bbs.emath.ac.cn/images/common/back.gif
...
顺便说下
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的检验结果
应该程序错误。128比特长度应该够用了，用Windows计算器都能够算出来结果不对:lol
我想了，我那种预先处理到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