该楼层疑似违规已被系统折叠 隐藏此楼查看此楼

// 素数定理li(x)计算值

#include

#include

#include

#include

#include "my_data.h"

typedef long double extended; // 定义扩展精度类型

// 欧拉常数 γ = 0.57721566490153286060651209...

const extended GAMMA = 0.5772156649015328606L;

const extended __li2 = 1.045163780117492785L;

extended fact(uint32);

extended Ei(extended,long,long *); //标准算法

extended li(extended,long,long *); //校对算法

void liDemo()

{

clock_t t1,t2;

long i, c1, c2;

long v1 = -100L;

extended d1 = 10.0L;

extended e1,e2,u1;

for (i = 1; i <= 200; i++)

{

u1 = logl(d1);

t2 = clock();

e1 = Ei(u1,v1,&c1);

e2 = li(d1,v1,&c2);

t1 = clock();

printf("10^%3ld = %.18Le  %.18Le",i,e1,e2);

printf("  %4ld  %4ld  %2ld\n",c1,c2,t1-t2);

d1 *= 10.0L;

}

} // 演示素数定理li(10^N)计算值(N<=200)

int main()

{

liDemo();

system("PAUSE");

return 0;

}

extended fact(uint32 n)

{

if (n == 0) return 1.0L;

else return n * fact(n - 1);

} // 计算n的浮点阶乘

extended Ei(extended u, long v, long * c)

// Ei(u) := γ + lnu + sum{}

// sum {1 to inf} u^n / n*n!

{

extended d0 = powl(10,v);

extended d1 = logl(u);

extended e0,e1,e2,f0;

*c = 0;e0 = 0.0L;

do {

f0 = e0; *c += 1;

e1 = expl(*c * d1);

e2 = *c * fact(*c);

e0 += e1 / e2;

} while ((e0 != f0) &&

(fabsl(e0 - f0) >= d0));

return GAMMA + d1 + e0;

} // 计算Ei(u)即li(e^u),存在误差

extended li(extended x, long v, long * c)

{

// li(x) := γ + ln(ln(x)) + sqrt(x) * expr1

// expr1 := sum{n:=1_to_inf}(-1^n-1)ln(x)^n

// /n!2^n-1*sum{k:=1_to_floor((n-1)/2)}1/(2k+1)

extended d0 = powl(10,v);

extended d1 = logl(logl(x));

extended e0,e1,e2,e3,f0;

*c = 0; e0 = 0.0L;

do

{

f0 = e0; *c += 1;

e1 = expl(*c * d1);

if (!(*c & 1)) e1 = -e1;

e2 = fact(*c)*powl(2,*c-1);

e3 = 0.0; long k;

for (k = 0; k <= (*c-1)>>1; k++)

{

e3 += (extended)1/((k<<1)+1);

}

e0 += e1 / e2 * e3;

} while ((e0 != f0) &&

(fabsl(e0 - f0) >= d0));

return GAMMA + d1 + e0 * sqrtl(x);

} //计算li(x)的一个更快的级数,存在误差