Arbitrary precision Euler-Mascheroni constant via a Brent-McMillan algorithm with no math module 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())
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)