  The generalized Stieltjes constants $\gamma\_n(v)$ are, up to a simple
scaling factor, the Laurent series coefficients of the Hurwitz zeta function
$\zeta(s,v)$ about its unique pole $s = 1$. In this work, we devise an
efficient algorithm to compute these constants to arbitrary precision with
rigorous error bounds, for the first time achieving this with low complexity
with respect to the order~$n$. Our computations are based on an integral
representation with a hyperbolic kernel that decays exponentially fast. The
algorithm consists of locating an approximate steepest descent contour and then
evaluating the integral numerically in ball arithmetic using the Petras
algorithm with a Taylor expansion for bounds near the saddle point. An
implementation is provided in the Arb library. We can, for example, compute
$\gamma\_n(1)$ to 1000 digits in a minute for any $n$ up to $n=10^{100}$. We
also provide other interesting integral representations for $\gamma\_n(v)$,
$\zeta(s)$, $\zeta(s,v)$, some polygamma functions and the Lerch transcendent.
