#!/usr/local/bin/pythonw
# coding: iso-8859-1
# Gamma.py, copyright 2009 by Amos Newcombe.
# 
# This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by the Free Software Foundation, either version 3 of the License, or (at your option) any later version.
# 
# This program is distributed in the hope that it will be useful, but WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License for more details.
# 
# You should have received a copy of the GNU General Public License along with this program.  If not, see <http://www.gnu.org/licenses/>.

'''the Gamma functions: complete and incomplete
'''

from operator import add, div
from math import log, exp, fabs
# import from cmath to get complex version of these functions

#	Gamma function: Gamma(n) = (n-1)!

def LnGamma(x):
	'''Returns ln(Gamma(x)).
	
	Follows _Numerical Recipes in C_, pp 213-214, but fixes a couple of errors.
	"You can see that this is sort of a take-off on Stirling's approximation, but with a series of corrections that take into account the first few poles in the left complex plane."
	'''
	if x <= 0: raise ValueError, 'Gamma(x): x = %s <= 0' % x
	x = x-1	# 1st error: this step inserted
	xT = x+5.5
	xT = xT - (x+0.5)*log(xT)
	lxCof = [
		 76.18009172947146,
		-86.50532032941677,
		 24.01409824083091,
		 -1.231739572450155,
		  0.1208650973866179e-2,
		 -0.5395239384953e-5
	]
	lxDenom = [x+n+1 for n in range(len(lxCof))]
	lx = map(div, lxCof, lxDenom)
	xSer = 1.000000000190015 + reduce(add, lx)
	return -xT + log(2.5066282746310005 * xSer)	# 2nd error: expression corrected

def Gamma(x):
	'''Returns Gamma(x).
	
	If x is a positive integer then Gamma(x) = (x-1)!
	
	>>> [int(round(Gamma(n))) for n in range(1,12)]
	[1, 1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800]
	'''
	return exp(LnGamma(x))

#	Incomplete gamma function

def IGamma(a, x):
	'''Returns the incomplete gamma function P(a, x).
	
	Follows _Numerical Recipes in C_, pp 216-219, section 6.2.
	'''
	if a <= 0: raise ValueError, 'IGamma(a, x): (a == %s) <= 0' % a
	if x <  0: raise ValueError, 'IGamma(a, x): (x == %s) <  0' % x
	if x < a+1.: return IGammaSeries(a, x)
	else       : return IGammaCF(a, x)

def IGammaSeries(a, x):
	if x == 0: return 0.
	epsilon = 3.0e-7
	xA = a
	xAMax = 101 + xA
	xTerm = xSum = 1.0/a
	while xA < xAMax:
		xA = xA + 1
		xTerm = xTerm * x / xA
		xSum = xSum + xTerm
		if abs(xTerm) < abs(xSum) * epsilon:
			return xSum * exp(-x + a*log(x) - LnGamma(a))
	raise ValueError, 'IGammaSeries(%s, %s) fails' % (a, x)

def IGammaCF(a, x):
	iMax = 101
	epsilon = 3.0e-7
	xTiny = 1.0e-30
	b = x + 1.0 - a
	c = 1.0 / xTiny
	h = d = 1.0 / b
	i = 1
	while i < iMax:
		an = -i*(i-a)
		b = b + 2.0
		d = an * d + b
		if abs(d) < xTiny: d = xTiny
		c = b + an / c
		if abs(c) < xTiny: c = xTiny
		d = 1.0/d
		t = d * c
		h = h * t
		if abs(t - 1.0) < epsilon:
			return 1.0 - h * exp(-x + a*log(x) - LnGamma(a))
		i = i+1
	raise ValueError, 'IGammaCF(%s, %s) fails' % (a, x)

def InverseIGamma(a, y, xEpsilon=1e-7):
# 	print 'InverseIGamma(%g, %g)' % (y, xEpsilon)
	iMax = 101
	a = float(a)
	x = a - 1.
	if x <= 0.: x = xEpsilon
	i = 1
	while i < iMax:
		d = y - IGamma(a, x)
		if fabs(d) < xEpsilon: break
		m = pow(x, a-1.0) * exp(-x) / Gamma(a)
		#print x, d, m
		x = x + d / m
		if x <= 0.: x = xEpsilon
		i = i+1
	return x

if __name__ == '__main__':

	import doctest
	tp = doctest.testmod()
	print '%d tests failed out of %d' % tp