无心人 发表于 2008-12-20 22:15:42

可以发现上面几个例子的结果是两个不同的输出
测试发现正确的输出是
149832529!!!!
但计算器计算发现,似乎此时的对数值精度并不高

无心人 发表于 2008-12-20 22:23:15

重新修改成
#include <gmp.h>
#include <mpfr.h>
#include <stdio.h>
#include <math.h>
#include <sys/time.h>

#define DefaultBitLength 128

mpfr_t ln10, cA, cT;
//=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_init2(s, DefaultBitLength);

    mpfr_div(s, cT, a, GMP_RNDD); //s = 1/12a = 1/12n   
    mpfr_sub(s, s, a, GMP_RNDD); //s = s - 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, cA, GMP_RNDD); //最终结果s=1/2ln(2*pi) + (n + 1/2)ln(n) - n + 1/12n

    mpfr_mul(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, used_time, LT2;
int find = 0, test1 = 0;
long i = 0, j = 0;
struct timeval start_time, end_time;

mpfr_init2(ln10, DefaultBitLength);
mpfr_init2(cA, DefaultBitLength);
mpfr_init2(cT, DefaultBitLength);
mpfr_set_ui(ln10, 10, GMP_RNDD);
mpfr_set_ui(cT, 1, GMP_RNDD);

mpfr_log(ln10, ln10, GMP_RNDD);
mpfr_div(ln10, cT, ln10, GMP_RNDD);

mpfr_const_pi(cA, GMP_RNDD);
mpfr_mul_ui(cA, cA, 2, GMP_RNDD);
mpfr_log(cA, cA, GMP_RNDD);
mpfr_div_ui(cA, cA, 2, GMP_RNDD);

mpfr_div_ui(cT, cT, 12, GMP_RNDD);

printf("请输入起始搜索数字:");
scanf("%lf", &start);
printf("请输入上界: ");
scanf("%lf", &hi);
printf("请输入下界: ");
scanf("%lf", &low);

LT2 = log10(2);
hi = log10(hi);
low = log10(low);

printf("起始位置:%.0f\n当前界限[%.17f, %.17f]\n", start, low, hi);
a = start;
b = 1.0;
while (a >= 10.0)
{
    a /= 10.0;
    b *= 10.0; ///////////////////////////////////
}

lf = 2.0;
gettimeofday(&start_time, NULL);
for (i = 0; i < 1000; i ++)
   lf += stirling(lf + 100000000.0);
gettimeofday(&end_time, NULL);
used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec)) / 1000.0;
printf("高精度stirling函数执行时间(毫秒): %8.6lf, 测试值:%.8lf\n", used_time, lf);

lf = stirling(start);
printf("起始对数: %.17f\n", lf);

gettimeofday(&start_time, NULL);

while (1)
{
    start += 1.0;
    a = log2(start/b) * LT2;///////////////////////////
   
    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 >= 1000)
    {
//      lf = stirling(start);
      i = 0;
   j ++;
   if (j >= 1000)
   {
if (test1 == 0)
{
          gettimeofday(&end_time, NULL);
          used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec))/1000000.0;
          printf("每10^6次查找执行时间(秒): %.6lf\n", used_time);
   test1 = 1;
      }
   }
    }
}

if (find)
{
    printf("\n当前界限下的阶乘%.0f!, 对应的对数%.15f", start, lf);
}

mpfr_clear(ln10);
mpfr_clear(cA);
mpfr_clear(cT);
return 1;
}

编译命令
gccfactPi1.c -o factpi1.exe -lmpfr -lgmp -lm
======================================================
请输入起始搜索数字:100000000
请输入上界: 3.1415927
请输入下界: 3.1415926
起始位置:100000000
当前界限
高精度stirling函数执行时间(毫秒): 30.043000, 测试值:493.74586181
起始对数: 0.20876475177585671
每10^6次查找执行时间(秒): 0.300432

当前界限下的阶乘149832529!, 对应的对数0.497149873638380
============================================
计算器得到的值
0.4971498735933375434005
============================================
直接使用10的对数
其他条件不变的输出
即使用    a = log10(start/b);// * LT2;取代a = log2(start / b) * LT2;
============================================
请输入起始搜索数字:100000000
请输入上界: 3.1415927
请输入下界: 3.1415926
起始位置:100000000
当前界限
高精度stirling函数执行时间(毫秒): 30.044000, 测试值:493.74586181
起始对数: 0.20876475177585671
每10^6次查找执行时间(秒): 0.450648

当前界限下的阶乘149832529!, 对应的对数0.497149873595150
===========================================
精度大些,但消耗时间多了
即使使用log10(2)的高精度值直接乘而不是用对数函数
也不能增大精度,所以似乎库函数存在某些优化????

无心人 发表于 2008-12-20 22:44:22

程序可以确定为62#程序但
使用a = log10(start/b);// * LT2;取代a = log2(start / b) * LT2;
此时
不同的优化结果如下
=========================================
不使用任何优化
每10^6次查找执行时间(秒): 0.470677
当前界限下的阶乘149832529!, 对应的对数0.4971 4987 3595 150
=========================================
带-O2
每10^6次查找执行时间(秒): 0.450648
当前界限下的阶乘149832529!, 对应的对数0.4971 4987 3593 083
===========================================
带-O2 -ffast-math
每10^6次查找执行时间(秒): 0.160230
当前界限下的阶乘149832529!, 对应的对数0.4971 4987 3593 496
===========================================
对比计算器得到的数值0.4971 4987 3593 3375 4340 05
==========================================
反而是最后一种情况又快, 精度又高
呵呵

无心人 发表于 2008-12-21 21:37:29

修改了代码
使用gcc -O2 -ffast-math -march=pentium4 factpi1.c -o factpi1.exe -lmpfr -lgmp
#include <gmp.h>
#include <mpfr.h>
#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <sys/time.h>

#define DefaultBitLength 128
#define TEST

mpfr_t ln10, cA, cT;

FILE * fp;

//=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_init2(s, DefaultBitLength);

    mpfr_div(s, cT, a, GMP_RNDD); //s = 1/12a = 1/12n   
    mpfr_sub(s, s, a, GMP_RNDD); //s = s - 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, cA, GMP_RNDD); //最终结果s=1/2ln(2*pi) + (n + 1/2)ln(n) - n + 1/12n

    mpfr_mul(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, used_time;
int find = 0, test1 = 0;
long i, j = 0;
struct timeval start_time, end_time;

mpfr_init2(ln10, DefaultBitLength);
mpfr_init2(cA, DefaultBitLength);
mpfr_init2(cT, DefaultBitLength);
mpfr_set_ui(ln10, 10, GMP_RNDD);
mpfr_set_ui(cT, 1, GMP_RNDD);

mpfr_log(ln10, ln10, GMP_RNDD);
mpfr_div(ln10, cT, ln10, GMP_RNDD);

mpfr_const_pi(cA, GMP_RNDD);
mpfr_mul_ui(cA, cA, 2, GMP_RNDD);
mpfr_log(cA, cA, GMP_RNDD);
mpfr_div_ui(cA, cA, 2, GMP_RNDD);

mpfr_div_ui(cT, cT, 12, GMP_RNDD);

printf("请输入起始搜索数字:");
scanf("%lf", &start);
printf("请输入上界: ");
scanf("%lf", &hi);
printf("请输入下界: ");
scanf("%lf", &low);
fp = fopen("factpi.log", "a+t");
printf("起始位置:%.0f\n当前界限[%.17f, %.17f]\n", start, low, hi);
fprintf(fp, "\n///////////////////////////////////////////////////////////////\n");
fprintf(fp, "起始位置:%.0f\n当前界限[%.17f, %.17f]\n", start, low, hi);

hi = log10(hi);
low = log10(low);

printf("当前对数界限[%.17f, %.17f]\n", low, hi);
fprintf(fp, "当前对数界限[%.17f, %.17f]\n", low, hi);
a = start;
b = 1.0;
while (a >= 10.0)
{
    a /= 10.0;
    b *= 10.0;
}

lf = 2.0;
gettimeofday(&start_time, NULL);
for (i = 0; i < 1000; i ++)
   lf += stirling(lf + 100000000.0);
gettimeofday(&end_time, NULL);
used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec)) / 1000.0;
printf("高精度stirling函数执行时间(毫秒): %8.6lf, 测试值:%.8lf\n", used_time, lf);
fprintf(fp, "高精度stirling函数执行时间(毫秒): %8.6lf, 测试值:%.8lf\n", used_time, lf);

lf = stirling(start);
printf("起始对数: %.17f\n", lf);
fprintf(fp, "起始对数: %.17f\n", lf);
fflush(fp);
gettimeofday(&start_time, NULL);

i = 0;
while (1)
{   
    if ((a = log10((start += 1.0)/b)) >= 1.0)
    {
      b *= 10.0;
      a -= 1.0;
    }
      
    if ((lf += a) >= 1.0)
      lf -= 1.0;

//    printf("[%.0f, %.15f]", start, lf);
   
    if ((lf >= low) && (lf < hi))
   {
       find = 1;
       break;
   }
   
    if (( ++i) >= 1000)
    {
      lf = stirling(start);
      i = 0;
   j ++;
#ifdef TEST
   if (j >= 1000)
   {
if (test1 == 0)
{
          gettimeofday(&end_time, NULL);
          used_time = ((end_time.tv_sec - start_time.tv_sec)*1000000.0 + (end_time.tv_usec - start_time.tv_usec))/1000000.0;
          printf("每10^6次查找执行时间(秒): %.6lf\n", used_time);
   test1 = 1;
      }
   }
#else
   if (j >= 100000000)
   {
j = 0;
fprintf(fp, "当前搜索数值: %.0f\n", start);
fflush(fp);
   }
#endif
    }
}

if (find)
{
    printf("\n当前界限下的阶乘%.0f!, 对应的对数%.15f", start, lf);
    fprintf(fp, "\n当前界限下的阶乘%.0f!, 对应的对数%.15f", start, lf);
}

mpfr_clear(ln10);
mpfr_clear(cA);
mpfr_clear(cT);
fclose(fp);
return 1;
}

==========================================================
请输入起始搜索数字:100000000
请输入上界: 3.1415927
请输入下界: 3.1415926
起始位置:100000000
当前界限
当前对数界限
高精度stirling函数执行时间(毫秒): 30.044000, 测试值:493.74586181
起始对数: 0.20876475177585671
每10^6次查找执行时间(秒): 0.170244

当前界限下的阶乘149832529!, 对应的对数0.497149873593338
========================================================
经过查询汇编代码,程序已经被gcc优化到了很好的程度
相对以前,执行时间已经缩小了很多了
使用了P6以上的fcomi和直接的log指令, 另外,其最终的对数值精度提高了
==========================================================
P4 2.0 /1G DDR/XP SP3

无心人 发表于 2008-12-22 07:43:57

d = 314159265358,n = 1154318938997!, 对应的对数0.497149872693353,首数字3141592653584138
=============================================================
通过PARI验证
(09:52) gp > n = 1154318938997
%38 = 1154318938997
(07:44) gp > (n+1/2)*log(n) + 1/2*log(2*pi) - n + 1/(12*n) - 1/(360*n^3)
%39 = 30906348934959.82738115133341720201016520726795684547468467
(07:44) gp > %39 / log(10)
%40 = 13422456798229.49714987269335219103203453332675692373542468
(07:44) gp > %40 - floor(%40)
%41 = 0.497149872693352191032034533326756923735424682255
(07:44) gp > 10^(%41)
%42 = 3.14159265358413885452821840174232968895743730673

无心人 发表于 2008-12-22 08:39:24

用来计算的机器是
Core 2 XEON 1.6G / 1G DDR2
下一个将在10万亿附近产生
目前增加到10000个数字校正一次
每10 ^8的时间是7.8125秒
如果还是每1000校正一次
则是9.375秒
==========================
差距是两天
呵呵
==========================
依据是对最新求得的结果1154318938997
使用10000间隔校正
从之前的11543000000重新搜索
结果对数是0.497149872693352
而GP得到的结果是0.497149872693353
===========================
目前的对数范围

差距是13822
但因为所有代码均被gcc优化成汇编代码,而gcc默认是63位浮点
所以实际上的差距是1382200
我想足够10000间隔校正而不丢失分辨精度吧
=============================
呵呵,求出下个结果,就不再计算了
再下个至少需要100天

无心人 发表于 2008-12-22 11:00:46

在这十天的空隙里
先贴点以10^n开始的
d = 1, n = 5
d = 10, n = 27
d = 100, n = 197
d = 1000, n = 197
d = 10000, n = 7399
////////////////////////////////////////////
呵呵,下面的似乎是错误的
看来需要校正结果

无心人 发表于 2008-12-23 08:00:46

竟然昨天下午5:37就得到了数据
=============================
请输入起始搜索数字:1154318938997
请输入上界: 3.141592653590
请输入下界: 3.141592653589
起始位置:1154318938997
当前界限
当前对数界限
高精度stirling函数执行时间(毫秒): 15.625000, 测试值:493.74586181
起始对数: 0.49714987269335220
每10^6次查找执行时间(秒): 0.078125

当前界限下的阶乘1544607046599!, 对应的对数0.497149872694156
=================================================
GP的验证
(11:04) gp > stirling1(n) = (n+1/2)*log(n) + 1/2*log(2*pi) - n + 1/(12*n) - 1/(
360*n^3)
(11:04) gp > stirling(n) = stirling1(n) / log(10) - floor(stirling1(n) / log(10
))
(11:06) gp > stirling(1544607046599)
%57 = 0.497149872694155502964270816144690329633515005612
(07:54) gp > 10^%57
%58 = 3.14159265358994983986193148792491264569132697520

gxqcn 发表于 2008-12-23 09:23:25

看兄弟比较孤寂的,回复一个。

可惜没有这么好的软硬件条件,所以只能欣赏而已。

无心人 发表于 2008-12-23 09:30:31

呵呵
这个情况
导致了我还有继续计算下去的必要
再算一个,直到结果超越10万亿
页: 1 2 3 4 5 6 [7] 8 9 10 11 12 13 14
查看完整版本: 阶乘和圆周率