# Calculate the Euler-Mascheroni constant without the math module

The below code was written to generate γ, for educational purposes.

Single threaded, no functional zeroes required, no binary splitting(which can all be used to compute competitively like y-cruncher, that version in works). Uses the Arithmetic Geometric Mean to calculate large logarithms quickly. Uses decimal module for precision management.

I've computed 3000 digits in a few hours with it, and 200 in about a minute. I am happy calculating.

import decimal

D = decimal.Decimal

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

def pi_agm():       #Pi via AGM and lemniscate
a, b, t, p, pi, k = 1, D(2).sqrt()/2, 1/D(2), 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) * (a + b) / (2*t)
k    += 1
if pi == piold:
break
return pi

def factorial(x):       #factorial fast loop
x = int(x)
factorial = D(1)
for i in range(1, x+1):
factorial *= i
return factorial

def lntwo():        #Fast converging Ln 2
logsum, logold, n = D(0), D(0), 0
while 1:
logold = logsum
logsum += D(1/((D(961**n))*((2*n)+1)))
n += 1
if logsum == logold:
logsum1 = (D(14)/D(31))*logsum
break
logsum, logold, n = D(0), D(0), 0
while 1:
logold = logsum
logsum += D(1/((D(25921**n))*((2*n)+1)))
n += 1
if logsum == logold:
logsum2 = (D(6)/D(161))*logsum
break
logsum, logold, n = D(0), D(0), 0
while 1:
logold = logsum
logsum += D(1/((D(2401**n))*((2*n)+1)))
n += 1
if logsum == logold:
logsum3 = (D(10)/D(49))*logsum
break
ln2 = logsum1 + logsum2 + logsum3
return ln2

def lnagm(x):   #Natural log via AGM,
try:
if int(x) == 1:
return 0
if int(x) == 2:
return lntwo()
except:
pass
m = prec*2
ln2 = lntwo()
decimal.getcontext().prec = m
pi = D(pi_agm())
twoprec = D(2**(2-D(m)))/D(x)
den = agm(1, twoprec)*2
diff = m*ln2
result = (D(pi/den) - D(diff))
logr = D(str(result)[:m//2])
decimal.getcontext().prec = prec
return logr

def gamma():   #Compute Gamma from Digamma Expansion
print('Computing Gamma!')
k = D(prec/2)
print('Calculating Logarithms...')
lnk = lnagm(k)
logsum = D(0)
upper = int((12*k)+2)
print('Summing...')
for r in range(1, upper):
logsum += D((D(-1)**D(r-1))*D(k**D(r+1)))/D(factorial(r-1)*D(r+1))
if r%1000==0:
print(str((D(r)/D(upper))*100)[:5], '% ; Sum 1 of 2')
logsum1 = D(0)
print('...')
for r in range(1, upper):
logsum1 += D((D(-1)**D(r-1))*(k**D(r+1)))/D(factorial(r-1)*D(D(r+1)**2))
if r%1000==0:
print(str((D(r)/D(upper))*100)[:5], '% ; Sum 2 of 2')
twofac = D(2)**(-k)
gammac = str(D(1)-(lnk*logsum)+logsum1+twofac)
return D(gammac[:int(prec//6.66)])

#Calling Gamma
prec = int(input('Precision for Gamma: '))*8
decimal.getcontext().prec = prec
gam = gamma()
print(gam)

• From purely mathematical point of view, you might be interested in this post, which has a collection of different algorithms for this constant – Yuriy S Sep 5 at 16:35
• That is excellent, thank you! I need to learn to write latex. The top answer is indeed the algorithm I wrote here :) I believe the only way I can go faster, is to use the algorithm involving a(log(a)-1)=1 and binary splitting.(ignoring its python). – TheHoyt Sep 5 at 16:41

• Computation of logsum and logsum1 in gamma() are suboptimal. You do costly operations of raising to power, and recompute factorial on each iteration (the latter invokes the quadratic time complexity BTW). Notice that in the $$\\sum \dfrac{(-1)^{r-1} k^{r+1}}{(r+1)(r-1)!}\$$ a consecutive term can be expressed via the previous one, as$$\T_{r+1} = -k\dfrac{r+1}{(r+2)r} T_n\$$. Instead of computing each term from scratch, use this recurrence, and enjoy a significant performance boost.
• Tree loops in lntwo cry to become a function.