Calculate the Euler-Mascheroni constant without the math module

Question

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)

• Yuriy S: From purely mathematical point of view, you might be interested in this post, which has a collection of different algorithms for this constant
• FodderOverflow: 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).

Answer 1

• 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.

  Converting the summation into a Horner schema would likely improve the accuracy, that is to achieve the desired number of digits you'd need less number of terms.

• Tree loops in lntwo cry to become a function.

• FodderOverflow: Thank you! I will implement everything you have said. However I believe using the recurrence and/or the Horner schema will prevent me from multi-threading in the future, just as an FYI to readers.