[Python] Calculate gamma functions in python without external packages (besides 'math')

I want to make a python script to put in the TI-84, that calculates gamma functions without any outside packages like scipy, numpy. I can import a select few math packages as they are supported by the TI-84 calculator, but the other functions from math such as gamma aren't supported. I also know that the calculator has it's own functions built in to calculate gamma. I have some code that works, but it isn't always accurate especially with the more complex functions. I had some help from AI. And I am calculating for t distributions if that helps

from math import pi, sin, pow, exp, floor
def isfinite(x):
return not (x == float('inf') or x == -float('inf') or x != x)
def abs(x):
return x if x >= 0 else -x
def nan():
return float('nan')
def inf():
return float('inf')
def tgamma(x):
# Constants used in the Lanczos approximation
lanczos_g = 6.024680040776729583740234375
lanczos_g_minus_half = lanczos_g - 0.5
lanczos_num_coeffs = [
23531376880.410759688572007674451636754734846804940,
42919803642.649098768957899047001988850926355848959,
35711959237.355668049440185451547166705960488635843,
17921034426.037209699919755754458931112671403265390,
6039542586.3520280050642916443072979210699388420708,
1439720407.3117216736632230727949123939715485786772,
248874557.86205415651146038641322942321632125127801,
31426415.585400194380614231628318205362874684987640,
2876370.6289353724412254090516208496135991145378768,
186056.26539522349504029498971604569928220784236328,
8071.6720023658162106380029022722506138218516325024,
210.82427775157934587250973392071336271166969580291,
2.5066282746310002701649081771338373386264310793408
]
lanczos_den_coeffs = [
0.0, 39916800.0, 120543840.0, 150917976.0, 105258076.0, 45995730.0,
13339535.0, 2637558.0, 357423.0, 32670.0, 1925.0, 66.0, 1.0
]
gamma_integral = [
1.0, 1.0, 2.0, 6.0, 24.0, 120.0, 720.0, 5040.0, 40320.0, 362880.0,
3628800.0, 39916800.0, 479001600.0, 6227020800.0, 87178291200.0,
1307674368000.0, 20922789888000.0, 355687428096000.0,
6402373705728000.0, 121645100408832000.0, 2432902008176640000.0,
51090942171709440000.0, 1124000727777607680000.0,
]
def lanczos_sum(x):
num = 0.0
den = 0.0
if x < 5.0:
for i in range(len(lanczos_num_coeffs) - 1, -1, -1):
num = num * x + lanczos_num_coeffs
den = den * x + lanczos_den_coeffs
else:
for i in range(len(lanczos_num_coeffs)):
num = num / x + lanczos_num_coeffs
den = den / x + lanczos_den_coeffs
return num / den
if not isfinite(x):
if x != x or x > 0.0:
return x
else:
return nan()
if x == 0.0:
return inf()
if x == floor(x):
if x < 0.0:
return nan()
if x <= len(gamma_integral):
return gamma_integral[int(x) - 1]
absx = abs(x)
if absx < 1e-20:
return 1.0 / x
if absx > 200.0:
if x < 0.0:
return 0.0 / sin(pi * x)
else:
return inf()
y = absx + lanczos_g_minus_half
if absx > lanczos_g_minus_half:
q = y - absx
z = q - lanczos_g_minus_half
else:
q = y - lanczos_g_minus_half
z = q - absx
z = z * lanczos_g / y
if x < 0.0:
r = -pi / sin(pi * absx) / absx * exp(y) / lanczos_sum(absx)
r -= z * r
if absx < 140.0:
r /= pow(y, absx - 0.5)
else:
sqrtpow = pow(y, absx / 2.0 - 0.25)
r /= sqrtpow
r /= sqrtpow
else:
r = lanczos_sum(absx) / exp(y)
r += z * r
if absx < 140.0:
r *= pow(y, absx - 0.5)
else:
sqrtpow = pow(y, absx / 2.0 - 0.25)
r *= sqrtpow
r *= sqrtpow
return r
def tcdf(x, v):
"""
x (float): The t-value.
v (float): The degrees of freedom.
"""
if v <= 0:
return nan()
if x == 0:
return 0.5
beta = tgamma(v / 2.0) * tgamma(0.5) / tgamma((v + 1) / 2.0)
integral = (1 + (x ** 2) / v) ** (-(v + 1) / 2.0)
if x > 0:
return 1 - 0.5 * beta * integral
else:
return 0.5 * beta * integral
print(tcdf(0.345, 30))
>>>0.7829983071979724


This is the code that I have, but when it runs a complex function it usually outputs something, but it is sometimes off by a bit compared to the calculator, or just completely wrong. For example:

from scipy.stats import t
print(t.cdf(0.345, 30))
>>>0.6337492192570378


As you can see the out puts are clearly different, and the scipy is way more accurate as it is what is calculated on the TI-84 Plus CE. If you could help fix the code, or explain how the gamma function works in a way that is pretty simple.

Continue reading...