4
\$\begingroup\$

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)
\$\endgroup\$
  • \$\begingroup\$ From purely mathematical point of view, you might be interested in this post, which has a collection of different algorithms for this constant \$\endgroup\$ – Yuriy S Sep 5 at 16:35
  • \$\begingroup\$ 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). \$\endgroup\$ – TheHoyt Sep 5 at 16:41
6
\$\begingroup\$
  • 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.

\$\endgroup\$
  • \$\begingroup\$ 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. \$\endgroup\$ – TheHoyt Sep 5 at 17:04

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.