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.