gamma_quad(t) = quadgk(x -> x^(t-1)*exp(-x), zero(t), inf(t), reltol=100*eps(t))[1]
# generalize to handle arrays of inputs:
@vectorize_1arg Real gamma_quad