RosettaCodeData/Task/Gamma-function/Stata/gamma-function-2.stata

32 lines
589 B
Plaintext

from mpmath import mp
mp.dps = 50
def gamma_coef(n):
a = [mp.mpf(1), mp.mpf(mp.euler)]
for k in range(3, n + 1):
s = sum((-1)**j * mp.zeta(j) * a[k - j - 1] for j in range(2, k))
a.append((s - a[1] * a[k - 2]) / (1 - k * a[0]))
return a
def horner(a, x):
y = 0
for s in reversed(a):
y = y * x + s
return y
gc = gamma_coef(30)
def gamma_approx(x):
y = mp.mpf(1)
while x < 0.5:
y /= x
x += 1
while x > 1.5:
x -= 1
y *= x
return y / horner(gc, x - 1)
for x in gc:
print(mp.nstr(x, 25))