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[i]
                den = den * x + lanczos_den_coeffs[i]
        else:
            for i in range(len(lanczos_num_coeffs)):
                num = num / x + lanczos_num_coeffs[i]
                den = den / x + lanczos_den_coeffs[i]
        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.

It's not perfect, but you could use numerical integration to get the gamma function.

First transform the basic definition:

with the substitution u=exp(-t) to get an integral over a finite range [0,1]:

Numerically you would have to evaluate this with something like the mid-ordinate rule, to avoid the fact that the integrand becomes infinite at the lower limit. Alternatively, you could revert to the original form (with infinite upper limit) by truncating at some high, but finite, limit and putting an upper bound on the part omitted.

I found it helpful to reduce the basic integration to the range 1<=x<2 by using one of the inductive formulae

where necessary.

Obviously, the gamma function isn't defined at negative integers or zero.

import math

def tg(x):
    if x >= 2: return ( x - 1 ) * tg( x - 1 )     # Bring into argument range 1 <= x < 2
    if x <  1: return tg( x + 1 ) / x             #

    N = 10000
    result = 0
    for i in range( N ):
        result += math.pow( math.log( N / ( i + 0.5 ) ), x - 1 )     # mid-ordinate rule
    return result / N

for x in [ -2.5, -0.5, -0.1, 0.1, 0.3, 1, 2, 3.5, 10, 80 ]:
    print( x, tg(x), math.gamma(x), sep="  " )

Output (x, approximation, math.gamma)

-2.5  -0.9453033012805948  -0.9453087204829417
-0.5  -3.5448873798022307  -3.544907701811032
-0.1  -10.686015248516997  -10.686287021193193
0.1  9.513514929597076  9.513507698668732
0.3  2.991563821703771  2.991568987687591
1  1.0  1.0
2  1.0  1.0
3.5  3.3233319185645915  3.323350970447842
10  362880.0  362880.0
80  8.946182130782973e+116  8.946182130782976e+116