2
\$\begingroup\$

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())
\$\endgroup\$
1
\$\begingroup\$

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)
\$\endgroup\$

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.