gammal.c

linux下用PCMCIA无线网卡虚拟无线AP的程序源码

		- 2.68132716049382716049E-3L) * w
		+ 3.47222222222222222222E-3L) * w
		+ 8.33333333333333333333E-2L) * w
		+ 1.0L;
else
	w = 1.0L + w * polevll( w, STIR, 8 );
y = expl(x);
if( x > MAXSTIR )
	{ /* Avoid overflow in pow() */
	v = powl( x, 0.5L * x - 0.25L );
	y = v * (v / y);
	}
else
	{
	y = powl( x, x - 0.5L ) / y;
	}
y = SQTPI * y * w;
return( y );
}



long double gammal(x)
long double x;
{
long double p, q, z;
int i;

sgngaml = 1;
#ifdef NANS
if( isnanl(x) )
	return(NANL);
#endif
#ifdef INFINITIES
if(x == INFINITYL)
	return(INFINITYL);
#ifdef NANS
if(x == -INFINITYL)
	goto gamnan;
#endif
#endif
q = fabsl(x);

if( q > 13.0L )
	{
	if( q > MAXGAML )
		goto goverf;
	if( x < 0.0L )
		{
		p = floorl(q);
		if( p == q )
			{
gamnan:
#ifdef NANS
			mtherr( "gammal", DOMAIN );
			return (NANL);
#else
			goto goverf;
#endif
			}
		i = p;
		if( (i & 1) == 0 )
			sgngaml = -1;
		z = q - p;
		if( z > 0.5L )
			{
			p += 1.0L;
			z = q - p;
			}
		z = q * sinl( PIL * z );
		z = fabsl(z) * stirf(q);
		if( z <= PIL/MAXNUML )
			{
goverf:
#ifdef INFINITIES
			return( sgngaml * INFINITYL);
#else
			mtherr( "gammal", OVERFLOW );
			return( sgngaml * MAXNUML);
#endif
			}
		z = PIL/z;
		}
	else
		{
		z = stirf(x);
		}
	return( sgngaml * z );
	}

z = 1.0L;
while( x >= 3.0L )
	{
	x -= 1.0L;
	z *= x;
	}

while( x < -0.03125L )
	{
	z /= x;
	x += 1.0L;
	}

if( x <= 0.03125L )
	goto small;

while( x < 2.0L )
	{
	z /= x;
	x += 1.0L;
	}

if( x == 2.0L )
	return(z);

x -= 2.0L;
p = polevll( x, P, 7 );
q = polevll( x, Q, 8 );
return( z * p / q );

small:
if( x == 0.0L )
	{
	  goto gamnan;
	}
else
	{
	if( x < 0.0L )
		{
		x = -x;
		q = z / (x * polevll( x, SN, 8 ));
		}
	else
		q = z / (x * polevll( x, S, 8 ));
	}
return q;
}



/* A[]: Stirling's formula expansion of log gamma
 * B[], C[]: log gamma function between 2 and 3
 */


/* log gamma(x) = ( x - 0.5 ) * log(x) - x + LS2PI + 1/x A(1/x^2)
 * x >= 8
 * Peak relative error 1.51e-21
 * Relative spread of error peaks 5.67e-21
 */
#if UNK
static long double A[7] = {
 4.885026142432270781165E-3L,
-1.880801938119376907179E-3L,
 8.412723297322498080632E-4L,
-5.952345851765688514613E-4L,
 7.936507795855070755671E-4L,
-2.777777777750349603440E-3L,
 8.333333333333331447505E-2L,
};
#endif
#if IBMPC
static short A[] = {
0xd984,0xcc08,0x91c2,0xa012,0x3ff7, XPD
0x3d91,0x0304,0x3da1,0xf685,0xbff5, XPD
0x3bdc,0xaad1,0xd492,0xdc88,0x3ff4, XPD
0x8b20,0x9fce,0x844e,0x9c09,0xbff4, XPD
0xf8f2,0x30e5,0x0092,0xd00d,0x3ff4, XPD
0x4d88,0x03a8,0x60b6,0xb60b,0xbff6, XPD
0x9fcc,0xaaaa,0xaaaa,0xaaaa,0x3ffb, XPD
};
#endif
#if MIEEE
static long A[21] = {
0x3ff70000,0xa01291c2,0xcc08d984,
0xbff50000,0xf6853da1,0x03043d91,
0x3ff40000,0xdc88d492,0xaad13bdc,
0xbff40000,0x9c09844e,0x9fce8b20,
0x3ff40000,0xd00d0092,0x30e5f8f2,
0xbff60000,0xb60b60b6,0x03a84d88,
0x3ffb0000,0xaaaaaaaa,0xaaaa9fcc,
};
#endif

/* log gamma(x+2) = x B(x)/C(x)
 * 0 <= x <= 1
 * Peak relative error 7.16e-22
 * Relative spread of error peaks 4.78e-20
 */
#if UNK
static long double B[7] = {
-2.163690827643812857640E3L,
-8.723871522843511459790E4L,
-1.104326814691464261197E6L,
-6.111225012005214299996E6L,
-1.625568062543700591014E7L,
-2.003937418103815175475E7L,
-8.875666783650703802159E6L,
};
static long double C[7] = {
/* 1.000000000000000000000E0L,*/
-5.139481484435370143617E2L,
-3.403570840534304670537E4L,
-6.227441164066219501697E5L,
-4.814940379411882186630E6L,
-1.785433287045078156959E7L,
-3.138646407656182662088E7L,
-2.099336717757895876142E7L,
};
#endif
#if IBMPC
static short B[] = {
0x9557,0x4995,0x0da1,0x873b,0xc00a, XPD
0xfe44,0x9af8,0x5b8c,0xaa63,0xc00f, XPD
0x5aa8,0x7cf5,0x3684,0x86ce,0xc013, XPD
0x259a,0x258c,0xf206,0xba7f,0xc015, XPD
0xbe18,0x1ca3,0xc0a0,0xf80a,0xc016, XPD
0x168f,0x2c42,0x6717,0x98e3,0xc017, XPD
0x2051,0x9d55,0x92c8,0x876e,0xc016, XPD
};
static short C[] = {
/*0x0000,0x0000,0x0000,0x8000,0x3fff, XPD*/
0xaa77,0xcf2f,0xae76,0x807c,0xc008, XPD
0xb280,0x0d74,0xb55a,0x84f3,0xc00e, XPD
0xa505,0xcd30,0x81dc,0x9809,0xc012, XPD
0x3369,0x4246,0xb8c2,0x92f0,0xc015, XPD
0x63cf,0x6aee,0xbe6f,0x8837,0xc017, XPD
0x26bb,0xccc7,0xb009,0xef75,0xc017, XPD
0x462b,0xbae8,0xab96,0xa02a,0xc017, XPD
};
#endif
#if MIEEE
static long B[21] = {
0xc00a0000,0x873b0da1,0x49959557,
0xc00f0000,0xaa635b8c,0x9af8fe44,
0xc0130000,0x86ce3684,0x7cf55aa8,
0xc0150000,0xba7ff206,0x258c259a,
0xc0160000,0xf80ac0a0,0x1ca3be18,
0xc0170000,0x98e36717,0x2c42168f,
0xc0160000,0x876e92c8,0x9d552051,
};
static long C[21] = {
/*0x3fff0000,0x80000000,0x00000000,*/
0xc0080000,0x807cae76,0xcf2faa77,
0xc00e0000,0x84f3b55a,0x0d74b280,
0xc0120000,0x980981dc,0xcd30a505,
0xc0150000,0x92f0b8c2,0x42463369,
0xc0170000,0x8837be6f,0x6aee63cf,
0xc0170000,0xef75b009,0xccc726bb,
0xc0170000,0xa02aab96,0xbae8462b,
};
#endif

/* log( sqrt( 2*pi ) ) */
static long double LS2PI  =  0.91893853320467274178L;
#define MAXLGM 1.04848146839019521116e+4928L


/* Logarithm of gamma function */


long double lgaml(x)
long double x;
{
long double p, q, w, z, f, nx;
int i;

sgngaml = 1;
#ifdef NANS
if( isnanl(x) )
	return(NANL);
#endif
#ifdef INFINITIES
if( !isfinitel(x) )
	return(INFINITYL);
#endif
if( x < -34.0L )
	{
	q = -x;
	w = lgaml(q); /* note this modifies sgngam! */
	p = floorl(q);
	if( p == q )
		{
#ifdef INFINITIES
		mtherr( "lgaml", SING );
		return (INFINITYL);
#else
		goto loverf;
#endif
		}
	i = p;
	if( (i & 1) == 0 )
		sgngaml = -1;
	else
		sgngaml = 1;
	z = q - p;
	if( z > 0.5L )
		{
		p += 1.0L;
		z = p - q;
		}
	z = q * sinl( PIL * z );
	if( z == 0.0L )
		goto loverf;
/*	z = LOGPI - logl( z ) - w; */
	z = logl( PIL/z ) - w;
	return( z );
	}

if( x < 13.0L )
	{
	z = 1.0L;
	nx = floorl( x +  0.5L );
	f = x - nx;
	while( x >= 3.0L )
		{
		nx -= 1.0L;
		x = nx + f;
		z *= x;
		}
	while( x < 2.0L )
		{
		if( fabsl(x) <= 0.03125 )
			goto lsmall;
		z /= nx +  f;
		nx += 1.0L;
		x = nx + f;
		}
	if( z < 0.0L )
		{
		sgngaml = -1;
		z = -z;
		}
	else
		sgngaml = 1;
	if( x == 2.0L )
		return( logl(z) );
	x = (nx - 2.0L) + f;
	p = x * polevll( x, B, 6 ) / p1evll( x, C, 7);
	return( logl(z) + p );
	}

if( x > MAXLGM )
	{
loverf:
#ifdef INFINITIES
	return( sgngaml * INFINITYL );
#else
	mtherr( "lgaml", OVERFLOW );
	return( sgngaml * MAXNUML );
#endif
	}

q = ( x - 0.5L ) * logl(x) - x + LS2PI;
if( x > 1.0e10L )
	return(q);
p = 1.0L/(x*x);
q += polevll( p, A, 6 ) / x;
return( q );


lsmall:
if( x == 0.0L )
	goto loverf;
if( x < 0.0L )
	{
	x = -x;
	q = z / (x * polevll( x, SN, 8 ));
	}
else
	q = z / (x * polevll( x, S, 8 ));
if( q < 0.0L )
	{
	sgngaml = -1;
	q = -q;
	}
else
	sgngaml = 1;
q = logl( q );
return(q);
}