# Mathematica's algorithm for computing $\gamma$

Please note that this question is not a duplicate of What is the fastest/most efficient algorithm for estimating Euler's Constant $\gamma$?

I was looking for algorithms for computing $$\gamma$$ (Euler's constant) and various representations of $$\gamma$$, and I noticed something strange:

On this page (https://functions.wolfram.com/Constants/EulerGamma/09/), look all the way down; you'll find this formula: $$\gamma=\lim_{x\to\infty}\left(\frac{1}{\sum_{k=0}^\infty \frac{x^k}{k!^2}}\sum_{k=0}^\infty \sum_{i=1}^k \frac{x^k}{ik!^2}-\frac{\log x}{2}\right)$$ Allegedly, the above formula provides an algorithm for computing $$\gamma$$ to arbitrary precision. They say

The above formula is used for the numerical computation of Euler-Mascheroni constant in Mathematica. The algorithm which is based on this formula is the fastest known algorithm for computing this constant.

Now this was quite astounding for me when I first read it.

Sure, it approaches $$\gamma$$, but how does one use it for computing $$\gamma$$? If you take the above formula and choose some $$x$$, how can you bound the error between $$\gamma$$ and your approximation?

I tried squeezing the formula for some "alternating character" that would somehow give a sequence of lower bounds and a sequence of upper bounds, but I failed.

• Apparently , this series converges very quickly so in practice we will have to run $k$ to only a moderate number to get a very precise result. Commented Apr 13 at 9:59
• I also wonder how it can be proven that this strange sum (Is Ramanujan still alive?) actually converges to $\gamma$. Commented Apr 13 at 10:00

You may find this formula mentioned in the MathWorld article for the Euler-Mascheroni constant; it is Equation $$(44)$$. According to the article, it is based on Equation $$(38)$$, and the cited references are:

Brent, R. P. and McMillan, E. M. "Some New Algorithms for High-Precision Computation of Euler's Constant." Math. Comput. 34, 305-312, 1980

Trott, M. The Mathematica GuideBook for Programming. New York: Springer-Verlag, 2004. http://www.mathematicaguidebooks.org/

I should also point out that in many cases, Mathematica leverages efficient sub-algorithms in its internal computations. So for instance, if the algorithm for $$\gamma$$ requires a computation of harmonic numbers, rather than a specific implementation, it could simply call one that was programmed for standalone use.

• I'm reading Brent's paper, page 307. Can you explain what does "The iteration (10) is terminated when, to the working precision, $U_k=U_{k-1}$ and $V_k=V_{k-1}$." mean in mathematical terms? (As clearly $U_k\ne U_{k-1}$ and $V_k\ne V_{k-1}$ in the strict sense). How do you bound the error between $U_k(n)$ and $U(n)$, etc.? Commented Apr 13 at 10:43
• @Nomas2 I haven't read the paper as I don't have time to study it in detail presently, but the idea presumably is to stop computing additional steps in a recursion when the difference between the current step and the previous step falls below a certain threshold, much like how if you wanted to compute a series representation, say $$e = \sum_{k=0}^\infty \frac{1}{k!}$$ to some number of decimal digits of precision, say $10^{-6}$, you would stop once the denominator $k!$ grew large enough. Commented Apr 13 at 10:46
• So the "working precision" here is a kind of algorithmically imposed error tolerance. It could be a machine precision bound, so in the case of floating-point arithmetic, it might be around the order of $10^{-16}$; but because the algorithm is intended to be used to compute $\gamma$ to arbitrary precision, the "working precision" here is a function of the number of correct digits of $\gamma$ you want to get. Commented Apr 13 at 10:50