阶乘和圆周率

无心人 · 发表于 2008-12-17 15:14:22

发现50#用的程序是每10^7计数一次
所以以前说的每4秒 10^8错误
是40秒

无心人 · 发表于 2008-12-18 08:00:01

当前界限下的阶乘77773146302!, 对应的对数0.497149872686531
开始数字3141592653534
检验：
(15:09) gp > n = 77773146302
%26 = 77773146302
(07:55) gp > (n+1/2)*log(n) + 1/2*log(2*pi) - n + 1/(12*n) - 1/(360*n^3)
%27 = 1872548868987.616634238479898278998546965437108456887088346
(07:55) gp > %27 / log(10)
%28 = 813237640895.4971498726865435859265735865657796251777623316
(07:55) gp > %28 - floor(%28)
%29 = 0.497149872686543585926573586565779625177762331512
(07:56) gp > 10^%29
%30 = 3.14159265353488687304651737046701462136950721154

gxqcn · 发表于 2008-12-18 08:33:02

原帖由 无心人 于 2008-12-17 14:08 发表
另外一个问题，似乎目前的库，对long double都不支持
我修改程序用long double，取不到对数值

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

真正的 long double 有 80 bits，至少需要 10 字节存储，
所以一般编译器不支持。

但可以通过内嵌汇编得到80 bits精度的浮点数，
可惜FPU设计成栈模式，手工编写调试汇编代码非常累人（我之前曾写过一些）！

以下是 FPU 中一些可能对本主题讨论话题有帮助的指令：

FLDL2E Load the log base 2 of e (Napierian constant)

FLDL2T Load the log base 2 of Ten

FLDLG2 Load the log base 10 of 2 (common log of 2)

FLDLN2 Load the log base e of 2 (natural log of 2)

FLDPI    Load the value of PI

FSQRT    Square root of ST(0)

FYL2X    Y*Log2X

FYL2XP1 Y*Log2(X+1)

无心人 · 发表于 2008-12-18 09:06:09

写汇编代码在 gcc里并不方便
在windows下编译GMP和MPFR也不方便

呵呵
两难

无心人 · 发表于 2008-12-18 09:17:27

下一个选择就是5000亿左右了，估计是极限了
仅仅在两个界限的对数值最后四位有差别
[2780, 4162]
每次加1000个，再调整一次

估计差别值在范围上很小
有错过的可能
=============================
刚才加大了输出的精度
发现是6位的差别
我再去查下double的具体数值

无心人 · 发表于 2008-12-20 21:10:13

今天有时间测试了下程序的时间

#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;
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);
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 = 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 >= 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;
}

复制代码

===============================================
E:\math>factpi1
请输入起始搜索数字：100000000
请输入上界: 3.1415927
请输入下界: 3.1415926
起始位置：100000000
当前界限[0.49714986528586863, 0.49714987910989134]
高精度stirling函数执行时间(毫秒): 40.058000, 测试值：493.74586181
起始对数: 0.20876475177585671
每10^6次查找执行时间(秒): 0.460662

当前界限下的阶乘149832529!, 对应的对数0.497149873593380
=================================================
对应机器P4 2.0/1G DDR400
软件XP
编译器MinGW + GCC 4.3.2
使用GMP 4.2.4 + MPFR 2.3.2，编译为静态包
===================================================
可以看到MPFR下，每1000次stirling函数计算仅0.04秒
而double下每1000*1000次叠加对数值是需要0.46秒
可计算出高精度浮点的运算时间仅是double叠加对数的100倍

所以在目前情况下，普通函数计算称为瓶颈
另外，根据执行时间判断
每次计算需要800-1000时钟周期

无心人 · 发表于 2008-12-20 21:50:52

使用-ffast-math选项编译后的输出

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

当前界限下的阶乘143984543!, 对应的对数0.497149871271711

无心人 · 发表于 2008-12-20 21:55:52

57#同样条件下
但修改程序如下

#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;
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);
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 = 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 >= 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;
}

复制代码

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

当前界限下的阶乘149832529!, 对应的对数0.497149873593496

无心人 · 发表于 2008-12-20 21:57:48

可以看到
乘法和除法在浮点运算下
对时间的影响是很大的

还有一个可以修改优化的地方是对取10的对数
改为取2的对数再乘以log10(2)
但不知道标准函数有取2的对数的函数么？
还是需要自己写汇编

无心人 · 发表于 2008-12-20 22:13:45

去掉-ffast-math，同时增加了对求10的对数的优化

#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;
}

复制代码

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

当前界限下的阶乘149832529!, 对应的对数0.497149877671004

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

[讨论] 阶乘和圆周率

评分