#!/usr/bin/env python
from math import gamma,e
# Continued Fraction Computation
# 6.5.31 Handbook of Mathematical Functions, page 263
# Recursive implementation
def upper_incomplete_gamma(a,x,d=0,iterations=100):
if d == iterations:
if ((d % 2) == 1):
return 1.0 # end iterations
else:
m = d/2
return x + (m-a)
if d == 0:
result = ((x**a) * (e**(-x)))/upper_incomplete_gamma(a,x,d=d+1)
return result
elif ((d % 2) == 1):
m = 1.0+((d-1.0)/2.0)
return x+ ((m-a)/(upper_incomplete_gamma(a,x,d=d+1)))
else:
m = d/2
return 1+(m/(upper_incomplete_gamma(a,x,d=d+1)))
# 6.5.31 Handbook of Mathematical Functions, page 263
# Recursive implementation
def upper_incomplete_gamma2(a,x,d=0,iterations=100):
if d == iterations:
return 1.0
if d == 0:
result = ((x**a) * (e**(-x)))/upper_incomplete_gamma2(a,x,d=d+1)
return result
else:
m = (d*2)-1
return (m-a)+x+ ((d*(a-d))/(upper_incomplete_gamma2(a,x,d=d+1)))
def lower_incomplete_gamma(a,x,d=0,iterations=100):
if d == iterations:
if ((d % 2) == 1):
return 1.0 # end iterations
else:
m = d/2
return x + (m-a)
if d == 0:
result = ((x**a) * (e**(-x)))/lower_incomplete_gamma(a,x,d=d+1)
return result
elif ((d % 2) == 1):
m = d - 1
n = (d-1.0)/2.0
return a + m - (((a+n)*x)/lower_incomplete_gamma(a,x,d=d+1))
else:
m = d-1
n = d/2.0
return a+m+((n*x)/(lower_incomplete_gamma(a,x,d=d+1)))
def lower_incomplete_gamma2(a,x):
return gamma(a)-upper_incomplete_gamma2(a,x)
def complimentary_incomplete_gamma(a,x):
return 1.0-upper_incomplete_gamma(a,x)
# Scipy name mappings
def gammainc(a,x):
return lower_incomplete_gamma(a,x)/gamma(a)
def gammaincc(a,x):
return upper_incomplete_gamma(a,x)/gamma(a)
