/* lgamma.c */

/* Computes logarithm of gamma function (generalized factorial). Basic
* formula is given in Abramowitz and Stegun 6.1.40. Algorithm folows 
* that of IBM's Scientific Subroutine Package (FORTRAN). 
* S E Peterson 3/92 */

#define _LIBM_SOURCE

#include <math.h>
#include <stdlib.h>

double lgamma(x)
double x;
{
	double z,b,lgm,ans;
	static double coef[12] = {8.3333333333333333333E-2,
						  -2.7777777777777777778E-3,
						   7.9365079365079365079E-4,
						  -5.9523809523809523810E-4,
						   8.4175084175084175084E-4,
						  -1.9175269175269175269E-3,
						   6.4102564102564102564E-3,
						  -2.9550653594771241830E-2,	
						   1.7964437236883057316E-1,
						  -1.3924322169059011164E0,
						   1.3402864044168391994E1,
						  -1.5684828462600201731E2};
	int i;
	double term = 1.0;
	z=x;
	if (x<18.0){
		for(i=1;i<=((int)(19-x));i++){
			term *= z;
			z += 1.0;
		}
	}		
	b = 1.0/(z*z);
	ans =  coef[11];
	for (i = 10; i >= 0; i--)
		ans = coef[i] + b * ans;
	ans /= z;	
	lgm = (z - 0.5)*log(z) - z + 0.91893853320467274178 - log(term) + ans;		
	
	return (lgm);

}