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))