/*
	C program for floating point log Gamma function

	lgamma(x) computes the log of the absolute
	value of the Gamma function.
	The sign of the Gamma function is returned in the
	external quantity signgam.

	The coefficients for expansion around zero
	are #5243 from Hart & Cheney; for expansion
	around infinity they are #5404.

	Calls log, floor and sin.
*/

mesa double Libm←sin(), Libm←cos(), Libm←tan();
mesa double Libm←asin(), Libm←acos(), Libm←atan(), Libm←atan2();
mesa double Libm←sinh(), Libm←cosh(), Libm←tanh();
mesa double Libm←asinh(), Libm←acosh(), Libm←atanh();
mesa double Libm←log(), Libm←log1p();
mesa double Libm←exp(), Libm←frexp(), Libm←ldexp(), Libm←expm1(), Libm←pow();
mesa double Libm←fmod(), Libm←modf();
mesa double Libm←floor(), Libm←ceil(), Libm←rint();
mesa double Libm←cabs(), Libm←hypot();
mesa double Libm←sqrt(), Libm←cbrt(); 
mesa double Libm←j0(), Libm←j1(), Libm←jn();
mesa double Libm←fabs();
mesa int Libm←abs();

#define sin(x) (Libm←sin(x))
#define cos(x) (Libm←cos(x))
#define tan(x) (Libm←tan(x))
#define asin(x) (Libm←asin(x))
#define acos(x) (Libm←acos(x))
#define atan(x) (Libm←atan(x))
#define atan2(x,y) (Libm←atan2(x,y))
#define sinh(x) (Libm←sinh(x))
#define cosh(x) (Libm←cosh(x))
#define tanh(x) (Libm←tanh(x))
#define asinh(x) (Libm←asinh(x))
#define acosh(x) (Libm←acosh(x))
#define atanh(x) (Libm←atanh(x))
#define log(x) (Libm←log(x))
#define log1p(x) (Libm←log1p(x))
#define exp(x) (Libm←exp(x))
#define frexp(x,i) (Libm←frexp(x,i))
#define ldexp(x,e) (Libm←ldexp(x,e))
#define expm1(x) (Libm←expm1(x))
#define pow(x,y) (Libm←pow(x,y))
#define fmod(x,y) (Libm←fmod(x,y))
#define modf(d,i) (Libm←modf(d,i))
#define floor(d) (Libm←floor(d))
#define ceil(d) (Libm←ceil(d))
#define rint(x) (Libm←rint(x))
#define cabs(z) (Libm←cabs(z))
#define hypot(x,y) (Libm←hypot(x,y))
#define sqrt(x) (Libm←sqrt(x))
#define cbrt(x) (Libm←cbrt(x))
#define j0(a) (Libm←j0(a))
#define j1(a) (Libm←j1(a))
#define jn(n,x) (Libm←jn(n,x))
#define fabs(d) (Libm←fabs(d))
#define abs(i) (Libm←abs(i))

#define HUGE	1.701411733192644270e38
f

asm ("        export Libm");

int	signgam = 0;
static double goobie	= 0.9189385332046727417803297;	/* log(2*pi)/2 */
static double pi	= 3.1415926535897932384626434;

#define M 6
#define N 8
static double p1[] = {
	0.83333333333333101837e-1,
	-.277777777735865004e-2,
	0.793650576493454e-3,
	-.5951896861197e-3,
	0.83645878922e-3,
	-.1633436431e-2,
};
static double p2[] = {
	-.42353689509744089647e5,
	-.20886861789269887364e5,
	-.87627102978521489560e4,
	-.20085274013072791214e4,
	-.43933044406002567613e3,
	-.50108693752970953015e2,
	-.67449507245925289918e1,
	0.0,
};
static double q2[] = {
	-.42353689509744090010e5,
	-.29803853309256649932e4,
	0.99403074150827709015e4,
	-.15286072737795220248e4,
	-.49902852662143904834e3,
	0.18949823415702801641e3,
	-.23081551524580124562e2,
	0.10000000000000000000e1,
};

asm ("        exportproc ←lgamma, Libm");
double
lgamma(arg)
double arg;
{
	double log(), pos(), neg(), asym();

	signgam = 1.;
	if(arg <= 0.) return(neg(arg));
	if(arg > 8.) return(asym(arg));
	return(log(pos(arg)));
}

static double
asym(arg)
double arg;
{
	double log();
	double n, argsq;
	int i;

	argsq = 1./(arg*arg);
	for(n=0,i=M-1; i>=0; i--){
		n = n*argsq + p1[i];
	}
	return((arg-.5)*log(arg) - arg + goobie + n/arg);
}

static double
neg(arg)
double arg;
{
	double t;
	double log(), sin(), floor(), pos();

	arg = -arg;
     /*
      * to see if arg were a true integer, the old code used the
      * mathematically correct observation:
      * sin(n*pi) = 0 <=> n is an integer.
      * but in finite precision arithmetic, sin(n*PI) will NEVER
      * be zero simply because n*PI is a rational number.  hence
      *	it failed to work with our newer, more accurate sin()
      * which uses true pi to do the argument reduction...
      *	temp = sin(pi*arg);
      */
	t = floor(arg);
	if (arg - t  > 0.5e0)
	    t += 1.e0;				/* t := integer nearest arg */
	signgam = (int) (t - 2*floor(t/2));	/* signgam =  1 if t was odd, */
						/*            0 if t was even */
	signgam = signgam - 1 + signgam;	/* signgam =  1 if t was odd, */
						/*           -1 if t was even */
	t = arg - t;				/*  -0.5 <= t <= 0.5 */
	if (t < 0.e0) {
	    t = -t;
	    signgam = -signgam;
	}
	return(-log(arg*pos(arg)*sin(pi*t)/pi));
}

static double
pos(arg)
double arg;
{
	double n, d, s;
	register i;

	if(arg < 2.) return(pos(arg+1.)/arg);
	if(arg > 3.) return((arg-1.)*pos(arg-1.));

	s = arg - 2.;
	for(n=0,d=0,i=N-1; i>=0; i--){
		n = n*s + p2[i];
		d = d*s + q2[i];
	}
	return(n/d);
}