# Arbitrary precision Euler-Mascheroni constant via a Brent-McMillan algorithm with no math module

Utilizing the below relation, I am able to compute the Euler constant to great precision on a single thread quickly and simply.

My process is to compute the natural log via the AGM(utilizing Pi an the Ln(2)), and apply the result to the BMM algorithm.

$$\\gamma=\lim_{n\to\infty}\mathscr{G}_n=\lim_{n\to\infty}\frac{\sum\limits_{k=0}^\infty \left(\frac{n^k}{k!}\right)^2 (H_k-\log\,n)}{\sum\limits_{k=0}^\infty \left(\frac{n^k}{k!}\right)^2}\$$

where $$\H_k=\sum\limits_{j=1}^k \frac1{j}\$$ is a Harmonic number.

This process is much more efficient than my previous post.

import decimal
from tqdm import tqdm

def agm(a, b):
"""Arithmetic Geometric Mean"""
a, b = D(a),D(b)
for x in tqdm(range(prec)):
a, b = (a + b) / 2, (a * b).sqrt()
return a

def pi_agm():
"""Pi via AGM and lemniscate"""
print('Computing Pi...')
a, b, t = 1, D(0.5).sqrt(), 1/D(2)
p, pi, k = 2, 0, 0
while 1:
an    = (a+b)/2
b     = (a*b).sqrt()
t    -= p*(a-an)**2
a, p  = an, 2**(k+2)
piold = pi
pi    = ((a+b)**2)/(2*t)
k    += 1
if pi == piold:
break
return pi

def lnagm(x):
"""Natural log of via AGM"""
if x == D(1):
return 0
if x == D(2):
return lntwo()
m = prec*2
ln2 = lntwo()
decimal.getcontext().prec = m
pi = pi_agm()
print('Computing Ln({0})...'.format(x))
twoprec = (D(2)**(2-m))/x
den = agm(1, twoprec)*2
diff = m*ln2
result = ((pi/den) - diff)
logr = D(str(result)[:m//2])
decimal.getcontext().prec = prec
return logr

def lntwo():        #Fast converging Ln 2
print('Computing Ln(2)...')
def lntwosum(n, d, b):
logsum, logold, e = D(0), D(0), 0
while True:
logold = logsum
logsum += (1/((D(b**e))*((2*e)+1)))
e += 1
if logsum == logold:
return (D(n)/D(d))*logsum
logsum1 = lntwosum(14, 31, 961)
logsum2 = lntwosum(6, 161, 25921)
logsum3 = lntwosum(10, 49, 2401)
ln2 = logsum1 + logsum2 + logsum3
return ln2

def gamma():
n = D(prec)
a = u = -lnagm(n)
b = v = D(1)
i = D(1)
while True:
k = (n/i)**2
a *= k
b *= k
a += (b/i)
if (u+a) == u or (v+b) == v:
break
u += a
v += b
i += 1
return D(str(u/v)[:prec-3])

if __name__ == "__main__":
"""Setting precision plus 5 to accomodate error"""
prec = int(input('Precision: '))+5
decimal.getcontext().prec = prec
D = decimal.Decimal
print(gamma())


The single largest optimization I was able to make was to the logarithms. Once the two terms of a AGM are the same, they will not differ again. Therefore you can be sure you have converged. The argument of oscillating convergence does not apply.

I also noticed that, code aside, y=2^x+8; AGM(1, y) converges the fastest out of the numbers. I don't know the exact Math behind AGM convergence rates but I noticed that.(therefore in the script, put precision as 2^x -1 for ultra fast Ln convergence)

Since this is my own answer, I'm fine with me putting gmpy2 in places to make the computations faster.

I computed 212,143 digits of Gamma in 3.5 hours with the script. Usable precision into the thousands is done in less than breath :) .

import decimal
from tqdm import tqdm
from gmpy2 import mpc
import gmpy2

def agm(a, b):
"""Arithmetic Geometric Mean"""
a, b = D(a),D(b)
for i in tqdm(range(prec)):
a, b = (a + b) / 2, (a * b).sqrt()
if a == b:
break
return a

def pi_agm():
"""Pi via AGM and lemniscate"""
print('Computing Pi...')
a, b, t = 1, D(0.5).sqrt(), 1/D(2)
p, pi, k = 2, 0, 0
while 1:
an    = (a+b)/2
b     = (a*b).sqrt()
t    -= p*(a-an)**2
a, p  = an, 2**(k+2)
piold = pi
pi    = ((a+b)**2)/(2*t)
k    += 1
if pi == piold:
break
return pi

def lnagm(x):
"""Natural log of via AGM"""
if x == D(1):
return 0
if x == D(2):
return lntwo()
m = prec*2
ln2 = lntwo()
decimal.getcontext().prec = m
pi = pi_agm()
print('Coverging on Ln({0})...'.format(x))
twoprec = (D(2)**(2-m))/x
den = agm(1, twoprec)*2
diff = m*ln2
result = ((pi/den) - diff)
logr = D(str(result)[:m//2])
decimal.getcontext().prec = prec
return logr

def lntwo():
"""Fast converging Ln2"""
print('Computing Ln(2)...')
def lntwosum(n, d, b):
logsum, logold = mpc(0), mpc(0),
e, n, d = mpc(0), mpc(n), mpc(d)
while True:
logold = logsum
logsum += (1/((b**e)*((2*e)+1)))
e += 1
if logsum == logold:
return (n/d)*logsum
logsum1 = lntwosum(14, 31, 961)
logsum2 = lntwosum(6, 161, 25921)
logsum3 = lntwosum(10, 49, 2401)
ln2 = logsum1 + logsum2 + logsum3
return D(str(ln2.real)[:-4])

def gamma():
n = D(prec)
a = u = mpc(str(-lnagm(n)))
print('Computing Gamma...')
b = v = mpc(1)
i = mpc(1)
while True:
k = (n/i)**2
a *= k
b *= k
a += (b/i)
if (u+a) == u or (v+b) == v:
break
u += a
v += b
i += 1
return D(str(u/v)[:prec-3])

if __name__ == "__main__":
"""Setting precision plus 5 to accomodate error"""
prec = int(input('Precision: '))+5
gmpy2.get_context().precision=(prec*4)
decimal.getcontext().prec = prec
D = decimal.Decimal
gam = gamma()
print(gam)