阶乘和圆周率

无心人 · 发表于 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;
}

复制代码

编译命令
gcc factPi1.c -o factpi1.exe -lmpfr -lgmp -lm
======================================================
请输入起始搜索数字：100000000
请输入上界: 3.1415927
请输入下界: 3.1415926
起始位置：100000000
当前界限[0.49714986528586863, 0.49714987910989134]
高精度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
当前界限[0.49714986528586863, 0.49714987910989134]
高精度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
当前界限[3.14159260000000010, 3.14159269999999990]
当前对数界限[0.49714986528586863, 0.49714987910989134]
高精度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
===========================
目前的对数范围
[0.49714987269402422 , 0.49714987269416244]
差距是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
当前界限[3.14159265358900000, 3.14159265359000010]
当前对数界限[0.49714987269402422, 0.49714987269416244]
高精度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万亿

账号		自动登录	找回密码
密码			欢迎注册

[讨论] 阶乘和圆周率