/*
 * (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
 * See the copyright notice in the ACK home directory, in the file "Copyright".
 *
 * Author: Ceriel J.H. Jacobs
 */

/* $Header$ */

#include <math.h>
#include <errno.h>

extern int errno;

static double
smallpos_gamma(x)
	double x;
{
	/*	Approximation of gamma function using
		gamma(x) = P(x-2) / Q(x-2) for x in [2,3]
	*/
	/* Hart & Cheney # 5251 */

	static double p[11] = {
		-0.2983543278574342138830437659e+06,
		-0.2384953970018198872468734423e+06,
		-0.1170494760121780688403854445e+06,
		-0.3949445048301571936421824091e+05,
		-0.1046699423827521405330650531e+05,
		-0.2188218110071816359394795998e+04,
		-0.3805112208641734657584922631e+03,
		-0.5283123755635845383718978382e+02,
		-0.6128571763704498306889428212e+01,
		-0.5028018054416812467364198750e+00,
		-0.3343060322330595274515660112e-01
	};

	static double q[9] = {
		-0.2983543278574342138830438524e+06,
		-0.1123558608748644911342306408e+06,
		 0.5332716689118142157485686311e+05,
		 0.8571160498907043851961147763e+04,
		-0.4734865977028211706556819770e+04,
		 0.1960497612885585838997039621e+03,
		 0.1257733367869888645966647426e+03,
		-0.2053126153100672764513929067e+02,
		 0.1000000000000000000000000000e+01
	};

	double result = 1.0;

	while (x > 3) {
		x -= 1.0;
		result *= x;
	}
	while (x < 2) {
		result /= x;
		x += 1.0;
	}

	x -= 2.0;

	return result * POLYNOM10(x, p) / POLYNOM8(x, q);
}

#define log_sqrt_2pi 0.91893853320467274178032973640561763

int	signgam;

static double
bigpos_loggamma(x)
	double x;
{
	/*	computes the log(gamma(x)) function for big arguments
		using the Stirling form
		  log(gamma(x)) = (x - 0.5)log(x) - x + log(sqrt(2*pi)) + fi(x)
		where fi(x) = (1/x)*P(1/(x*x))/Q(1/(x*x)) for x in [12,1000]
	*/
	/* Hart & Cheney # 5468 */

	static double p[4] = {
		 0.12398282342474941538685913e+00,
		 0.67082783834332134961461700e+00,
		 0.64507302912892202513890000e+00,
		 0.66662907040200752600000000e-01
	};

	static double q[4] = {
		 0.14877938810969929846815600e+01,
		 0.80995271894897557472821400e+01,
		 0.79966911236636441947720000e+01,
		 0.10000000000000000000000000e+01
	};

	double rsq = 1.0/(x*x);
	extern double log();

	return (x-0.5)*log(x)-x+log_sqrt_2pi+POLYNOM3(rsq, p)/(x*POLYNOM3(rsq, q));
}

static double
neg_loggamma(x)
	double x;
{
	/*	compute the log(gamma(x)) function for negative values of x,
		using the rule:
			-x*gamma(x)*gamma(-x) = pi/sin(z*pi)
	*/
	extern double sin(), log();
	double sinpix;

	x = -x;
	sinpix = sin(M_PI * x);
	if (sinpix == 0.0) {
		errno = EDOM;
		return HUGE;
	}
	if (sinpix < 0) sinpix = -sinpix;
	else signgam = -1;
	return log(M_PI/(x * smallpos_gamma(x) * sinpix));
}

double
gamma(x)
	double x;
{
	/*	Wrong name; Actually computes log(gamma(x))
	*/
	extern double log();

	signgam = 1;
	if (x <= 0) {
		return neg_loggamma(x);
	}
	if (x > 12.0) {
		return bigpos_loggamma(x);
	}
	return log(smallpos_gamma(x));
}