acb_dirichlet_stieltjes

The function uses an integral representation that allows fast computation for extremely large $n$ [JB2018]. If $n$ is moderate and the precision high enough, it instead expands the Hurwitz zeta function in a power series at $s = 1$ and reads off the last coefficient. For $a = 1$, the power series of the deflated zeta function $\zeta(1+x) - 1/x = \sum_n (-1)^n \gamma_n x^n / n!$ comes from acb_dirichlet_zeta_jet, which at very high precision uses the approximate functional equation (acb_dirichlet_l_fmpq_afe_jet) in time $\tilde O(p^{1.5})$. To compute all of $\gamma_0(a), \ldots, \gamma_n(a)$, expanding $\zeta(s, a)$ once (for instance with acb_poly_zeta_series) is faster than calling this function repeatedly, unless $n$ is at least several hundred.