## 1. Introduction

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

(1) |

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.^{1}^{1}1Its 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 ,

(2) |

which extends representation (5) from [5] to the generalized Stieltjes constants. The above expression is similar to the Hermite formula

(3) |

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 [14]. Section 4 contains benchmark results.

### 1.1. Background

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 [12], Liang and Todd [21], Jensen’s calculations of were based on the limit representation

(4) |

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 [8] 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 [5], 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 [2] rediscovered the case of (3) and were able to compute up to using Gaussian quadrature.

Keiper [18] 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.

Kreminski [20]

used a version of the limit formula combined with Newton-Cotes quadrature to estimate the resulting sums, and computed accurate values of

up to for and up to for various rational . More recently, Johansson [13] 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 [1]

have used probability densities for binomial processes to obtain new rapidly convergent series for

in 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.^{2}^{2}2The formula (3)
has also been used by Johansson for a numerical implementation of Stieltjes constants
in the mpmath library [16],
but the algorithm as implemented in mpmath
loses accuracy for large
(for example,
but mpmath 1.0 computes ).

The fast Euler-Maclaurin method [13]
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 .^{3}^{3}3It 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 [9].

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 [10]
and Berndt [3],^{4}^{4}4For 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 [23], [27, p. 3]. Four years later he derived the complete asymptotic expansion

(5) |

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

(6) |

satisfying and ,
where and are the gamma and digamma
functions respectively, see [24, pp. 49–50].^{5}^{5}5Matsuoka’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 [19] presented a simpler asymptotic
formula^{6}^{6}6If we put
in Matsuoka’s expansion (5), we retrieve, after some calculations
and several approximations, the same result as Knessl and Coffey (7).

(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 [11] has given an alternative asymptotic formula with similar accuracy to (7). Paris [25] 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 [2], 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

(8) |

###### Theorem 1.

The Hurwitz zeta function and the generalized Stieltjes constants may be represented by the following integrals

(11) | |||||

and

(13) | |||||

respectively. All formulas hold for complex and such that
and .^{7}^{7}7In 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.^{8}^{8}8Note 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 [5], or [4] or [22]),
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.

###### Proof.

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 ,

(14) |

On the contour the last integral may be bounded as follows:

(15) |

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

(16) |

Thus, accounting for the symmetry of about , we deduce that

(17) | |||||

From the inequality

it follows that

(18) |

and since is large, exponential terms on both sides may be neglected. Thus at .^{9}^{9}9Another
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

(20) |

if . Hence, making , equality (14) becomes

(21) |

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:

(22) | |||

Equating (21) with the last result yields

(23) |

Splitting the interval of integration in two parts and and recalling that

(24) |

the latter expression may also be written as

(26) | |||||

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 [5] when . ∎

###### Corollary 1.

For complex and such that and , the Hurwitz zeta function and the generalized Stieltjes constants admit the representations

(28) | |||||

and

(29) |

respectively.

###### Proof.

Let be such that as . Then, by integration by parts one has

(30) |

provided the convergence of both integrals and the existence of . Putting straightforwardly yields (28). By virtue of

(31) |

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). ∎

###### Remark 1.

For real , our formulas for the Stieltjes constants may be simplified to

(32) |

and to

(33) |

respectively.

###### Remark 2.

Using similar techniques one may obtain many other integral formulas with kernels decreasing exponentially fast, for instance:

(35) | |||||

(37) | |||||

(38) |

(39) |

(40) |

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.^{10}^{10}10Some other integral representations
with the kernels decreasing exponentially fast for and the polygamma functions
may also be found in [4] and [5]. Also, various relationships
between and the polygamma functions are given and discussed in [4],
[5] and [7].

It is similarly possible to derive integral representations for the Lerch transcendent , for example

valid for , or