An asymptotic expansion for the error term in the Brent-McMillan algorithm for Euler's constant

The Brent-McMillan algorithm is the fastest known procedure for the high-precision computation of Euler's constant $\gamma$ and is based on the modified Bessel functions $I_0(2x)$ and $K_0(2x)$. An error estimate for this algorithm relies on the optimally truncated asymptotic expansion for the product $I_0(2x) K_0(2x)$ when $x$ assumes large positive integer values. An asymptotic expansion for this optimal error term is derived by exploiting the techniques developed in hyperasymptotics, thereby enabling more precise information on the error term than recently obtained bounds and estimates.


Introduction
The Brent-McMillan algorithm [2] (when implemented with binary splitting) is the fastest known method of high-precision computation of Euler's constant γ. This relies on the formula [6, (10.31. 3 where throughout we take x to be a positive integer, I 0 (x), K 0 (x) are the standard modified Bessel functions and H n (n!) 2 ( 1 2 x) 2n , H n := 1 + 1 2 + · · · + 1 n .
For large x, the final term in (1.1) is O(e −4x ). Greater precision can be achieved following the suggestion made in [2] of optimally truncating the asymptotic expansion corresponding to the truncation index k = 2x, followed by computing the term K 0 (2x)/I 0 (2x) from I 0 (2x)K 0 (2x)/(I 0 (2x)) 2 . In [3], a bound was obtained for the remainder term in the optimally truncated expansion (1.2) given by 24e −8x , thereby providing rigour to the algorithm. More recently, Demailly [4] established the leading large-x behaviour of this remainder, together with an error bound, in the form This leads to the error in the optimally truncated expansion of the final term in (1.1) given by −5 √ 2πe −8x /(12x 1/2 ) to leading order. The problem with using the well-known asymptotic expansions of the modified Bessel functions is that the positive real axis is a Stokes line for I 0 (x) (but not for K 0 (x)). The standard expansion [6, (10.40.5)], [9, p. 203] is clearly inadequate, since this predicts a purely imaginary exponentially small contribution as x → +∞ when clearly it must be real. The correct form of expansion of I 0 (x) for x → +∞ that takes into account the Stokes phenomenon on the positive x-axis has been considered in [7]. In this paper we derive an asymptotic expansion for the remainder in the optimally truncated expansion (1.2) by applying the first stage of the hyperasymptotic expansion process (also known as exponential improvement) to a suitable integral representation for the product I 0 (x)K 0 (x). A discussion of the new theory of hyperasymptotics (initiated by Berry [1]) can be found in the book [8,Ch. 6] in the context of the confluent hypergeometric functions; see also [6,Section 2.11]. We present some numerical results to illustrate the accuracy of the expansion so obtained.
Throughout this paper we shall restrict x to be a positive integer in keeping with the strategy of the Brent-McMillan algorithm, although the analysis can be developed for complex x. The integrand has simple poles situated at s = 0, 1, 2, . . . and double poles at s = − 1 2 , − 3 2 , . . . . Displacement of the integration path to the right over the first N poles, where for the moment N is an arbitrary integer, together with the fact that the residue of Γ(−s) at s = k is (−1) k−1 /k!, then shows that with L N denoting the rectilinear path (−c+N −∞i, −c+N +∞i) and 0 < c < 1.
We now choose N to be the optimal truncation index of the expansion in (2.1) (corresponding to truncation at, or near, the least term), which is easily verified to be N = x. As a consequence, since x → +∞ the variable s in the quotient of gamma functions in (2.2) is uniformly large on the displaced path L N . From Lemma 2.2 in [8, p. 39] we have the inverse factorial expansion We now introduce the so-called terminant function T ν (z) defined 1 as a multiple of the incomplete gamma function Γ(a, z) by From the formula connecting Γ(a, ze ±πi ) given in [ The Mellin-Barnes integral representation of this function is [8, (6 provided ν = 0, −1, −2, . . . . Then, if we make the change of variable s → s + N in the integrals appearing in (2.5), write cos 2πs in terms of exponentials, and use (2.7) (when it is supposed that M < 2N ) these integrals can be written as upon application of (2.6). This then yields the expansion for R N (x) given by It now remains to exploit the known asymptotic expansions of the terminant function T ν (x) when ν ∼ x as x → +∞, which we do in the next section.
3. An asymptotic expansion for R N (x) The asymptotic expansion of the terminant function T ν (z) for large ν and complex z, when ν ∼ |z|, has been discussed in detail by Olver in [5]; see also [6, Section 2.11] and the detailed account in [8, pp. 259-265]. By expressing T ν (z) in terms of the Laplace integral Olver established by application of the saddle-point method that when µ ∼ x (and bounded integer j) where A 0,j = 1 (j ≥ 0) and On the negative real axis, where a saddle point and a simple pole become coincident in the above Laplace integral, we have the expansion where the coefficients G k,j result from the expansion The branch of w(τ ) is chosen such that w ∼ τ − 1 as τ → 1. Upon reversion of the w-τ mapping to yield τ = 1 + w + 1 3 w 2 + 1 36 w 3 − 1 270 w 4 + 1 4320 w 5 + · · · , it is found with the help of Mathematica that the first five even-order coefficients G 2k,j ≡ 6 −2kĜ for x → +∞. Then we obtain the following theorem: the coefficients c k are specified in (2.4) and the quantities A k,j and G k,j ≡ ( 1 2 ) k 2 k G k,j are defined in (3.2) and (3.5).
Since µ = 2N and the variables in the terminant functions in (2.8) involve 2x, we have from This produces the expansion as x → +∞, which is the main result of the paper. In Table 1 we present values of the absolute relative error in the computation of the above expansion for R N (x) compared with the exact evaluation from (2.1) for different x and truncation index k.  1.020 × 10 −9 4 3.555 × 10 −9 1.137 × 10 −10 1.510 × 10 −11

Concluding remarks
In the Brent-McMillan algorithm we have, with N = x, In [4], Demailly defined his remainder function ∆(x) using the optimal truncation index k = 2N , instead of k = 2N − 1, and wrote The connection with our R 2N (2x) is consequently given by Then, from (4.1), the error resulting from K 0 (2x)/I 0 (2x) in (1.1) at optimal truncation has the expansion Appendix: Determination of the coefficients c j in the expansion (2.3) Use of the duplication formula Γ(2s) = π −1/2 2 2s−1 Γ(s)Γ(s + 1 2 ) shows that the inverse factorial expansion (2.3) can be written as  Use of the well-known expansion [6, (5.11