阶乘和圆周率

无心人

20# 为检索方便 记录下20#的结果的链接 有时间考虑下上面提出的每隔1024个数字，或者是1000个数字用striling公式计算精确结果 再进行累加 的方法 ========================================================== 还有个前提，数字必须转化成小于等于10的小数，否则会影响精度

无心人

#include #include #include #include #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 #include #include #include #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