The Hurwitz zeta function is defined for all complex and by analytic continuation for all complex except for the point , at which it has a simple pole. The Laurent series in a neighborhood of this unique pole is usually written as
The coefficients are known as the generalized Stieltjes constants. The ordinary Stieltjes constants , appearing in the analogous expansion of the Riemann zeta function , are also known as the generalized Euler constants and include the Euler-Mascheroni constant as a special case.111Its generalized analog includes the digamma function , namely , see e.g. [5, Eq. (14)].
This work presents an original method to compute rigorously to arbitrary precision, with the property of remaining fast for arbitrarily large . Such an algorithm has never been published (even in the case of ), despite an extensive literature dedicated to the Stieltjes constants. At the heart of the method is Theorem 1, given below in Section 2, which provides computationally viable integral representations for and . In particular, for and ,
which extends representation (5) from  to the generalized Stieltjes constants. The above expression is similar to the Hermite formula
see e.g. [5, Eq. 13], but more convenient to use for computations since the integrand with the hyperbolic kernel does not possess a removable singularity at Additionally, Section 2 provides some other integral representations for and , some of which may also be suitable for computations (see, in particular, Corollary 1 and Remark 1).
Section 3 describes a robust numerical integration strategy for computing . A crucial step is to determine an approximate steepest descent contour that avoids catastrophic oscillation for large . This is combined with validated integration to ensure that the computation is accurate (indeed, yielding proven error bounds). An open source implementation is available in the Arb library . Section 4 contains benchmark results.
The numbers with were first computed to nine decimal places by Jensen in 1887. Many authors have followed up on this work using an array of techniques. Fundamentally, any method to compute or can be adapted to compute or respectively by taking derivatives. For example, as discussed by Gram , Liang and Todd , Jensen’s calculations of were based on the limit representation
which follows from the Euler–Maclaurin summation formula, so that
while Gram expressed in terms of derivatives of the Riemann function which he evaluated using integer zeta values. Liang and Todd proposed computing either via the Euler-Maclaurin summation formula for , or, as an alternative, via the application of Euler’s series transformation to the alternating zeta function. Bohman and Fröberg  later refined the limit formula technique.
Formula (3) is a differentiated form of Hermite’s integral representation for , which can be interpreted as the Abel-Plana summation formula applied to the series for . As discussed by Blagouchine , the formula (3) has been rediscovered several times in various forms (the case should be credited to Jensen and Franel and dates back to the end of the XIXth century). Ainsworth and Howell  rediscovered the case of (3) and were able to compute up to using Gaussian quadrature.
Keiper  proposed an algorithm based on the approximate functional equation for computing the Riemann function, which upon differentiation yields derivatives of as integrals involving Jacobi theta functions. The Stieltjes constants are then recovered by a power series transformation.
used a version of the limit formula combined with Newton-Cotes quadrature to estimate the resulting sums, and computed accurate values ofup to for and up to for various rational . More recently, Johansson  combined the Euler-Maclaurin formula with fast power series arithmetic for computing , proved rigorous error bounds for this method, and performed the most extensive computation of Stieltjes constants to date resulting in 10000-digit values of for all .
Even more recently, Adell and Lekuona 
have used probability densities for binomial processes to obtain new rapidly convergent series forin terms of Bernoulli numbers.
The drawback of the previous methods is that the complexity to compute is at least linear in . In most cases, the complexity is actually at least quadratic in since the formulas tend to have a high degree of cancellation necessitating use of -digit arithmetic. For the same reason, the space complexity is also usually quadratic in , at least in the most efficient forms of the algorithms. For numerical integration of (3), the difficulty for large lies in the oscillation of the integrand which leads to slow convergence and catastrophic cancellation.222The formula (3) has also been used by Johansson for a numerical implementation of Stieltjes constants in the mpmath library , but the algorithm as implemented in mpmath loses accuracy for large (for example, but mpmath 1.0 computes ).
The fast Euler-Maclaurin method  does allow computing simultaneously to a precision of bits in time if , which is quasi-optimal. However, this is not ideal if we only need or a single .333It is of course also interesting to consider the complexity of computing a single value to variable accuracy . For , the complexity is with all known methods (although the fast Euler-Maclaurin method amortizes this to per coefficient when computing values simultaneously). The exception is which can be computed in time by exploiting its role as a hypergeometric connection constant .
This leads to the question of whether we can compute quickly for any ; ideally, in time depending only polynomially on . If we assume that the accuracy goal is fixed, then any asymptotic formula where is an easily computed function should do the job.
Various asymptotic estimates and bounds for the Stieltjes constants have been published, going back at least to Briggs  and Berndt ,444For the more complete history, see [6, Sect. 3.4]. but the explicit computations by Kreminski and others showed that these estimates were far from precise.
A breakthrough came in 1984, when Matsuoka succeeded in obtaining the first–order asymptotics for the Stieltjes constants , [27, p. 3]. Four years later he derived the complete asymptotic expansion
where , and are the sequences of numbers defined by
and and are the functions defined as
The pair , , is the unique solution of the equation
satisfying and , where and are the gamma and digamma functions respectively, see [24, pp. 49–50].555Matsuoka’s Lemma 1 may be written in our form (6) if we notice that equations (2) and (3) [24, p. 49] actually represent one single equation in which real and imaginary parts were written separately with , and then recall that , where is the complex conjugate of The Matsuoka expansion (5) accurately predicts the behavior of , but is very cumbersome to use. In 2011 Knessl and Coffey  presented a simpler asymptotic formula666If we put in Matsuoka’s expansion (5), we retrieve, after some calculations and several approximations, the same result as Knessl and Coffey (7).
in terms of the slowly varying functions
where is the unique solution of
In (7), the “” symbol signifies asymptotic equality as long as the cosine factor is bounded away from zero. The factor captures the overall growth rate of while the cosine factor explains the local oscillations (and semi-regular sign changes).
More recently, Fekih-Ahmed  has given an alternative asymptotic formula with similar accuracy to (7). Paris  has also generalized (7) to and extended the result to an asymptotic series with higher order correction terms, permitting the determination of several digits for moderately large .
The Matsuoka, Knessl-Coffey, Fekih-Ahmed and Paris formulas were obtained using the standard asymptotic technique of applying saddle point analysis to a suitable contour integral. From a computational point of view, these formulas still have three drawbacks. First, being asymptotic in nature, they only provide a fixed level of accuracy for a fixed , so a different method must be used for small and high precision . Second, the terms in Matsuoka’s expansion and the high-order terms in Paris’s expansion are quite complicated to compute. Third, explicit error bounds are not currently available.
A natural approach to construct an algorithm with the desired properties is to take a similar integral representation and perform numerical integration instead of developing an asymptotic expansion symbolically. The integral representations behind the previous asymptotic formulas do not appear to be convenient for this purpose, since they involve nonsmooth functions (periodic Bernoulli polynomials) or require a summation over several integrals. We therefore use integrals with exponentially decreasing kernels, as in the previous computational work by Ainsworth and Howell , but with the addition of saddle point analysis (which is necessary to handle large ) and a rigorous treatment of error bounds.
2. Integral representations
We obtain the following formulas in terms of elementary integrands that are rapidly decaying and analytic on the path of integration. Although restricted to , they permit computation on the whole and domains through application of the recurrence relations
The Hurwitz zeta function and the generalized Stieltjes constants may be represented by the following integrals
respectively. All formulas hold for complex and such that and .777In these formulas “” signifies that either sign can be taken. Throughout this paper when several “” or “” are encountered in the same formula, it signifies that either the upper signs are used everywhere or the lower signs are used everywhere (but not the mix of them).
In order to prove the above formulas, we will use the contour integration method.888Note that since many formulas with the kernels decaying exponentially fast were already obtained in the past by Legendre, Poisson, Binet, Malmsten, Jensen, Hermite, Lindelöf and many others (see e.g. a formula for the digamma function on p. 541 , or  or ), it is possible that formulas similar or equivalent to those we derive in this section might appear in earlier sources of which we are not aware. In particular, after the publication of the second draft version of this work, we learnt that a formula equivalent to our (11) appears in two books by Srivastava and Choi, [28, p. 92, Eq. (23)] and [29, p. 160, Eq. (23)] respectively. In both sources it appears without proof and without references to other sources.
Consider the following line integral taken along a contour consisting of the interval on the real axis and a semicircle of the radius in the upper half-plane, denoted ,
On the contour the last integral may be bounded as follows:
where we denoted , and
for the purpose of brevity. It can be shown that as tends to infinity and remains integer the integral tends to zero as . For this aim, we first remark that
Since and are both real, except for the case when . Hence
except perhaps at . But at the latter point, since is integer,
Therefore remains always bounded for integer (see also Fig. 1),
and when we have
Thus, accounting for the symmetry of about , we deduce that
From the inequality
it follows that
and since is large, exponential terms on both sides may be neglected. Thus at .999Another
way to obtain the same result is to recall that the integral (18) may be evaluated in terms of the modified
Bessel function of the first kind and the modified Struve function . Using the asymptotic expansions of these special functions
we obtain even a more exact result, namely
if . Hence, making , equality (14) becomes
where the latter integral is taken around an infinitely large semicircle in the upper half-plane. The integrand is not a holomorphic function: it has the poles of the second order at , , due to the hyperbolic secant, and a branch point at due to the term in the numerator. If , the branch point lies outside the integration contour and we may use the Cauchy residue theorem:
Equating (21) with the last result yields
Splitting the interval of integration in two parts and and recalling that
the latter expression may also be written as
Setting in our formulas for , we immediately retrieve our (11)–(11). From the principle of analytic continuation it also follows that above integral formulas are valid for all complex and (because of the branch point which should not lie inside the integration contour). Note that at our formulas (11)–(11) reduce to Jensen’s formulas for the function [5, Eqs. (88)].
Now, in order to get the corresponding formulas for the generalized Stieltjes constant we proceed as follows. The function is holomorphic on the entire complex –plane, and hence, may be expanded into a Taylor series. The latter expansion about reads
But also admits integral representations (23) and (26). Expanding them into the Taylor series in a neighborhood of and equating coefficients in produces formulas (13)–(13). As a particular case of these formulas we obtain formula (5) from  when . ∎
For complex and such that and , the Hurwitz zeta function and the generalized Stieltjes constants admit the representations
Let be such that as . Then, by integration by parts one has
provided the convergence of both integrals and the existence of . Putting straightforwardly yields (28). By virtue of
we also obtain (28). We remark that at formulas (28)–(28) reduce to yet another formula of Jensen for the function [5, Eqs. (88)] and its differentiated form. Formula (29) is obtained analogously from integral (13). ∎
For real , our formulas for the Stieltjes constants may be simplified to
Using similar techniques one may obtain many other integral formulas with kernels decreasing exponentially fast, for instance:
where the latter formulas hold for and For the case, one should remove the term from the last formula. The previous formulas for and also give rise to corresponding expressions for and . For example,
where and are the trigamma and tetragamma functions respectively.101010Some other integral representations with the kernels decreasing exponentially fast for and the polygamma functions may also be found in  and . Also, various relationships between and the polygamma functions are given and discussed in ,  and .
It is similarly possible to derive integral representations for the Lerch transcendent , for example
valid for , or