import numberdb.sage as numberdb
from sage.arith.misc import factorial
from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
from sage.rings.rational_field import QQ

R = PolynomialRing(QQ, 'x')
x = R.gen()

def mul(a, b, order):
    out = [R.zero() for _ in range(order + 1)]
    for i, ai in enumerate(a):
        for j, bj in enumerate(b[:order + 1 - i]):
            out[i + j] += ai * bj
    return out

def base_series(order):
    a = [QQ(0) for _ in range(order + 1)]
    a[0] = QQ(1)
    for r in range(2, order + 2):
        total = QQ(0)
        for m in range(2, r + 1):
            total += a[r - m] * QQ((-1)**(m + 1)) / QQ(factorial(m))
        a[r - 1] = -total
    return a

def binomial_x_plus_1(m):
    value = R.one()
    for j in range(m):
        value *= x + QQ(1 - j)
        value /= QQ(j + 1)
    return value

def stirling_polynomial(k):
    base = [R(c) for c in base_series(k)]
    base[0] = R.zero()
    series = [R.zero() for _ in range(k + 1)]
    series[0] = R.one()
    power = [R.zero() for _ in range(k + 1)]
    power[0] = R.one()
    for m in range(1, k + 1):
        power = mul(power, base, k)
        factor = binomial_x_plus_1(m)
        for j in range(k + 1):
            series[j] += factor * power[j]
    return series[k] * QQ(factorial(k))

stirling_polynomial(31)      # the next one after this table

The generator computes the ordinary power-series coefficients of $t/(1-e^{-t})$ from the identity $\frac{t}{1-e^{-t}}(1-e^{-t})=t$, then raises that series to the polynomial exponent $x+1$ over Sage's rational polynomial ring. Every division is made in Sage's rational field.