Fredrik Johansson
email: fredrik.johansson@gmail.com

Iaroslav V Blagouchine
email: iaroslav.blagouchine@univ-tln.fr

COMPUTING STIELTJES CONSTANTS USING COMPLEX INTEGRATION

Keywords: XXX and, in revised form, XXX. 2010 Mathematics Subject Classification. Primary 11M35, 65D20; Secondary 65G20 Stieltjes constants, Hurwitz zeta function, Riemann zeta function, integral representation, complex integration, numerical integration, complexity, arbitrary-precision arithmetic, rigorous error bounds

The generalized Stieltjes constants γn(v) are, up to a simple scaling factor, the Laurent series coefficients of the Hurwitz zeta function ζ(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 γn(1) to 1000 digits in a minute for any n up to n = 10 100 . We also provide other interesting integral representations for γn(v), ζ(s), ζ(s, v), some polygamma functions and the Lerch transcendent.

Introduction

The Hurwitz zeta function ζ(s, v) = ∞ k=0 (k + v) -s is defined for all complex v = 0, -1, -2, . . . and by analytic continuation for all complex s except for the point s = 1, at which it has a simple pole. The Laurent series in a neighborhood of this unique pole is usually written as [START_REF] Adell | Fast computation of the Stieltjes constants[END_REF] ζ(s, v)

= 1 s -1 + ∞ n=0 (-1) n n! γ n (v)(s -1) n , s ∈ C \ {1}.
The coefficients γ n (v) are known as the generalized Stieltjes constants. The ordinary Stieltjes constants γ n = γ n (1), appearing in the analogous expansion of the Riemann zeta function ζ(s) = ζ(s, 1), are also known as the generalized Euler constants and include the Euler-Mascheroni constant γ 0 = γ = 0.5772156649 . . . as a special case. 1 This work presents an original method to compute γ n (v) rigorously to arbitrary precision, with the property of remaining fast for arbitrarily large n. Such an algorithm has never been published (even in the case of v = 1), 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 ζ(s, v) and γ n (v). In particular, for n ∈ N 0 and Re(v) > 1 2 ,

(2)

γ n (v) = - π 2(n + 1) ˆ+∞ -∞ log n+1 (v -1 2 + ix) cosh 2 πx dx,
which extends representation [START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF] from [START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF] to the generalized Stieltjes constants. The above expression is similar to the Hermite formula

(3) γ n (v) = 1 2v - log v n + 1 log n v -i ˆ∞ 0 dx e 2πx -1 log n (v-ix) v -ix - log n (v+ix) v + ix ,
see e.g. [START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF]Eq. 13], but more convenient to use for computations since the integrand with the hyperbolic kernel sech 2 πx does not possess a removable singularity at x = 0. Additionally, Section 2 provides some other integral representations for ζ(s, v) and γ n (v), 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 γ n (v). A crucial step is to determine an approximate steepest descent contour that avoids catastrophic oscillation for large n. 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 [START_REF] Johansson | Arb: efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]. Section 4 contains benchmark results.

1.1. Background. The numbers γ n with n ≤ 8 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 ζ(s) or ζ(s, v) can be adapted to compute γ n or γ n (v) respectively by taking derivatives. For example, as discussed by Gram [START_REF] Gram | Note sur le calcul de la fonction ζ(s) de Riemann Oversigt[END_REF], Liang and Todd [START_REF] Liang | The Stieltjes constants[END_REF], Jensen's calculations of γ n were based on the limit representation [START_REF] Ia | Rediscovery of Malmsten's integrals, their evaluation by contour integration methods and some related results[END_REF] ζ(s) = lim

N →∞ N k=1 1 k s - N 1-s 1 -s , Re(s) > 0 ,
which follows from the Euler-Maclaurin summation formula, so that

γ n = lim N →∞ N k=1 log n k k - log n+1 N n + 1 , n ∈ N 0 ,
while Gram expressed γ n in terms of derivatives of the Riemann ξ function which he evaluated using integer zeta values. Liang and Todd proposed computing γ n either via the Euler-Maclaurin summation formula for ζ(s), or, as an alternative, via the application of Euler's series transformation to the alternating zeta function.

Bohman and Fröberg [START_REF] Bohman | The Stieltjes function -definition and properties[END_REF] later refined the limit formula technique. Formula (3) is a differentiated form of Hermite's integral representation for ζ(s, v), which can be interpreted as the Abel-Plana summation formula applied to the series for ζ(s, v). As discussed by Blagouchine [START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF], the formula (3) has been rediscovered several times in various forms (the v = 1 case should be credited to Jensen and Franel and dates back to the end of the XIXth century). Ainsworth and Howell [START_REF] Ainsworth | An integral representation of the generalized Euler-Mascheroni constants[END_REF] rediscovered the v = 1 case of (3) and were able to compute γ n up to n = 2000 using Gaussian quadrature. 1 Its generalized analog γ 0 (v) includes the digamma function Ψ(v), namely γ 0 (v) = -Ψ(v), see e.g. [5, Eq. ( 14)].

Keiper [START_REF] Keiper | Power series expansions of Riemann's ξ function[END_REF] proposed an algorithm based on the approximate functional equation for computing the Riemann ξ function, which upon differentiation yields derivatives of ξ(s) as integrals involving Jacobi theta functions. The Stieltjes constants are then recovered by a power series transformation.

Kreminski [START_REF] Kreminski | Newton-Cotes integration for approximating Stieltjes (generalized Euler) constants[END_REF] used a version of the limit formula combined with Newton-Cotes quadrature to estimate the resulting sums, and computed accurate values of γ n (v) up to n = 3200 for v = 1 and up to n = 600 for various rational v. More recently, Johansson [START_REF] Johansson | Rigorous high-precision computation of the Hurwitz zeta function and its derivatives[END_REF] combined the Euler-Maclaurin formula with fast power series arithmetic for computing γ n (v), proved rigorous error bounds for this method, and performed the most extensive computation of Stieltjes constants to date resulting in 10000-digit values of γ n for all n ≤ 10 5 .

Even more recently, Adell and Lekuona [START_REF] Adell | Fast computation of the Stieltjes constants[END_REF] have used probability densities for binomial processes to obtain new rapidly convergent series for γ n in terms of Bernoulli numbers.

The drawback of the previous methods is that the complexity to compute γ n is at least linear in n. In most cases, the complexity is actually at least quadratic in n since the formulas tend to have a high degree of cancellation necessitating use of Ω(n)-digit arithmetic. For the same reason, the space complexity is also usually quadratic in n, at least in the most efficient forms of the algorithms. For numerical integration of (3), the difficulty for large n lies in the oscillation of the integrand which leads to slow convergence and catastrophic cancellation. 2 The fast Euler-Maclaurin method [START_REF] Johansson | Rigorous high-precision computation of the Hurwitz zeta function and its derivatives[END_REF] does allow computing γ 0 , . . . , γ n simultaneously to a precision of p bits in n 2+o (1) time if p = Θ(n), which is quasi-optimal. However, this is not ideal if we only need p = O(1) or a single γ n . 3 This leads to the question of whether we can compute γ n quickly for any n; ideally, in time depending only polynomially on log(n). If we assume that the accuracy goal p is fixed, then any asymptotic formula γ n ∼ G(n) where G(n) 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 [START_REF] Briggs | Some constants associated with the Riemann zeta-function[END_REF] and Berndt [START_REF] Berndt | On the Hurwitz zeta-function[END_REF], 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 firstorder asymptotics for the Stieltjes constants [START_REF] Matsuoka | Generalized Euler constants associated with the Riemann zeta function[END_REF], [27, p. 3]. Four years later he derived the complete asymptotic expansion (5)

γ n ∼ n! e g(b) π N k=0 |h 2k | 2 k+ 1 2 Γ(k + 1 2 ) g ′′ (b) 2 + f ′′ (b) 2 1 2 k+ 1 4 × × cos f (b) -(k + 1 2 ) arctan f ′′ (b) g ′′ (b) + arctan v 2k u 2k
, N = 0, 1, 2, . . . 2 The formula (3) has also been used by Johansson for a numerical implementation of Stieltjes constants in the mpmath library [START_REF] Johansson | mpmath: a Python library for arbitrary-precision floating-point arithmetic[END_REF], but the algorithm as implemented in mpmath loses accuracy for large n (for example, γ 10 4 ≈ -2.21 • 10 6883 but mpmath 1.0 computes -1.25 • 10 6800 ). 3 It is of course also interesting to consider the complexity of computing a single γn value to variable accuracy p. For n ≥ 1, the complexity is Ω(p 2 ) with all known methods (although the fast Euler-Maclaurin method amortizes this to p 1+o (1) per coefficient when computing n = Θ(p) values simultaneously). The exception is γ 0 which can be computed in time p 1+o (1) by exploiting its role as a hypergeometric connection constant [START_REF] Brent | Some new algorithms for high-precision computation of Euler's constant[END_REF]. 4 For the more complete history, see [START_REF] Ia | Expansions of generalized Euler's constants into the series of polynomials in π -2 and into the formal enveloping series with rational coefficients only[END_REF]Sect. 3.4].

where h k , v k and u k are the sequences of numbers defined by

∞ k=0 h k (y -b) k = exp φ(a + iy) -φ(a + ib) + 1 2 φ ′′ (a + ib)(y -b) 2 , u k ≡ Re(h k ) , v k ≡ Im(h k ) ,
and f (b) and g(b) are the functions defined as

g(y) ≡ Re φ(a + iy) , f (y) ≡ Im φ(a + iy) , φ(z) = -(n + 1) log z -z log 2πi + log Γ(z) .
The pair a = Re z, b = Im z, is the unique solution of the equation ( 6)

dφ(z) dz = - n + 1 z -log 2πi + Ψ(z) = 0 ,
satisfying 0 < Im z < Re z and √ n < Re z < n, where Γ(z) and Ψ(z) are the gamma and digamma functions respectively, see [24, pp. 49-50]. 5 The Matsuoka expansion (5) accurately predicts the behavior of γ n , but is very cumbersome to use. In 2011 Knessl and Coffey [START_REF] Knessl | An effective asymptotic formula for the Stieltjes constants[END_REF] presented a simpler asymptotic formula6 (7)

γ n ∼ B √ n e nA cos(an + b)
in terms of the slowly varying functions

A = 1 2 log(α 2 + β 2 ) - α α 2 + β 2 , B = 2 2π(α 2 + β 2 ) 4 (α + 1) 2 + β 2 , a = arctan β α + β α 2 + β 2 , b = arctan β α - 1 2 arctan β α + 1 ,
where β is the unique solution of 2π exp(β tan β) = n cos β β , with 0 < β < 1 2 π, α = β tan β .

In [START_REF] Ia | Three notes on Ser's and Hasse's representations for the zeta-functions[END_REF], the "∼" symbol signifies asymptotic equality as long as the cosine factor is bounded away from zero. The factor Be nA / √ n captures the overall growth rate of γ n while the cosine factor explains the local oscillations (and semi-regular sign changes).

More recently, Fekih-Ahmed [START_REF] Fekih-Ahmed | A new effective asymptotic formula for the Stieltjes constants[END_REF] has given an alternative asymptotic formula with similar accuracy to [START_REF] Ia | Three notes on Ser's and Hasse's representations for the zeta-functions[END_REF]. Paris [START_REF] Paris | An asymptotic expansion for the Stieltjes constants[END_REF] has also generalized [START_REF] Ia | Three notes on Ser's and Hasse's representations for the zeta-functions[END_REF] to γ n (v) and extended the result to an asymptotic series with higher order correction terms, permitting the determination of several digits for moderately large n.

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 n, so a different method must be used for small n and high precision p. 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 [START_REF] Ainsworth | An integral representation of the generalized Euler-Mascheroni constants[END_REF], but with the addition of saddle point analysis (which is necessary to handle large n) and a rigorous treatment of error bounds.

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 Re(v) > 1 2 , they permit computation on the whole (s, v) and (n, v) domains through application of the recurrence relations

(8) ζ(s, v) = ζ(s, v + 1) + 1 v s , γ n (v) = γ n (v + 1) + log n v v .
Theorem 1. The Hurwitz zeta function ζ(s, v) and the generalized Stieltjes constants γ n (v) may be represented by the following integrals

ζ(s, v) = π 2(s -1) ˆ+∞ -∞ (v -1 2 ± ix) 1-s cosh 2 πx dx (9) = π 2(s -1) ˆ∞ 0 (v -1 2 -ix) 1-s + (v -1 2 + ix) 1-s cosh 2 πx dx (10) = π 2(s -1) ˆ∞ 0 cos (s -1) arctan 2x 2v-1 (v 2 -v + 1 4 + x 2 ) 1 2 (s-1) cosh 2 πx dx (11) and γ n (v) = - π 2(n + 1) ˆ+∞ -∞ log n+1 (v -1 2 ± ix) cosh 2 πx dx (12) = - π 2(n + 1) ˆ∞ 0 log n+1 (v -1 2 -ix) + log n+1 (v -1 2 + ix) cosh 2 πx dx (13)
respectively. All formulas hold for complex v and s such that Re(v) > 1 2 and s = 1. 7 In order to prove the above formulas, we will use the contour integration method. 8 7 In 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).

8 Note 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 [START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF], or [START_REF] Ia | Rediscovery of Malmsten's integrals, their evaluation by contour integration methods and some related results[END_REF] or [START_REF] Lindelöf | Le calcul des résidus et ses applications à la théorie des fonctions[END_REF]), 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 [START_REF] Fekih-Ahmed | A new effective asymptotic formula for the Stieltjes constants[END_REF] appears in two books by Srivastava Proof. Consider the following line integral taken along a contour C consisting of the interval [-R, +R] on the real axis and a semicircle of the radius R in the upper half-plane, denoted C R , ( 14)

‰ C (a -iz) 1-s cosh 2 πz dz = ˆ+R -R (a -ix) 1-s cosh 2 πx dx + ˆCR (a -iz) 1-s cosh 2 πz dz .
On the contour C R the last integral may be bounded as follows:

ˆCR (a -iz) 1-s cosh 2 πz dz = R ˆπ 0 (a -iRe iϕ ) 1-s e iϕ cosh 2 (πRe iϕ ) dϕ ≤ ≤ R max ϕ∈[0,π] (a -iRe iϕ ) 1-s • I R ≤ R max ϕ∈[0,π] |a| 2 +2R (a x sin ϕ -a y cos ϕ) + R 2 1 2 Re (1-s) e π|Im(1-s)| I R (15)
where we denoted a x ≡ Re(a), a y ≡ Im(a) and

I R ≡ ˆπ 0 dϕ | cosh(πRe iϕ )| 2 , R > 0 ,
for the purpose of brevity. It can be shown that as R tends to infinity and remains integer the integral I R tends to zero as 1/R. For this aim, we first remark that 1

| cosh(πRe iϕ )| 2 = 2 cosh(2πR cos ϕ) + cos(2πR sin ϕ) .
Since R and ϕ are both real, cosh(2πR cos ϕ) > 1 except for the case when cos ϕ = 0. Hence cosh(2πR cos ϕ) + cos(2πR sin ϕ) > 0 , except perhaps at ϕ = 1 2 π. But at the latter point, since R is integer, cosh(2πR cos ϕ) + cos(2πR sin ϕ)

ϕ= 1 2 π = 1 + cos(2πR) = 2 .
Therefore | cosh(πRe iϕ )| -2 remains always bounded for integer R (see also Fig. 1), and when R → ∞ we have [START_REF] Johansson | mpmath: a Python library for arbitrary-precision floating-point arithmetic[END_REF] 1

| cosh(πRe iϕ )| 2 =      O(e -2πR cos ϕ ) , ϕ ∈ [0, 1 2 π] , O(e +2πR cos ϕ ) , ϕ ∈ [ 1 2 π, π] .
Thus, accounting for the symmetry of | cosh(πRe iϕ )| -2 about ϕ = 1 2 π, we deduce that

I R = ˆπ 0 2 dϕ cosh(2πR cos ϕ) + cos(2πR sin ϕ) = ˆπ 2 0 4 dϕ cosh(2πR cos ϕ) + cos(2πR sin ϕ)
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. clearly displays the boundness of the latter (R integer). Note also that at large R the contribution of the point ϕ = 1 2 π to the integral IR becomes infinitely small (its height equals one, while the width tends to zero).

= O ˆπ 2 0 e -2πR cos ϕ dϕ = O ˆπ 2 0 e -2πR sin ϑ dϑ , R → ∞ . (17) From the inequality 2ϑ π ≤ sin ϑ ≤ ϑ , ϑ ∈ [0, 1 2 π] , it follows that (18) 1 -e -π 2 R 2πR ≤ ˆπ 2 0 e -2πR sin ϑ dϑ ≤ 1 -e -2πR 4R ,
and since R is large, exponential terms on both sides may be neglected. Thus 9 Inserting this result into (15), we obtain

I R = O(1/R) at R → ∞.
ˆCR (a -iz) 1-s cosh 2 πz dz → 0 as R → ∞ , R ∈ N , ( 20
)
9 Another way to obtain the same result is to recall that the integral (18) may be evaluated in terms of the modified Bessel function In(z) of the first kind and the modified Struve function Ln(z). Using the asymptotic expansions of these special functions we obtain even a more exact result, namely [START_REF] Knessl | An effective asymptotic formula for the Stieltjes constants[END_REF] ˆπ

2 0 e -2πR sin ϑ dϑ = π 2 I 0 (2πR) -L 0 (2πR) ∼ 1 2πR , R → ∞ ,
i.e. the integral asymptotically tends to the left bound [START_REF] Keiper | Power series expansions of Riemann's ξ function[END_REF].

if Re(s) > 1. Hence, making R → ∞, equality ( 14) becomes ( 21)

ˆ+∞ -∞ (a -ix) 1-s cosh 2 πx dx = ‰ C (a -iz) 1-s cosh 2 πz dz
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 z = z n ≡ i n + 1 2 , n ∈ N 0 , due to the hyperbolic secant, and a branch point at z = -ia due to the term in the numerator. If Re(a) > 0, the branch point lies outside the integration contour and we may use the Cauchy residue theorem:

‰ C (a -iz) 1-s cosh 2 πz dz = 2πi ∞ n=0 res z=zn (a -iz) 1-s cosh 2 πz = (22) = - 2i π ∞ n=0 ∂ ∂z (a -iz) 1-s z=i(n+ 1 2 ) = = 2(s -1) π ∞ n=0 (a + 1 2 + n) -s = 2(s -1) π ζ(s, a + 1 2 ) .
Equating [START_REF] Liang | The Stieltjes constants[END_REF] with the last result yields

(23) ζ(s, a + 1 2 ) = π 2(s -1) ˆ+∞ -∞ (a -ix) 1-s cosh 2 πx dx , Re(a) > 0 .
Splitting the interval of integration in two parts (-∞, 0] and [0, +∞] and recalling that (24) (a + ix) s + (a -ix) s = 2(a 2 + x 2 ) s 2 cos s arctan x a the latter expression may also be written as

ζ(s, a + 1 2 ) = π 2(s -1) ˆ∞ 0 (a + ix) 1-s + (a -ix) 1-s cosh 2 πx dx (25) = π s -1 ˆ∞ 0 cos (s -1) arctan x a (a 2 + x 2 ) 1 2 (s-1) cosh 2 πx dx , Re(a) > 0 . ( 26
)
Setting a = v -1 2 in our formulas for ζ(s, a + 1 2 ), we immediately retrieve our ( 9)- [START_REF] Fekih-Ahmed | A new effective asymptotic formula for the Stieltjes constants[END_REF]. From the principle of analytic continuation it also follows that above integral formulas are valid for all complex s = 1 and Re(a) > 0 (because of the branch point which should not lie inside the integration contour). Note that at v = 1 our formulas ( 9 Now, in order to get the corresponding formulas for the generalized Stieltjes constant γ n (v) we proceed as follows. The function (s-1)ζ(s, a+ 1 2 ) is holomorphic on the entire complex s-plane, and hence, may be expanded into a Taylor series. The latter expansion about s = 1 reads

(s -1)ζ(s, a + 1 2 ) = 1 + ∞ n=0 (-1) n γ n (a + 1 2 ) n! (s -1) n+1 , s ∈ C \{1}. But (s-1)ζ(s, a+ 1
2 ) also admits integral representations [START_REF] Matsuoka | Generalized Euler constants associated with the Riemann zeta function[END_REF] and [START_REF] Paris | An asymptotic expansion for the Stieltjes constants[END_REF]. Expanding them into the Taylor series in a neighborhood of s = 1 and equating coefficients in (s -1) n+1 produces formulas ( 12)-( 13). As a particular case of these formulas we obtain formula (5) from [START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF] when v = 1.

Corollary 1. For complex v and s such that Re(v) > 1 2 and s = 1, the Hurwitz zeta function ζ(s, v) and the generalized Stieltjes constants γ n (v) admit the representations

ζ(s, v) = (v -1 2 ) 1-s s -1 + i ˆ∞ 0 (v -1 2 -ix) -s -(v -1 2 + ix) -s e 2πx + 1 dx (27) = (v -1 2 ) 1-s s -1 -2 ˆ∞ 0 sin s arctan 2x 2v-1 (v 2 -v + 1 4 + x 2 )
s 2

• dx e 2πx + 1 [START_REF] Srivastava | Series Associated with the Zeta and Related Functions[END_REF] and ( 29)

γ n (v) = - log n+1 (v -1 2 ) n + 1 + + i ˆ∞ 0 log n (v -1 2 -ix) v -1 2 -ix - log n (v -1 2 + ix) v -1 2 + ix dx e 2πx + 1 respectively.
Proof. Let f (x) be such that f (x) = o(e 2πx ) as x → ∞. Then, by integration by parts one has

(30) ˆ∞ 0 f (x) cosh 2 πx dx = f (0) π + 2 π ˆ∞ 0 f ′ (x) e 2πx + 1 dx ,
provided the convergence of both integrals and the existence of f (0). Putting f (x) = (a + ix) 1-s + (a -ix) 1-s straightforwardly yields [START_REF] Saad-Eddin | On two problems concerning the Laurent-Stieltjes coefficients of Dirichlet L-series[END_REF]. By virtue of (31

) (a + ix) s -(a -ix) s = 2i(a 2 + x 2 ) s 2 sin s arctan x a ,
we also obtain [START_REF] Srivastava | Series Associated with the Zeta and Related Functions[END_REF]. We remark that at v = 1 formulas ( 27)- [START_REF] Srivastava | Series Associated with the Zeta and Related Functions[END_REF] reduce to yet another formula of Jensen for the ζ function [5, Eqs. (88)] and its differentiated form. Formula ( 29) is obtained analogously from integral [START_REF] Johansson | Rigorous high-precision computation of the Hurwitz zeta function and its derivatives[END_REF].

Remark 1. For real v > 1 2 , our formulas for the Stieltjes constants may be simplified to

(32) γ n (v) = - π n + 1 Re ˆ∞ 0 log n+1 (v -1 2 ± ix) cosh 2 πx dx
and to

γ n (v) = - log n+1 (v -1 2 ) n + 1 ± 2 Im ˆ∞ 0 log n (v -1 2 ± ix) v -1 2 ± ix • dx e 2πx + 1 (33) respectively.
Remark 2. Using similar techniques one may obtain many other integral formulas with kernels decreasing exponentially fast, for instance:

ζ(s) = 1 1 -2 s-1 1 2 + 1 2i ˆ+∞ -∞ (1 -ix) -s dx sinh πx = (34) = 1 1 -2 s-1 1 2 + ˆ∞ 0 sin (s arctan x) (1 + x 2 ) s 2 sinh πx dx , (35) ζ(s, v) = ± i π 2 (s -1)(s -2) ˆ+∞ -∞ (v -1 2 ± ix) 2-s sinh πx cosh 3 πx dx = (36) = 2 π 2 (s -1)(s -2) ˆ∞ 0 sin (s -2) arctan 2x 2v-1 (v 2 -v + 1 4 + x 2 ) s 2 -1 • sinh πx cosh 3 πx dx , ( 37
) (38) ζ(s, v) = 3π 3 (s -1)(s -2)(s -3) ˆ+∞ -∞ 1 cosh 4 πx - 2 3 cosh 2 πx dx (v -1 2 ± ix) s-3 , ( 39
) γ 1 = π 2 24 - γ 2 2 - log 2 2 2 + log 2 π 2 -log 2 • log π + ˆ∞ 0 arctan x • log(1 + x 2 ) sinh πx dx , γ n (v) = (-1) n+1 4π n(n -1)ζ (n-2) (3, v) + 3nζ (n-1) (3, v) + 2ζ (n-2) (3, v) - - 3π 4(n + 1) ˆ+∞ -∞ log n+1 (v -1 2 ± ix) cosh 4 πx dx , (40)
where the latter formulas hold for Re(v) > 1 2 and n = 2, 3, 4, . . . For the n = 1 case, one should remove the n(n -1)ζ (n-2) (3, v) term from the last formula. The previous formulas for ζ(s, v) and γ n (v) also give rise to corresponding expressions for Ψ(v) and log Γ(v). For example,

Ψ(v) = - Ψ 2 (v) 4π 2 + 3π 4 ˆ∞ 0 log (v 2 -v + 1 4 + x 2 ) cosh 4 πx dx , log Γ(v) = 1 2 log 2π + v - 1 2 Ψ(v) -1 -π ˆ∞ 0 x arctan 2x 2v-1 cosh 2 πx dx , log Γ(v) = 1 2 log 2π + v - 1 2 Ψ(v) + Ψ 2 (v) 4π 2 -1 - - Ψ 1 (v) 4π 2 - 3π 2 ˆ∞ 0 x arctan 2x 2v-1 cosh 4 πx dx ,
where Ψ 1 (v) and Ψ 2 (v) are the trigamma and tetragamma functions respectively. 10 It is similarly possible to derive integral representations for the Lerch transcen-

dent Φ(z, s, v) = ∞ n=0 z n (n + v) -s , for example Φ(z, s, v) = v -s 2 + i 2 ˆ∞ 0 (-z) ix (v + ix) -s -(-z) -ix (v -ix) -s
sinh πx dx 10 Some other integral representations with the kernels decreasing exponentially fast for log Γ(v)

and the polygamma functions may also be found in [START_REF] Ia | Rediscovery of Malmsten's integrals, their evaluation by contour integration methods and some related results[END_REF] and [START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF]. Also, various relationships between log Γ(z) and the polygamma functions are given and discussed in [START_REF] Ia | Rediscovery of Malmsten's integrals, their evaluation by contour integration methods and some related results[END_REF], [START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF] and [START_REF] Ia | Three notes on Ser's and Hasse's representations for the zeta-functions[END_REF].

= v -s 2 + ˆ∞ 0 cos(x log z) sin s arctan x v -sin(x log z) cos s arctan x v (v 2 + x 2 ) s 2 tanh πx dx, valid for z > 0, or ζ(s, v) = v -s 2 + i 2 ˆ∞ 0 e -πx (v + ix) -s -e +πx (v -ix) -s sinh πx dx = v -s 2 + i 2 ˆ∞ 0 (v + ix) -s -(v -ix) -s tanh πx dx = v -s 2 + ˆ∞ 0 sin s arctan x v (v 2 + x 2 ) s 2 tanh πx dx ,
whose integrands are not of exponential decay, despite the presence of the hyperbolic cosecant. 11 At the same time, the above formula for Φ(z, s, v) is suitable for z > 0, while the same formula for the negative first argument reads

Φ(-z, s, v) = v -s 2 + i 2 ˆ∞ 0 z ix (v + ix) -s -z -ix (v -ix) -s sinh πx dx = v -s 2 + ˆ∞ 0 cos(x log z) sin s arctan x v -sin(x log z) cos s arctan x v (v 2 + x 2 ) s 2 sinh πx dx,
z > 0, and the integrand decreases exponentially fast.

Computation of γ n (v) by integration

For the computation of γ n (v), we use formulas ( 12)-( 13), (32). 12 For the purpose of brevity throughout this section, we write a for v -1 2 . We denote the integrand (with a and n as implicit parameters) and the half-line integral by

(41) f (z) ≡ log n+1 (a + iz) cosh 2 πz , I n (a) ≡ ˆ∞ 0 f (x) dx
respectively. After applying (8) as needed to ensure Re(a) > 0 (or better, Re(a) ≥ 1 2 to stay some distance away from the logarithmic branch point and avoid convergence issues during the numerical integration to follow), we may compute

(42) γ n (v) = - π (n + 1) •    2 Re(I n (a)), Im(a) = 0 I n (a) + I n (a), Im(a) = 0.
where " " stands for the complex conjugate. For a given accuracy goal of p bits, we aim to compute I n (a) with a relative error less than 2 -p . More precisely, we assume use of ball arithmetic [START_REF] Van Der Hoeven | Ball arithmetic[END_REF], and we aim to compute an enclosure with relative radius less than 2 -p . A first important observation is that the computations must be done with a working precision of about p + log 2 n bits for p-bit accuracy, due to the sensitivity of the integrand. In other words, we lose about log 2 n bits to the exponents of the floating-point numbers when evaluating exponentials. Heuristically, a few more guard bits in addition to this will be sufficient to account for all rounding errors, and the computed ball provides a final certificate.

A technical point is that we cannot make any a priori statements about the relative error of γ n (v) since we do not have lower bounds for |γ n (v)|. Cancellation is possible in the final addition (or extraction of the real part) in (42). This should roughly correspond to multiplying by the cosine factor in [START_REF] Ia | Three notes on Ser's and Hasse's representations for the zeta-functions[END_REF]; it is reasonable to set the accuracy goal with respect to the nonoscillatory factor Be nA / √ n. We primarily have in mind "small" parameters v (for example v = 1) such that |v|≪ n if n is large. The algorithm works for any complex v where γ n (v) is defined, but we do not specifically address optimization for large |v| which therefore may result in deteriorating efficiency and less precise output enclosures.

3.1. Estimation of the tail. We approximate I n (a), given by (41), by the truncated integral ´N 0 f (x) dx for some N > 0. The following theorem provides an upper bound for the tail T N .

Theorem 2. Let (43)

T N ≡ ˆ∞ N log n+1 (a + ix) cosh 2 πx dx and assume Re(a) > 0. Then, the following bound holds:

(44) |T N | < 0.934 e -2πN | log(a + N i)| n+1 , N ≥ n + 2 + |Im(a)| .
Proof. For x ≥ 0, using |log ′ (z)|= 1/|z| and the assumptions on a and N gives

| log(a + i(N + x))| n+1 = | log(a + N i)| n+1 1 + log(a + i(N + x)) -log(a + N i) log(a + N i) n+1 ≤ | log(a + N i)| n+1 1 + x |a + N i| | log(a + N i)| n+1 ≤ |log(a + N i)| n+1 exp (n + 1)x |a + N i| log|a + N i| ≤ |log(a + N i)| n+1 exp(2x).
Since sech 2 (x) < 4e -2x on the whole real line, we have

|T N | ≤ ˆ∞ N 4e -2πx |log (a + ix)| n+1 dx ≤ 4e -2πN |log(a + N i)| n+1 ˆ∞ 0 exp (-2πx + 2x) dx
and the last integral equals 1 2 (π -1) -1 .

We can select N by starting with N = n+ 2 + |Im(a)| and repeatedly doubling N until |T N |≤ 2 -p-20 , say. This bound does not need to be tight since the integration algorithm, described later, discards negligible segments cheaply through bisection.

Remark 3. The bound in the previous theorem can be made slightly sharper, although this does not matter for the algorithm. By using the same line of reasoning as above, the inequality | log(1 + z)| < |z| holding true for Re(z) ≥ 0 and the fact that the error function is always lesser than 1, one can obtain, for example,

|T N | < 0.637 1 + λ √ 2 e λ 2 8π e -2πN | log(a + N i)| n+1 , N ≥ 4(n + 1) 2 λ 2 + |Im(a)| ,
where λ is some positive parameter lesser than 2 (the smaller λ, the better this estimation; for λ < 0.649 this estimation outperforms (44), but N must be large with respect to n 2 ). Moreover, for 0 < λ ≪ 1 we may neglect the term O(λ) between the parenthesis, and hence obtain

|T N | 0.637 e -2πN | log(a + N i)| n+1 , N ≫ n 2 .
Both these bounds and the value 0.637, coming from 2/π, are in good agreement with the numerical results. Note that if (44) is suitable for cases in which N is comparable to n, the above estimations are suitable only for cases of large and extra-large N with respect to n 2 .

3.2.

Cancellation-avoiding contour. For small n, the integral ´N 0 f (x) dx can be computed directly. For large n, the integrand oscillates on the real line and a higher working precision must be used due to cancellation. At least for v = 1, the amount of cancellation can be calculated accurately by numerically computing the maximum value of |f (x)| on 0 ≤ x < ∞ and comparing this magnitude to the asymptotic formula [START_REF] Ia | Three notes on Ser's and Hasse's representations for the zeta-functions[END_REF]. For example, we need about 30 extra bits when n = 10 3 , 1740 bits when n = 10 6 , and 4 • 10 5 bits when n = 10 9 .

For n larger than about 10 3 , we shift the path to eliminate the cancellation problem. The integrand can be written as Assuming that n ≫ |a|, the function exp(g(z)) has a single saddle point in the right half-plane. The saddle point equation g ′ (ω) = 0 can be reduced to (47) (n + 1) + 2πi (a + iω) log (a + iω) = 0 which admits the closed-form solution

(48) ω = i a - u W 0 (u) , u = (n + 1)i 2π .
where W 0 (u) is the principal branch of the Lambert W function. Only the principal branch works, a fact which is not obvious from the symbolic form of the solution but which can be checked numerically.

We can now integrate along four segments

ˆN 0 f (x) dx = ˆM 0 f (z) dz + ˆM+Ci M f (z) dz + ˆN+Ci M+Ci f (z) dz + ˆN N +Ci f (z) dz
with the choice of vertical offset C = Im(ω) to (approximately) minimize the peak magnitude of f (t + Ci) on M ≤ t ≤ N . The left point M > 0 just serves to avoid the poles of the integrand on the imaginary axis and the nearby vertical branch cut of the logarithm; we can for instance take M = 10. Numerical tests (compare Fig. 2) confirm that there is virtually no cancellation with this contour (again, assuming that |v| is not too large). The path does not exactly pass through the saddle point of f (z), but since h(z) is exponentially close to a constant, the perturbation is negligible. The deviation between the straightline path through the saddle point and the actual steepest descent contour also has negligible impact on the numerical stability.

We note that the complex Lambert W function can be computed with rigorous error bounds [START_REF] Johansson | Computing the Lambert W function in arbitrary-precision complex interval arithmetic[END_REF]. However, it is not actually necessary to compute ω rigorously for this application since the integration follows a connected path and ball arithmetic will account for the actual cancellation; it is sufficient to use a floating-point approximation for ω with heuristic accuracy of about log 2 n bits. For example, an approximation of ω computed with 53-bit machine arithmetic is sufficient up to about n = 10 15 .

3.3.

Integration and bounds near the saddle point. The main task of integrating f (z) along one or four segments in the plane is not difficult in principle, since f (z) is analytic (and non-oscillatory) in a neighborhood of each segment. Constructing a reliable and fast algorithm, in particular for extremely large n, does nevertheless require some attention to detail.

Gauss-Legendre quadrature is a good option, and was already used by Ainsworth and Howell [START_REF] Ainsworth | An integral representation of the generalized Euler-Mascheroni constants[END_REF], who, however, did not prove any error bounds since "The integrand is much too complex to use the standard remainder terms". To obtain rigorous error bounds and ensure rapid convergence with a manageable level of manual error analysis, we use the self-validating Petras algorithm [START_REF] Petras | Self-validating integration and approximation of piecewise analytic functions[END_REF] which was recently adapted for arbitrary-precision ball arithmetic and implemented in the Arb library [START_REF] Johansson | Numerical integration in arbitrary-precision ball arithmetic[END_REF].

The Petras algorithm combines Gauss-Legendre quadrature with adaptive bisection. Given a segment [α, β], the algorithm first evaluates the direct enclosure (β -α)f ([α, β]) and uses this if the error is negligible (which in this application always occurs near the tail ends of the integral when p ≪ n). Otherwise, it bounds the error of d-point quadrature to ensure a relative error smaller than 2 -p , the segment [α, β] is bisected and the same procedure is applied recursively. The remaining issue is the evaluation of the integrand. The pointwise evaluations w k f (z k ) pose no problem: here we simply use (41) directly. It is slightly more complicated to compute good enclosures for f (z) on wide intervals representing z, which is needed both for the direct enclosures on subintervals f ([α, β]) and for the bounds on ellipses. 13 Bounding the integrand on wide ellipses (or enclosing rectangles) by evaluating (41) or ( 45)-( 46) directly in interval or ball arithmetic results at best in n 1/2+o(1) complexity as n → ∞. 14 The explanation for this phenomenon is that f (z) is a quotient of two functions f 1 (z) = log n+1 (a + iz), f 2 (z) = cosh 2 πz that individually vary rapidly near the saddle point, i.e.

ˆβ α f (z)dz ≈ d k=1 w k f (z k ) in terms of the magnitude |f (z)| on a Bernstein ellipse E around [α, β]: if f (z) is analytic on E with max z∈E |f (z)|≤ V ,
(49)

f 1 (z + ε) f 1 (z) ∼ f 2 (z + ε) f 2 (z) ∼ e 2πε
while f 1 /f 2 is nearly constant. Direct evaluation fails to account for this correlation, which is an example of the dependency problem in interval arithmetic. Therefore, although f (z) is nearly constant close to the saddle point, direct upper bounds for |f (z)| are exponentially sensitive to the width of input intervals, and this forces the integration algorithm to bisect down to subsegments of width O(1) around the saddle point. Since the Gaussian peak of the integrand around the saddle point has an effective width of O(n 1/2 ), the integration algorithm has to bisect down to O(n 1/2 ) subsegments before converging.

To solve this problem, we compute tighter bounds on wide intervals using the standard trick of Taylor expanding with respect to a symbolic perturbation ε. Taylor's theorem applied to exp(g(z)) gives

(51) exp (g(m + ε)) = exp (g(m)) exp g ′ (m)ε + ˆε 0 g ′′ (m + t)(ε -t)dt
for all |ε|≤ r.

To implement the bound (50), we compute g(m) and g ′ (m) in ball arithmetic (where m is an exact floating-point number), using the formula

g ′ (m) = i(n + 1) (a + im) log(a + im) -2π.
The behavior near the saddle point is now captured precisely by the cancellation in g ′ (m). At least log 2 n bits of precision must be used to evaluate exp(g(m)) (to ensure that the magnitude of the integrand near the peak is approximated accurately) 13 The complex ball arithmetic in Arb actually uses rectangles with midpoint-radius real and imaginary parts rather than complex disks, so ellipses will always be represented by enclosing rectangles (with up to a factor √ 2 overestimation), but this detail is immaterial to the principle of the algorithm. 14 In fact, the complexity becomes n 1+o (1) when using ball arithmetic with a fixed precision for the radii (30 bits in Arb). The n 1/2+o (1) estimate holds when the endpoints are tracked accurately.

and also to evaluate g ′ (m) (to ensure that the remainder after the catastrophic cancellation is evaluated accurately). Finally, to compute G, we evaluate g ′′ (z) =

(n + 1) 1 + 1 log t t 2 log t , t = a + iz directly over the complex ball representing z. As a minor optimization, we can compute lower bounds for |t| and |log t|. This completes the algorithm.

3.4. Asymptotic complexity. If the accuracy goal p is fixed (or grows sufficiently slowly compared to n), then we can argue heuristically that the bit complexity of computing γ n to p-bit accuracy with this algorithm is log 2+o (1) n. This estimate accounts for the bisection depth around the saddle point as well as the extra precision of log 2 n bits. The logarithmic complexity agrees well with the actual timings (presented in the next section).

We stop short of attempting to prove a formal complexity result, which would require more detailed calculations and careful accounting for the accuracy of the enclosures in ball arithmetic as well as details about the integration algorithm. We have delegated as much work as possible to a general-purpose integration algorithm in order to minimize the analysis necessary for a complete implementation. However, in future work, it would be interesting to pursue such analysis not just for this specific problem, but more generally for evaluating classes of parametric integrals using the combination of saddle point analysis and numerical integration.

If we on the other hand fix n and consider varying p, then the asymptotic bit complexity is of course p 2+o (1) since Gaussian quadrature uses O(p) evaluations of the integrand on a fixed segment and O(log p) segments are sufficient.

Implementation and benchmark results

The new integration algorithm has been implemented in Arb [START_REF] Johansson | Arb: efficient arbitrary-precision midpoint-radius interval arithmetic[END_REF]. 15 The method acb dirichlet stieltjes computes γ n (v), given a complex ball representing v, an arbitrary-size integer n, and a precision p. The working precision is set automatically so that the result will be accurate to about p bits, at least when v = 1. This method selects automatically between two internal methods:

• acb dirichlet stieltjes integral uses the new integration algorithm.

• acb dirichlet stieltjes em is a wrapper around the existing code for computing the Hurwitz zeta function using Euler-Maclaurin summation [START_REF] Johansson | Rigorous high-precision computation of the Hurwitz zeta function and its derivatives[END_REF]. For very small n, the integration algorithm is one-three orders of magnitude slower than Euler-Maclaurin summation, but the cost of the latter increases rapidly with n. Integration was found to be faster when n > max(100, p/2), and this automatic cutoff is used in the code. We remark that the Euler-Maclaurin code actually computes γ 0 (v), . . . , γ n (v) simultaneously and reads off the last entry. At this time, we do not have an implementation of the Euler-Maclaurin formula optimized for a single γ n (v) value, which would be significantly faster for n from about 10 to 10 3 .

Table 1 shows the time in seconds to evaluate the ordinary Stieltjes constants γ n to a target accuracy p of 64 bits (about 18 digits), 333 bits (about 100 digits) and 3333 bits (just more than 1000 digits) on an Intel Core i5-4300U CPU running 64-bit Ubuntu Linux. Here we only show the timing results for the Arb method acb dirichlet stieltjes integral, omitting use of Euler-Maclaurin summation. As expected, the running time of our algorithm only depends weakly on n. The performance is also reasonable for large p. The timings fluctuate slightly rather than increasing monotonically with n, which appears to be an artifact of the local adaptivity of the integration algorithm.

Mathematica returns incorrect answers for large n when using machine precision. At higher precision, the performance is consistent up to about n = 10 4 , but the running time then starts to increase rapidly. With n = 10 5 and 100-digit or 1000digit precision, Mathematica did not finish when left to run overnight.

Mathematica uses Keiper's algorithm according to the documentation [START_REF]Some notes on internal implementation[END_REF], but unfortunately we do not have details about the implementation. The timings and failures for large n are seemingly consistent with use of numerical integration in some form without the precautions we have taken against oscillation problems. We also mention that Maple is much slower than Mathematica, taking 0.1 seconds to compute γ 10 , a minute to compute γ 1000 and six minutes to compute γ 2000 to 10 digits. 4.1. Multi-evaluation. Table 2 compares the performance of Euler-Maclaurin summation and the integration method in Arb for computing γ 0 (v), . . . , γ n (v) simultaneously. With the integration algorithm, this means making n + 1 independent evaluations, while the Euler-Maclaurin algorithm only has to be executed once. Despite this, integration still wins for sufficiently large n, unless p also is large. As a sanity check, γ 10 5 agrees with the previous record Euler-Maclaurin computation [START_REF] Johansson | Rigorous high-precision computation of the Hurwitz zeta function and its derivatives[END_REF]. The value of γ n also agrees with the Knessl-Coffey formula [START_REF] Ia | Three notes on Ser's and Hasse's representations for the zeta-functions[END_REF] to about log 10 n digits, in perfect agreement with the error term in this asymptotic approximation being O(1/n).

For γ n (v) with a nonreal v, the computation time roughly doubles since two integrals are computed. With v = 1, we can for instance compute: These values similarly agree to log 10 n digits with the leading-order truncation of Paris's generalization [START_REF] Paris | An asymptotic expansion for the Stieltjes constants[END_REF] of the Knessl-Coffey formula, providing both a check on our implementation and an independent validation of Paris's results.

Discussion

A few possible optimizations of the integration algorithm are worth pointing out. The adaptive integration strategy in Arb can probably be improved, which should give a constant factor speedup. The working precision could also likely be reduced by a preliminary rescaling near the saddle point.

For evaluating a range of γ n (v) simultaneously, one could perform vector-valued integration and recycle the evaluations of log(a + iz) and cosh 2 πz. It would be interesting to compare this approach to simultaneous evaluation with the Euler-Maclaurin formula.

It would also be interesting to investigate use of double exponential quadrature instead of Gaussian quadrature.

The computational part of this study was done for two purposes: first, to develop working code for Stieltjes constants as part of the collection of rigorous special function routines in the Arb library, and second, to test the integration algorithm [START_REF] Petras | Self-validating integration and approximation of piecewise analytic functions[END_REF][START_REF] Johansson | Numerical integration in arbitrary-precision ball arithmetic[END_REF] for a family of integrals involving large parameters. We do not have a concrete application in mind for the code, but the Stieltjes constants are potentially useful in various types of analytic computations involving the Riemann zeta function, and large-n evaluation can be useful for testing the accuracy of asymptotic formulas for Stieltjes constants and related quantities.

The technique of evaluating parametric integrals by integrating numerically along a steepest descent contour is, of course, well established in the literature on computational methods for special functions, but such an algorithm has not previously been published for Stieltjes constants. The use of rigorous integration techniques in such a setting has also been explored very little in earlier work. The most important lesson learned here is that the heavy lifting can be done by the integration algorithm, requiring only an elementary pen-and-paper analysis of the integrand. The same technique should be effective for rigorously computing many other number sequences and special functions given by similar integral representations. On that note, it would be interesting to search for more integral representations similar to those obtained in Section 2. Many such representations with the integrands decreasing exponentially fast for log Γ(z) and for the polygamma functions may be found in [START_REF] Ia | Rediscovery of Malmsten's integrals, their evaluation by contour integration methods and some related results[END_REF] and [START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF].

Figure 1 .

 1 Figure1. 3D-plot of cosh πRe iϕ -2 for R ∈[START_REF] Adell | Fast computation of the Stieltjes constants[END_REF][START_REF] Johansson | mpmath: a Python library for arbitrary-precision floating-point arithmetic[END_REF] and ϕ ∈ [0, π] clearly displays the boundness of the latter (R integer). Note also that at large R the contribution of the point ϕ = 1 2 π to the integral IR becomes infinitely small (its height equals one, while the width tends to zero).

  )-[START_REF] Fekih-Ahmed | A new effective asymptotic formula for the Stieltjes constants[END_REF] reduce to Jensen's formulas for the ζ function[START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF] Eqs. (88)].

  (45) f (z) = exp (g(z)) h(z), where (46) g(z) = (n + 1) log (log (a + iz)) -2πz, h(z) = (1 + tanh πz) 2 .

Figure 2 .

 2 Figure 2. Real part of f (z) for z = x + 0i on the real line (left) and for z = x + Ci passing near the saddle point (right), here with parameters a = 1 2 , n = 500, where integrating along the real line results in about five digits of cancellation.

  the error is bounded by V c/ρ d where c and ρ only depend on E and [α, β]. If f (z) has poles or branch cuts on E or if the quadrature degree d determined by this bound would have to be larger than O(p)

Theorem 3 . 2

 32 If z is contained in a disk or rectangle Z with midpoint m and radius r, such that Re(Z) ≥ 1, and if max |u-m|≤r |g ′′ (u)|≤ G, then (50) |f (z)| < 4.015 |exp (g(m))| exp |g ′ (m)|r + 1 Gr 2 . Proof. We use the decomposition (45)-(46). If Re(z) ≥ 1, then |h(z)|< 4.015.

4. 2 .

 2 Numerical values. We show the computed values of a few large Stieltjes constants. The following significands are correctly rounded to 100 digits (with at most 0.5 ulp error):γ 10 5 ≈1.9919273063125410956582272431568589205211659777533113258 75975525936171259272227176914320666190965225 • 10 83432 , γ 10 10 ≈ 7.5883621237131051948224033799125486921750410324509700470 54093338492423974783927914992046654518550779 • 10 12397849705 , γ 10 15 ≈ 1.8441017255847322907032695598351364885675746553315587921 86085948502542608627721779023071573732022221 • 10 1452992510427658 , γ 10 100 ≈ 3.1874314187023992799974164699271166513943099108838469225 07106265983048934155937559668288022632306095 • 10 e , e = 2346394292277254080949367838399091160903447689869837 3852057791115792156640521582344171254175433483694.

γ 10 5

 5 (2 + 3i) ≈ (1.52933142489317896667092453331813941673604063614322663 9046917471026123822028695414669890818089958104 + 7.6266053170235392288 29846454534202735013368165330230700751870950104906000791927387438554979 23063058i) • 10 83440 , γ 10 100 (2 + 3i) ≈ (0.02447197253567132691871635713584630519276677767177878 733142765829147799303241971747565188937402242864 + 1.328114485458616967 078662312208319540579816973253179511750642930437359777538176731578318799 940692883i) • 10 e+10 .

Table 1 .

 1 Time in seconds to compute γn. The left columns show results for d digits in Mathematica using N[StieltjesGamma[n],d], or N[StieltjesGamma[n]] when d = 16 giving machine precision. The smallest results are omitted since the timer in Mathematica does not have sufficient resolution. The (wrong) entries signify that Mathematica returns an incorrect result. The (timeout) entries signify that Mathematica had not completed after several hours. The right columns show results for p-bit precision with the new integration algorithm in Arb.

			Mathematica		Arb (integration)
	n	d = 16	d = 100	d = 1000	p = 64	p = 333 p = 3333
	1			0.16	0.0011	0.0089	2.7
	10		0.016	0.39	0.0020	0.032	6.6
	10 2	0.016	0.16	2.7	0.0032	0.030	3.5
	10 3	0.031	0.16	3.3	0.0064	0.10	7.5
	10 4	(wrong)	0.41	4.5	0.0043	0.045	19.8
	10 5	(wrong)	(timeout)	(timeout)	0.0043	0.026	27.8
	10 6	(wrong)			0.0066	0.026	18.1
	10 10				0.0087	0.031	32.6
	10 15				0.014	0.061	7.0
	10 30				0.087	0.22	16.7
	10 60				0.26	0.86	30.9
	10 100				0.76	1.5	57.2

Table 2 .

 2 Time in seconds to compute γ0, . . . , γn simultaneously.

		Arb (Euler-Maclaurin)		Arb (integration)
	n	p = 64	p = 333	p = 3333	p = 64	p = 333	p = 3333
	1	0.000061	0.00026	0.012	0.012	0.12	18
	10	0.00035	0.0016	0.060	0.025	0.20	37
	10 2	0.0047	0.11	0.39	0.28	1.8	370
	10 3	0.69	0.87	5.5	4.3	23	4527
	10 4	1207	1210	1626	38	267
	The table also shows timings for Mathematica 11.0.0 for Microsoft Windows (64-
	bit) on an Intel Core i9-7900X CPU for comparison.

Matsuoka'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 z = x+iy, and then recall that z/|z| 2 = 1/z, where z is the complex conjugate of z.

If we put N = 0 in Matsuoka's expansion[START_REF] Ia | A theorem for the closed-form evaluation of the first generalized Stieltjes constant at rational arguments and some related summations[END_REF], we retrieve, after some calculations and several approximations, the same result as Knessl and Coffey[START_REF] Ia | Three notes on Ser's and Hasse's representations for the zeta-functions[END_REF].

Note that the second form of these expressions is obtained from the former one by a trivial simplification. Moreover, if we remark that coth πx = 1 + 2(e 2πx -1) -1 , we readily notice the relationship between these integrals and the Hermite and Jensen formulas for the ζ functions.

Note that formulas (29), (33) can also provide good computational results.

Acknowledgements

We thank Jacques Gélinas for pointing out the previous computations in [2] and Vladimir Reshetnikov for helping with some numerical verifications, and are especially greateful to Joseph Oesterlé for sharing many challenging ideas on the Stieltjes constants during his stay in St. Petersburg in June 2017.