1
\$\begingroup\$

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

My general methodology is as follows: Compute Gamma via the accepted answer's algorithm here.

In order to do this I need to compute the natural log of a large k.

  1. Compute Pi via the AGM as here.

  2. Compute the natural log of 2 via the last algorithm here.

  3. Compute the natural log of a large k via Gauss's AGM method here.

  4. Compute Gamma as described in the accepted answer here. I am very grateful to the community for sharing the knowledge and I have improved the speed by orders of magnitude compared to the previous version (because of the accepted answers improvements).

One can pip3 install tqdm to see progress, I've wrapped a few ranges with it.

import decimal
#non-builtin
from tqdm import tqdm

D = decimal.Decimal

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, 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)**2)/(2*t)
        k    += 1
        if pi == piold:
            break
    return pi

def factorial(x ,pfact, pfactprd):       
    x = int(x)
    if pfact == (x-1):
        return pfactprd*x
    else:
        factorial = D(1)
        for i in range(1, x+1):
            factorial *= i
        return factorial

def lntwo():        #Fast converging Ln 2
    print('Computing Ln(2)...')
    def lntwosum(n, d, b):
        logsum, logold, e = D(0), D(0), 0
        while 1:
            logold = logsum
            logsum += D(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 lnagm(x):   #Natural log of 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())
    print('Computing Ln(x)...')
    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)
    lnk = lnagm(k)
    upper = int((12*k)+2)
    print('Summing...')
    # First Sum
    logsum = D(0)
    pterm = D((k**2)/2)
    for r in tqdm(range(1, upper)):
        r = D(r)
        logsum += pterm
        nterm = D(((-1)*D(k)*D(r+1))/(r*(r+2)))*pterm
        pterm = nterm
    logsum1 = D(0)
    print('...')
    pfact, pfactprd = 1, 1
    for r in tqdm(range(1, upper)):
        calfact = factorial((r-1), pfact, pfactprd)
        pfact, pfactprd = (r-1), calfact
        logsum1 += D((D(-1)**D(r-1))*(k**D(r+1)))/D(calfact*D(D(r+1)**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('\n')
print(gam)
\$\endgroup\$
  • 2
    \$\begingroup\$ Please pardon my ignorance: What are “functional zeroes” and “binary splitting” (in this context)? – It might also help to understand (and judge) your code if you shortly present (or give links to) the underlying formulae. \$\endgroup\$ – Martin R Sep 6 at 11:41
  • \$\begingroup\$ Pardon my lack of explanation, I was referring to avoiding the Brent-McMillan being computed via the binary splitting method as seen on the last entry here numberworld.org/y-cruncher/internals/… \$\endgroup\$ – TheHoyt Sep 6 at 17:02
2
\$\begingroup\$

Firstly, this code has a couple of dozen PEP8 formatting violations. Following conventions generally helps readability.

Also on the subject of readability, comments providing references for the formulae used should be considered essential in mathematical software.


def agm(a, b):      #Arithmetic Geometric Mean

The conventional way of documenting a function is with a docstring. That also allows you to document the expected input types (decimal.Decimal? float?).


def pi_agm():       #Pi via AGM and lemniscate
    print('Computing Pi...')
    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)**2)/(2*t)
        k    += 1
        if pi == piold:
            break
    return pi

The while loop is unnecessarily ugly: while True would be preferable. I would say that it would be even better to use the loop index as a loop index with for k in itertools.count(), but actually that variable is wholly unnecessary.

If find it unhelpful to initialise six variables in one line where some of them are quite complicated. On the other hand, it could be more helpful to combine some of the updates in the loop body. Perhaps the happy medium is something like

    a, b, t = 1, D(0.5).sqrt(), 1
    p, pi = 1, 0
    while True:
        a, b, t = (a+b)/2, (a*b).sqrt(), t - p*(a-b)**2
        p, piold, pi = 2*p, pi, (a+b)**2 / t

I'm not entirely convinced by

        if pi == piold:
            break

Sometimes iterative approaches in finite data types oscillate around the solution rather than converging definitively. It might be more robust to track the last two or three values and, on finding a loop, return the average of the values in the loop.


def factorial(x ,pfact, pfactprd):       
    x = int(x)

Do you expect to pass a non-int? See previous note about using docstrings to document types.

    if pfact == (x-1):
        return pfactprd*x
    else:
        factorial = D(1)
        for i in range(1, x+1):
            factorial *= i
        return factorial

Is the else ever actually used? Might it be more maintainable to remove this function entirely?


    def lntwosum(n, d, b):
        logsum, logold, e = D(0), D(0), 0
        while 1:
            logold = logsum
            logsum += D(1/((D(b**e))*((2*e)+1)))
            e += 1
            if logsum == logold:
                return (D(n)/D(d))*logsum

Here the previous comment about itertools.count is relevant.

I'm confused as to why D is invoked where it is. Without any comments to justify it, it appears to be done at random.


def lnagm(x):   #Natural log of via AGM,
    try:
        if int(x) == 1:
            return 0
        if int(x) == 2:
            return lntwo()
    except:
        pass

??? Are you expecting int(x) to throw an exception?


    pi = D(pi_agm())

See previous comments about appearing to use D at random. Here, if pi_agm() returns a decimal.Decimal then it's unnecessary, and if it doesn't then surely that would be a bug because pi won't have the necessary precision? I don't see any further polishing of its error.


    decimal.getcontext().prec = m
    ...
    logr = D(str(result)[:m//2])
    decimal.getcontext().prec = prec
    return logr

Would the following work?

    decimal.getcontext().prec = m
    ...
    decimal.getcontext().prec = prec
    return D(result)

    k = D(prec//2)
    ...
    upper = int((12*k)+2)

Why not just upper = 6*prec + 2 with no need to coerce?


Other earlier comments also apply to gamma.


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

It's a Python best practice, which serves to make the file reusable as a library, to guard this with if __name__ == "__main__":.

A comment explaining the *8 would be useful. At a guess, prec is in bits?

\$\endgroup\$
  • \$\begingroup\$ Prec in bits- correct. Thank you for corrections. I do realize this is quite a pep felon. I appreciate you pointing out the type and functional redundancies. You have given me a good guide to clean up this code. \$\endgroup\$ – TheHoyt Sep 6 at 16:54

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.