#########################################################################
#
# This script implements the log-gamma function and the incomplete
# gamma function.
#
# Author: Tom Loredo, circa Apr 2001. Fudged by Igor Volobouev
# to remove unnecessary dependence on Numeric and to increase
# precision of the gammln function.
#
#########################################################################

from math import *

#============= Exceptions ===============

__max_iters = 'Too many iterations: '

#============= Global constants ===============

__gammln_cof = (76.18009172947146, -86.50532032941677,
                24.01409824083091, -1.231739572450155,
                0.1208650973866179e-2, -0.5395239384953e-5)
__gammln_stp = 2.5066282746310005

#============= Gamma, Incomplete Gamma ===========

def gammln(xx):
	"""Logarithm of the gamma function."""
        # Based on the Lanczos approximation
	x = xx - 1.
	tmp = x + 5.5
	tmp = (x + 0.5)*log(tmp) - tmp
	ser = 1.000000000190015
	for j in range(6):
		x = x + 1.
		ser = ser + __gammln_cof[j]/x
	return tmp + log(__gammln_stp*ser)


def __gser(a, x, itmax=700, eps=3.e-7):
	"""Series approx'n to the incomplete gamma function."""
	gln = gammln(a)
	if (x < 0.):
		raise bad_arg, x
	if (x == 0.):
		return(0.)
	ap = a
	sum = 1. / a
	delta = sum
	n = 1
	while n <= itmax:
		ap = ap + 1.
		delta = delta * x / ap
		sum = sum + delta
		if (abs(delta) < abs(sum)*eps):
			return (sum * exp(-x + a*log(x) - gln), gln)
		n = n + 1
	raise __max_iters, str((abs(delta), abs(sum)*eps))


def __gcf(a, x, itmax=200, eps=3.e-7):
	"""Continued fraction approx'n of the incomplete gamma function."""
	gln = gammln(a)
	gold = 0.
	a0 = 1.
	a1 = x
	b0 = 0.
	b1 = 1.
	fac = 1.
	n = 1
	while n <= itmax:
		an = n
		ana = an - a
		a0 = (a1 + a0*ana)*fac
		b0 = (b1 + b0*ana)*fac
		anf = an*fac
		a1 = x*a0 + anf*a1
		b1 = x*b0 + anf*b1
		if (a1 != 0.):
			fac = 1. / a1
			g = b1*fac
			if (abs((g-gold)/g) < eps):
				return (g*exp(-x+a*log(x)-gln), gln)
			gold = g
		n = n + 1
	raise __max_iters, str(abs((g-gold)/g))


# The following function is the regularized incomplete gamma.
# For more info see the "Numerical Recipes" book.
def gammp(a, x):
	"""Incomplete gamma function."""
	if (x < 0. or a <= 0.):
		raise ValueError, (a, x)
	if (x < a+1.):
		return __gser(a,x)[0]
	else:
		return 1.-__gcf(a,x)[0]


# This is the incomplete gamma function complement
def gammq(a, x):
	"""Complementary incomplete gamma function."""
	if (x < 0. or a <= 0.):
		raise ValueError, repr((a, x))
	if (x < a+1.):
		return 1.-__gser(a,x)[0]
	else:
		return __gcf(a,x)[0]