DeepAI

# Computing Stieltjes constants using complex integration

The Stieltjes constants γ_n are the coefficients appearing in the Laurent series of the Riemann zeta function at s = 1. We give a simple and efficient method to compute a p-bit approximation of γ_n with rigorous error bounds. Starting from an integral representation due to Blagouchine, we shift the contour to eliminate cancellation. The integral is then evaluated numerically in ball arithmetic using the Petras algorithm, with the use of a Taylor expansion for bounds near the saddle point. This appears to be the first algorithm for Stieltjes constants with uniformly low complexity with respect to both n and p. An implementation is provided in the Arb library. We can, for example, compute γ_n to 1000 digits in a minute for any n < 10^100.

02/22/2018

### Numerical integration in arbitrary-precision ball arithmetic

We present an implementation of arbitrary-precision numerical integratio...
06/18/2018

### Numerical Evaluation of Elliptic Functions, Elliptic Integrals and Modular Forms

We describe algorithms to compute elliptic functions and their relatives...
06/03/2019

### A library to compute the density of the distance between a point and a random variable uniformly distributed in some sets

In [3], algorithms to compute the density of the distance to a random va...
05/09/2017

### Computing the Lambert W function in arbitrary-precision complex interval arithmetic

We describe an algorithm to evaluate all the complex branches of the Lam...
05/21/2020

### An improved Neumann series of Bessel functions representation for solutions of perturbed Bessel equations

New representations for an integral kernel of the transmutation operator...
08/08/2022

### Numerical Evaluation of Mittag-Leffler Functions

The Mittag-Leffler function is computed via a quadrature approximation o...
10/10/2017

### Representations and evaluation strategies for feasibly approximable functions

A famous result due to Ko and Friedman (1982) asserts that the problems ...

## 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) ζ(s,v)=1s−1+∞∑n=0(−1)nn!γn(v)(s−1)n,s∈C∖{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.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 ,

 (2) γn(v)=−π2(n+1)∫+∞−∞logn+1(v−12+ix)cosh2πxdx,

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

 (3) γn(v)=(12v−logvn+1)lognv−i∫∞0dxe2πx−1{logn(v−ix)v−ix−logn(v+ix)v+ix},

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) ζ(s)=limN→∞{N∑k=11ks−N1−s1−s},Re(s)>0,

which follows from the Euler–Maclaurin summation formula, so that

 γn=limN→∞{N∑k=1lognkk−logn+1Nn+1},n∈N0,

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.222The 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 .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 [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],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 [23], [27, p. 3]. Four years later he derived the complete asymptotic expansion

 (5) γn∼n!eg(b)πN∑k=0|h2k|2k+12Γ(k+12)(g′′(b)2+f′′(b)2)12k+14××cos[f(b)−(k+12)arctanf′′(b)g′′(b)+arctanv2ku2k],N=0,1,2,…

where , and are the sequences of numbers defined by

 ∞∑k=0hk(y−b)k=exp[ϕ(a+iy)−ϕ(a+ib)+12ϕ′′(a+ib)(y−b)2], uk≡Re(hk),vk≡Im(hk),

and and are the functions defined as

 g(y)≡Reϕ(a+iy),f(y)≡Imϕ(a+iy), ϕ(z)=−(n+1)logz−zlog2πi+logΓ(z).

The pair , , is the unique solution of the equation

 (6) dϕ(z)dz=−n+1z−log2πi+Ψ(z)=0,

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 [19] 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).

 (7) γn∼B√nenAcos(an+b)

in terms of the slowly varying functions

 A=12log(α2+β2)−αα2+β2,B=2√2π(α2+β2)4√(α+1)2+β2,
 a=arctanβα+βα2+β2,b=arctanβα−12arctanβα+1,

where is the unique solution of

 2πexp(βtanβ)=ncosββ,with 0<β<12π,α=βtanβ.

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) ζ(s,v)=ζ(s,v+1)+1vs,γn(v)=γn(v+1)+lognvv.
###### Theorem 1.

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

 (11) ζ(s,v) =π2(s−1)∫+∞−∞(v−12±ix)1−scosh2πxdx =π2(s−1)∫∞0(v−12−ix)1−s+(v−12+ix)1−scosh2πxdx =π2(s−1)∫∞0cos[(s−1)arctan2x2v−1](v2−v+14+x2)12(s−1)cosh2πxdx

and

 (13) γn(v) =−π2(n+1)∫+∞−∞logn+1(v−12±ix)cosh2πxdx =−π2(n+1)∫∞0logn+1(v−12−ix)+logn+1(v−12+ix)cosh2πxdx

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 [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) \ointctrclockwiseC(a−iz)1−scosh2πzdz=∫+R−R(a−ix)1−scosh2πxdx+∫CR(a−iz)1−scosh2πzdz.

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

 ∣∣∣∫CR(a−iz)1−scosh2πzdz∣∣∣=R∣∣ ∣∣∫π0(a−iReiφ)1−seiφcosh2(πReiφ)dφ∣∣ ∣∣≤ ≤Rmaxφ∈[0,π]∣∣(a−iReiφ)1−s∣∣⋅IR (15) ≤Rmaxφ∈[0,π][|a|2+2R(axsinφ−aycosφ)+R2]12Re(1−s)eπ|Im(1−s)|IR

where we denoted , and

 IR≡∫π0dφ∣∣cosh(πReiφ)∣∣2,R>0,

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

 1∣∣cosh(πReiφ)∣∣2=2cosh(2πRcosφ)+cos(2πRsinφ).

Since and are both real, except for the case when . Hence

 cosh(2πRcosφ)+cos(2πRsinφ)>0,

except perhaps at . But at the latter point, since is integer,

 cosh(2πRcosφ)+cos(2πRsinφ)∫∣∣∣φ=12π=1+cos(2πR)=2.

and when we have

 (16) 1∣∣cosh(πReiφ)∣∣2={O(e−2πRcosφ),φ∈[0,12π],O(e+2πRcosφ),φ∈[12π,π].

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

 (17) IR =∫π02dφcosh(2πRcosφ)+cos(2πRsinφ) =∫π204dφcosh(2πRcosφ)+cos(2πRsinφ) =O(∫π20e−2πRcosφdφ)=O(∫π20e−2πRsinϑdϑ),R→∞.

From the inequality

it follows that

 (18) 1−e−π2R2πR≤∫π20e−2πRsinϑdϑ≤1−e−2πR4R,

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

(19)
i.e. the integral asymptotically tends to the left bound (18). Inserting this result into (15), we obtain

 (20) ∣∣∣∫CR(a−iz)1−scosh2πzdz∣∣∣→0asR→∞,R∈N,

if . Hence, making , equality (14) becomes

 (21) ∫+∞−∞(a−ix)1−scosh2πxdx=\ointctrclockwiseC(a−iz)1−scosh2πzdz

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) \ointctrclockwiseC(a−iz)1−scosh2πzdz=2πi∞∑n=0resz=zn(a−iz)1−scosh2πz= =−2iπ∞∑n=0∂∂z(a−iz)1−s∣∣∣z=i(n+12)= =2(s−1)π∞∑n=0(a+12+n)−s=2(s−1)πζ(s,a+12).

Equating (21) with the last result yields

 (23) ζ(s,a+12)=π2(s−1)∫+∞−∞(a−ix)1−scosh2πxdx,Re(a)>0.

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

 (24) (a+ix)s+(a−ix)s=2(a2+x2)s2cos(sarctanxa)

the latter expression may also be written as

 (26) ζ(s,a+12) =π2(s−1)∫∞0(a+ix)1−s+(a−ix)1−scosh2πxdx =πs−1∫∞0cos[(s−1)arctanxa](a2+x2)12(s−1)cosh2πxdx,Re(a)>0.

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

 (s−1)ζ(s,a+12)=1+∞∑n=0(−1)nγn(a+12)n!(s−1)n+1,s∈C∖{1}.

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) ζ(s,v) =(v−12)1−ss−1+i∫∞0(v−12−ix)−s−(v−12+ix)−se2πx+1dx =(v−12)1−ss−1−2∫∞0sin[sarctan2x2v−1](v2−v+14+x2)s2⋅dxe2πx+1

and

 (29) γn(v)=−logn+1(v−12)n+1++i∫∞0⎧⎨⎩logn(v−12−ix)v−12−ix−logn(v−12+ix)v−12+ix⎫⎬⎭dxe2πx+1

respectively.

###### Proof.

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

 (30) ∫∞0f(x)cosh2πxdx=f(0)π+2π∫∞0f′(x)e2πx+1dx,

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

 (31) (a+ix)s−(a−ix)s=2i(a2+x2)s2sin(sarctanxa),

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) γn(v)=−πn+1Re⎧⎨⎩∫∞0logn+1(v−12±ix)cosh2πxdx⎫⎬⎭

and to

 (33) γn(v) =−logn+1(v−12)n+1±2Im⎧⎨⎩∫∞0logn(v−12±ix)v−12±ix⋅dxe2πx+1⎫⎬⎭

respectively.

###### Remark 2.

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

 (35) ζ(s) =11−2s−1{12+12i∫+∞−∞(1−ix)−sdxsinhπx}= =11−2s−1{12+∫∞0sin(sarctanx)(1+x2)s2sinhπxdx},
 (37) ζ(s,v) =±iπ2(s−1)(s−2)∫+∞−∞(v−12±ix)2−ssinhπxcosh3πxdx= =2π2(s−1)(s−2)∫∞0sin[(s−2)arctan2x2v−1](v2−v+14+x2)s2−1⋅sinhπxcosh3πxdx,
 (38) ζ(s,v)=3π3(s−1)(s−2)(s−3)∫+∞−∞{1cosh4πx−23cosh2πx}dx(v−12±ix)s−3,
 (39) γ1=π224−γ22−log222+log2π2−log2⋅logπ+∫∞0arctanx⋅log(1+x2)sinhπxdx,
 γn(v)=(−1)n+14π{n(n−1)ζ(n−2)(3,v)+3nζ(n−1)(3,v)+2ζ(n−2)(3,v)}− (40) −3π4(n+1)∫+∞−∞logn+1(v−12±ix)cosh4πxdx,

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,

 Ψ(v)=−Ψ2(v)4π2+3π4∫∞0log(v2−v+14+x2)cosh4πxdx, logΓ(v)=12log2π+(v−12)(Ψ(v)−1)−π∫∞0xarctan2x2v−1cosh2πxdx, logΓ(v)=12log2π+(v−12)(Ψ(v)+Ψ2(v)4π2−1)− −Ψ1(v)4π2−3π2∫∞0xarctan2x2v−1cosh4πxdx,

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 [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

 Φ(z,s,v)=v−s2+i2∫∞0(−z)ix(v+ix)−s−(−z)−ix(v−ix)−ssinhπxdx =v−s2+∫∞0cos(xlogz)sin[sarctanxv]−sin(xlogz)cos[sarctanxv](v2+x2)s2tanhπxdx,

valid for , or

 ζ(s,v) =v−s2+i2∫∞0e−πx(v+ix)−s−e+πx(v−ix)−ssinh