cephes-old.c

随机数测试标准程序NIST

 *                   1       | |  -t  a-1
 *  igam(a,x)  =   -----     |   e   t   dt.
 *                  -      | |
 *                 | (a)    -
 *                           0
 *
 *
 * In this implementation both arguments must be positive.
 * The integral is evaluated by either a power series or
 * continued fraction expansion, depending on the relative
 * values of a and x.
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    IEEE      0,30       200000       3.6e-14     2.9e-15
 *    IEEE      0,100      300000       9.9e-14     1.5e-14
 */
/*                                                     igamc()
 *
 *     Complemented incomplete gamma integral
 *
 *
 *
 * SYNOPSIS:
 *
 * double a, x, y, igamc();
 *
 * y = igamc( a, x );
 *
 * DESCRIPTION:
 *
 * The function is defined by
 *
 *
 *  igamc(a,x)   =   1 - igam(a,x)
 *
 *                            inf.
 *                              -
 *                     1       | |  -t  a-1
 *               =   -----     |   e   t   dt.
 *                    -      | |
 *                   | (a)    -
 *                             x
 *
 *
 * In this implementation both arguments must be positive.
 * The integral is evaluated by either a power series or
 * continued fraction expansion, depending on the relative
 * values of a and x.
 *
 * ACCURACY:
 *
 * Tested at random a, x.
 *                a         x                      Relative error:
 * arithmetic   domain   domain     # trials      peak         rms
 *    IEEE     0.5,100   0,100      200000       1.9e-14     1.7e-15
 *    IEEE     0.01,0.5  0,100      200000       1.4e-13     1.6e-15
 */

/*
Cephes Math Library Release 2.0:  April, 1987
Copyright 1985, 1987 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

double igamc(double a, double x)
{
	double ans, ax, c, yc, r, t, y, z;
	double pk, pkm1, pkm2, qk, qkm1, qkm2;

	if( (x <= 0) || ( a <= 0) )
		return( 1.0 );

	if( (x < 1.0) || (x < a) )
		return( 1.e0 - igam(a,x) );

	ax = a * log(x) - x - lgam(a);
	if( ax < -MAXLOG ) {
		mtherr( "igamc", UNDERFLOW );
		return( 0.0 );
	}
	ax = exp(ax);

	/* continued fraction */
	y = 1.0 - a;
	z = x + y + 1.0;
	c = 0.0;
	pkm2 = 1.0;
	qkm2 = x;
	pkm1 = x + 1.0;
	qkm1 = z * x;
	ans = pkm1/qkm1;

	do {
		c += 1.0;
		y += 1.0;
		z += 2.0;
		yc = y * c;
		pk = pkm1 * z  -  pkm2 * yc;
		qk = qkm1 * z  -  qkm2 * yc;
		if( qk != 0 ) {
			r = pk/qk;
			t = fabs( (ans - r)/r );
			ans = r;
		}
		else
			t = 1.0;
		pkm2 = pkm1;
		pkm1 = pk;
		qkm2 = qkm1;
		qkm1 = qk;
		if( fabs(pk) > big ) {
			pkm2 *= biginv;
			pkm1 *= biginv;
			qkm2 *= biginv;
			qkm1 *= biginv;
		}
	} while( t > MACHEP );

	return( ans * ax );
}



/* left tail of incomplete gamma function:
 *
 *          inf.      k
 *   a  -x   -       x
 *  x  e     >   ----------
 *           -     -
 *          k=0   | (a+k+1)
 *
 */

double igam(double a, double x)
{
	double ans, ax, c, r;

	if( (x <= 0) || ( a <= 0) )
		return( 0.0 );

	if( (x > 1.0) && (x > a ) )
		return( 1.e0 - igamc(a,x) );

	/* Compute  x**a * exp(-x) / gamma(a)  */
	ax = a * log(x) - x - lgam(a);
	if( ax < -MAXLOG ) {
		mtherr( "igam", UNDERFLOW );
		return( 0.0 );
	}
	ax = exp(ax);

	/* power series */
	r = a;
	c = 1.0;
	ans = 1.0;

	do {
		r += 1.0;
		c *= x/r;
		ans += c;
	} while( c/ans > MACHEP );

	return( ans * ax/a );
}

/*                                                      gamma.c
 *
 *     Gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, gamma();
 * extern int sgngam;
 *
 * y = gamma( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns gamma function of the argument.  The result is
 * correctly signed, and the sign (+1 or -1) is also
 * returned in a global (extern) variable named sgngam.
 * This variable is also filled in by the logarithmic gamma
 * function lgam().
 *
 * Arguments |x| <= 34 are reduced by recurrence and the function
 * approximated by a rational function of degree 6/7 in the
 * interval (2,3).  Large arguments are handled by Stirling's
 * formula. Large negative arguments are made positive using
 * a reflection formula.
 *
 *
 * ACCURACY:
 *
 *                      Relative error:
 * arithmetic   domain     # trials      peak         rms
 *    DEC      -34, 34      10000       1.3e-16     2.5e-17
 *    IEEE    -170,-33      20000       2.3e-15     3.3e-16
 *    IEEE     -33,  33     20000       9.4e-16     2.2e-16
 *    IEEE      33, 171.6   20000       2.3e-15     3.2e-16
 *
 * Error for arguments outside the test range will be larger
 * owing to error amplification by the exponential function.
 *
 */
/*                                                      lgam()
 *
 *     Natural logarithm of gamma function
 *
 *
 *
 * SYNOPSIS:
 *
 * double x, y, lgam();
 * extern int sgngam;
 *
 * y = lgam( x );
 *
 *
 *
 * DESCRIPTION:
 *
 * Returns the base e (2.718...) logarithm of the absolute
 * value of the gamma function of the argument.
 * The sign (+1 or -1) of the gamma function is returned in a
 * global (extern) variable named sgngam.
 *
 * For arguments greater than 13, the logarithm of the gamma
 * function is approximated by the logarithmic version of
 * Stirling's formula using a polynomial approximation of
 * degree 4. Arguments between -33 and +33 are reduced by
 * recurrence to the interval [2,3] of a rational approximation.
 * The cosecant reflection formula is employed for arguments
 * less than -33.
 *
 * Arguments greater than MAXLGM return MAXNUM and an error
 * message.  MAXLGM = 2.035093e36 for DEC
 * arithmetic or 2.556348e305 for IEEE arithmetic.
 *
 *
 *
 * ACCURACY:
 *
 *
 * arithmetic      domain        # trials     peak         rms
 *    DEC     0, 3                  7000     5.2e-17     1.3e-17
 *    DEC     2.718, 2.035e36       5000     3.9e-17     9.9e-18
 *    IEEE    0, 3                 28000     5.4e-16     1.1e-16
 *    IEEE    2.718, 2.556e305     40000     3.5e-16     8.3e-17
 * The error criterion was relative when the function magnitude
 * was greater than one but absolute when it was less than one.
 *
 * The following test used the relative error criterion, though
 * at certain points the relative error could be much higher than
 * indicated.
 *    IEEE    -200, -4             10000     4.8e-16     1.3e-16
 *
 */

/*                                                      gamma.c */
/*      gamma function  */

/*
Cephes Math Library Release 2.2:  July, 1992
Copyright 1984, 1987, 1989, 1992 by Stephen L. Moshier
Direct inquiries to 30 Frost Street, Cambridge, MA 02140
*/

/* Gamma function computed by Stirling's formula.
 * The polynomial STIR is valid for 33 <= x <= 172.
 */
static double stirf(double x)
{
	double y, w, v;

	w = 1.0/x;
	w = 1.0 + w * polevl( w, STIR, 4 );
	y = exp(x);
	if( x > MAXSTIR ) { /* Avoid overflow in pow() */
		v = pow( x, 0.5 * x - 0.25 );
		y = v * (v / y);
	}
	else
		y = pow( x, x - 0.5 ) / y;
	y = SQTPI * y * w;

	return( y );
}



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

	sgngam = 1;
#ifdef NANS
	if( isnan(x) )
		return(x);
#endif

#ifdef INFINITIES
 #ifdef NANS
	if( x == INFINITY )
		return(x);
	if( x == -INFINITY )
		return(NAN);
 #else
	if( !isfinite(x) )
		return(x);
 #endif
#endif
	q = fabs(x);

	if( q > 33.0 ) {
		if( x < 0.0 ) {
			p = floor(q);
			if( p == q ) {
#ifdef NANS
				gamnan:
				mtherr( "gamma", DOMAIN );
				return (NAN);
#else
				goto goverf;
#endif
			}
			i = (int)p;
			if( (i & 1) == 0 )
				sgngam = -1;
			z = q - p;
			if( z > 0.5 ) {
				p += 1.0;
				z = q - p;
			}
			z = q * sin( PI * z );
			if( z == 0.0 ) {
#ifdef INFINITIES
				return( sgngam * INFINITY);
#else
				goverf:
				mtherr( "gamma", OVERFLOW );
				return( sgngam * MAXNUM);
#endif
			}
			z = fabs(z);
			z = PI/(z * stirf(q) );
		}
	else {
		z = stirf(x);
	}

	return( sgngam * z );
	}

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

	while( x < 0.0 ) {
		if( x > -1.E-9 )
			goto small;
		z /= x;
		x += 1.0;
	}

	while( x < 2.0 ) {
		if( x < 1.e-9 )
			goto small;
		z /= x;
		x += 1.0;
	}

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

	x -= 2.0;
	p = polevl( x, P, 6 );
	q = polevl( x, Q, 7 );
	return( z * p / q );

small:
	if( x == 0.0 ) {
#ifdef INFINITIES
 #ifdef NANS
		goto gamnan;
 #else
		return( INFINITY );
 #endif
#else
		mtherr( "gamma", SING );
		return( MAXNUM );
#endif
	}
else

	return( z/((1.0 + 0.5772156649015329 * x) * x) );
}



/* A[]: Stirling's formula expansion of log gamma
 * B[], C[]: log gamma function between 2 and 3
 */
#ifdef UNK
static double A[] = {
	 8.11614167470508450300E-4,
	-5.95061904284301438324E-4,
	 7.93650340457716943945E-4,
	-2.77777777730099687205E-3,
	 8.33333333333331927722E-2
};
static double B[] = {
	-1.37825152569120859100E3,
	-3.88016315134637840924E4,
	-3.31612992738871184744E5,
	-1.16237097492762307383E6,
	-1.72173700820839662146E6,
	-8.53555664245765465627E5
};
static double C[] = {
	/* 1.00000000000000000000E0, */
	-3.51815701436523470549E2,
	-1.70642106651881159223E4,
	-2.20528590553854454839E5,
	-1.13933444367982507207E6,
	-2.53252307177582951285E6,
	-2.01889141433532773231E6
};
/* log( sqrt( 2*pi ) ) */
static double LS2PI  =  0.91893853320467274178;
#define MAXLGM 2.556348e305
#endif

#ifdef DEC
static unsigned short A[] = {
	0035524,0141201,0034633,0031405,
	0135433,0176755,0126007,0045030,
	0035520,0006371,0003342,0172730,
	0136066,0005540,0132605,0026407,
	0037252,0125252,0125252,0125132
};
static unsigned short B[] = {
	0142654,0044014,0077633,0035410,
	0144027,0110641,0125335,0144760,
	0144641,0165637,0142204,0047447,
	0145215,0162027,0146246,0155211,
	0145322,0026110,0010317,0110130,
	0145120,0061472,0120300,0025363
};
static unsigned short C[] = {
	/*0040200,0000000,0000000,0000000*/
	0142257,0164150,0163630,0112622,
	0143605,0050153,0156116,0135272,
	0144527,0056045,0145642,0062332,
	0145213,0012063,0106250,0001025,
	0145432,0111254,0044577,0115142,
	0145366,0071133,0050217,0005122
};
/* log( sqrt( 2*pi ) ) */
static unsigned short LS2P[] = {040153,037616,041445,0172645,};
#define LS2PI *(double *)LS2P
#define MAXLGM 2.035093e36
#endif

#ifdef IBMPC
static unsigned short A[] = {
	0x6661,0x2733,0x9850,0x3f4a,
	0xe943,0xb580,0x7fbd,0xbf43,