#ifndef lint
static char sccsid[] = "@(#)lgamma.c	1.2	(ucb.beef)	3/1/89";
#endif	/* !defined(lint) */
/* 
 * This routine calculates the log(GAMMA) function for a positive real
 *   argument x.  Computation is based on an algorithm outlined in
 *   references 1 and 2.  The program uses rational functions that
 *   theoretically approximate log(GAMMA) to at least 18 significant
 *   decimal digits.  The approximation for x > 12 is from reference
 *   3, while approximations for x < 12.0 are similar to those in
 *   reference 1, but are unpublished.  The accuracy achieved depends
 *   on the arithmetic system, the compiler, the intrinsic functions,
 *   and proper selection of the machine-dependent constants.
 * 
 **********************************************************************
 **********************************************************************
 * 
 * Explanation of machine-dependent constants
 * 
 * beta   - radix for the floating-point representation
 * maxexp - the smallest positive power of beta that overflows
 * XBIG   - largest argument for which LN(GAMMA(x)) is representable
 *          in the machine, i.e., the solution to the equation
 *                  LN(GAMMA(XBIG)) = beta**maxexp
 * XINF   - largest machine representable floating-point number;
 *          approximately beta**maxexp.
 * EPS    - The smallest positive floating-point number such that
 *          1.0+EPS > 1.0
 * FRTBIG - Rough estimate of the fourth root of XBIG
 * 
 * 
 *     Approximate values for some important machines are:
 * 
 *                           beta      maxexp         XBIG
 * 
 * CRAY-1        (S.P.)        2        8191       9.62E+2461
 * Cyber 180/855
 *   under NOS   (S.P.)        2        1070       1.72E+319
 * IEEE (IBM/XT,
 *   SUN, etc.)  (S.P.)        2         128       4.08E+36
 * IEEE (IBM/XT,
 *   SUN, etc.)  (D.P.)        2        1024       2.55D+305
 * IBM 3033      (D.P.)       16          63       4.29D+73
 * VAX D-Format  (D.P.)        2         127       2.05D+36
 * VAX G-Format  (D.P.)        2        1023       1.28D+305
 *
 * 
 *                           XINF        EPS        FRTBIG
 * 
 * CRAY-1        (S.P.)   5.45E+2465   7.11E-15    3.13E+615
 * Cyber 180/855
 *   under NOS   (S.P.)   1.26E+322    3.55E-15    6.44E+79
 * IEEE (IBM/XT,
 *   SUN, etc.)  (S.P.)   3.40E+38     1.19E-7     1.42E+9
 * IEEE (IBM/XT,
 *   SUN, etc.)  (D.P.)   1.79D+308    2.22D-16    2.25D+76
 * IBM 3033      (D.P.)   7.23D+75     2.22D-16    2.56D+18
 * VAX D-Format  (D.P.)   1.70D+38     1.39D-17    1.20D+9
 * VAX G-Format  (D.P.)   8.98D+307    1.11D-16    1.89D+76
 * 
 ***************************************************************
 ***************************************************************
 * 
 * Error returns
 * 
 *  The program returns the value XINF for x <= 0.0 or when
 *     overflow would occur.  The computation is believed to 
 *     be free of underflow and overflow.
 * 
 * 
 * Intrinsic functions required are:
 * 
 *      log
 * 
 * 
 * References:
 * 
 *  1) W. J. Cody and K. E. Hillstrom, 'Chebyshev Approximations for
 *     the Natural Logarithm of the Gamma Function,' Math. Comp. 21,
 *     1967, pp. 198-203.
 * 
 *  2) K. E. Hillstrom, ANL/AMD Program ANLC366S, DGAMMA/DLGAMA, May,
 *     1969.
 * 
 *  3) Hart, Et. Al., Computer Approximations, Wiley and sons, New
 *     York, 1968.
 * 
 * 
 *  Authors: W. J. Cody and L. Stoltz
 *           Argonne National Laboratory
 * 
 *  Latest modification: June 16, 1988
 */
#if defined(__STDC__) || defined(__GNUC__)
extern double log(double);
#else
extern double log();
#endif
					/* Machine dependent parameters */
#if defined(vax) || defined(tahoe)
#define	XBIG	(double)2.05e36
#define	XINF	(double)1.70e38
#define	EPS	(double)1.39e-17
#define	FRTBIG	(double)1.20e9
#else	/* defined(vax) || defined(tahoe) */
#define	XBIG	(double)2.55e305
#define	XINF	(double)1.79e308
#define	EPS	(double)2.22e-16
#define	FRTBIG	(double)2.25e76
#endif	/* defined(vax) || defined(tahoe) */
					/* Mathematical constants */
#define	ONE	(double)1
#define	HALF	(double)0.5
#define	TWELVE	(double)12
#define	ZERO	(double)0
#define	FOUR	(double)4
#define	THRHAL	(double)1.5
#define	TWO	(double)2
#define	PNT68	(double)0.6796875
#define	SQRTPI	(double)0.9189385332046727417803297

/*
 * Numerator and denominator coefficients for rational minimax Approximation
 * over (0.5,1.5).
 */
static double D1 = -5.772156649015328605195174e-1;
static double P1[] = {
	4.945235359296727046734888e0,
	2.018112620856775083915565e2,
	2.290838373831346393026739e3,
	1.131967205903380828685045e4,
	2.855724635671635335736389e4,
	3.848496228443793359990269e4,
	2.637748787624195437963534e4,
	7.225813979700288197698961e3,
};
static double Q1[] = {
	6.748212550303777196073036e1,
	1.113332393857199323513008e3,
	7.738757056935398733233834e3,
	2.763987074403340708898585e4,
	5.499310206226157329794414e4,
	6.161122180066002127833352e4,
	3.635127591501940507276287e4,
	8.785536302431013170870835e3,
};

/*
 * Numerator and denominator coefficients for rational minimax Approximation
 * over (1.5,4.0).
 */
static double D2 = 4.227843350984671393993777e-1;
static double P2[] = {
	4.974607845568932035012064e0,
	5.424138599891070494101986e2,
	1.550693864978364947665077e4,
	1.847932904445632425417223e5,
	1.088204769468828767498470e6,
	3.338152967987029735917223e6,
	5.106661678927352456275255e6,
	3.074109054850539556250927e6,
};
static double Q2[] = {
	1.830328399370592604055942e2,
	7.765049321445005871323047e3,
	1.331903827966074194402448e5,
	1.136705821321969608938755e6,
	5.267964117437946917577538e6,
	1.346701454311101692290052e7,
	1.782736530353274213975932e7,
	9.533095591844353613395747e6,
};

/*
 * Numerator and denominator coefficients for rational minimax Approximation
 * over (4.0,12.0).
 */
static double D4 = 1.791759469228055000094023e0;
static double P4[] = {
	1.474502166059939948905062e4,
	2.426813369486704502836312e6,
	1.214755574045093227939592e8,
	2.663432449630976949898078e9,
	2.940378956634553899906876e10,
	1.702665737765398868392998e11,
	4.926125793377430887588120e11,
	5.606251856223951465078242e11,
};
static double Q4[] = {
	2.690530175870899333379843e3,
	6.393885654300092398984238e5,
	4.135599930241388052042842e7,
	1.120872109616147941376570e9,
	1.488613728678813811542398e10,
	1.016803586272438228077304e11,
	3.417476345507377132798597e11,
	4.463158187419713286462081e11,
};

/*
 * Coefficients for minimax approximation over (12, INF).
 */
static double C[] = {
	-1.910444077728e-03,
	 8.4171387781295e-04,
	-5.952379913043012e-04,
	 7.93650793500350248e-04,
	-2.777777777777681622553e-03,
	 8.333333333333333331554247e-02,
	 5.7083835261e-03,
};

double
#if defined(__STDC__) || defined(__GNUC__)
lgamma(double x)
#else
lgamma(x)
double x;
#endif
{
	register i;
	double res,corr,xden,xnum,dtmp;

#define	y x
	if (y > ZERO && y <= XBIG) {
		if (y <= EPS)
			res = -log(y);
		else if (y <= THRHAL) {			/*  EPS < x <= 1.5 */
#define	xm1 dtmp
			if (y < PNT68) {
				corr = -log(y);
				xm1 = y;
			}
			else {
				corr = ZERO;
				xm1 = y-HALF; xm1 -= HALF;
			}
			if (y <= HALF || y >= PNT68) {
				xden = ONE;
				xnum = ZERO;
				for (i = 0; i <= 7; i++) {
					xnum = xnum*xm1+P1[i];
					xden = xden*xm1+Q1[i];
				}
				res = xnum/xden; res = corr+xm1*(D1+xm1*res);
#undef	xm1
			}
			else {
#define	xm2 dtmp
				xm2 = y-HALF; xm2 -= HALF;
				xden = ONE;
				xnum = ZERO;
				for (i = 0; i <= 7; i++) {
					xnum = xnum*xm2+P2[i];
					xden = xden*xm2+Q2[i];
				}
				res = xnum/xden; res = corr+xm2*(D2+xm2*res);
			}
		}
		else if (y <= FOUR) {			/*  1.5 < x <= 4.0 */
			xm2 = y-TWO;
			xden = ONE;
			xnum = ZERO;
			for (i = 0; i <= 7; i++) {
				xnum = xnum*xm2+P2[i];
				xden = xden*xm2+Q2[i];
			}
			res = xnum/xden; res = xm2*(D2+xm2*res);
#undef	xm2
		}
		else if (y <= TWELVE) {			/*  4.0 < x <= 12.0 */
#define	xm4 dtmp
			xm4 = y-FOUR;
			xden = -ONE;
			xnum = ZERO;
			for (i = 0; i <= 7; i++) {
				xnum = xnum*xm4+P4[i];
				xden = xden*xm4+Q4[i];
			}
			res = xnum/xden; res = D4+xm4*res;
#undef	xm4
		}
		else {					/*  x >= 12.0 */
			res = ZERO;
			if (y <= FRTBIG) {
				res = C[6];
#define	ysq dtmp
				ysq = y*y;
				for (i = 0; i <= 5; i++)
					res = res/ysq+C[i];
#define	ysq dtmp
			}
			res /= y;
			corr = log(y);
			res += SQRTPI; res -= HALF*corr;
			res += y*(corr-ONE);
		}
#undef	y
	}
	else					/*  Return for bad arguments */
		res = XINF;
	return res;
}