#
#  e.py - Calculate Eular's number e
#
import sys
import math
import gmpy2
from gmpy2 import mpfr
from gmpy2 import mpz

#
# Constants used in Stirling's approximation
#
E       = float(2.718281828459045235360287)
pi      = float(3.141592653589793238462643)
C       = math.log10(2*pi) / 2

#
# Global Variables
#
count = 0
total = 0
old_p = 0

#
# Stirling's approximation
#
def logfactorial(n):
        return (C + math.log10(n)/2 + n*(math.log10(n)-math.log10(E))) ;

#
# Estimate how many terms in the serie sould be calculated.
#
def terms(digits):
        upper = 2;
        lower = 1;
        while (logfactorial(upper)<digits):
                upper <<= 1
        else:
                lower = upper/2;

        while ((upper-lower) > 1):
                n = (upper+lower)/2
                if (logfactorial(n) > digits):
                        upper = n;
                else:
                        lower = n;

        return n

#
# Show Progress
#
def progress():
        global count, old_p, total

        p = int(math.floor(1000*count/total+0.5))
        if (p > old_p):
                old_p = p
                g = int(math.floor(72.5*count/total+0.5))
                for c in range(72):
                        if (c<g):
                                print("H", sep="", end="")
                        else:
                                print("-", sep="", end="")
                print(" ", p/10, "%\r", sep="", end="", flush=True)

        if (count == total):
                print("\n", sep="", end="")

#
# Write digit string
#
def write_string(digit_string):
        fd = open("e-py.txt", mode="w")

        fd.write("  e = ")
        for c in range(len(digit_string)):
                if ((c != 1) and (c % 50 == 1)):
                        fd.write("\t")
                fd.write(digit_string[c])
                if (c == 0):
                        fd.write(".")
                elif ((c % 1000) == 0):
                        fd.write(" << ")
                        fd.write(str(c))
                        fd.write("\r\n")
                elif ((c % 500) == 0):
                        fd.write(" <\r\n")
                elif ((c % 50) == 0):
                        fd.write("\r\n")
                elif ((c % 5) == 0):
                        fd.write(" ")

        # Final new-line
        if ((c%50) != 0):
                fd.write("\r\n")

        fd.close()

#
# Recursive funcion.
#
def s(a, b):
        global count

        m = math.ceil((a + b) / 2)

        if (a == b):
                q = mpz(1)
                if (a == 0):
                        p = mpz(1)
                else:
                        p = mpz(0)
        elif (b - a == 1):
                if (a == 0):
                        p = mpz(2)
                        q = mpz(1)
                else:
                        p = mpz(1)
                        q = mpz(b)
        else:
                p1, q1 = s(a, m)
                p2, q2 = s(m, b)

                # Merge
                p = gmpy2.add(gmpy2.mul(p1, q2), p2)
                q = gmpy2.mul(q1, q2)

        count += 1
        progress()

        return p, q;

#
# Calculate e
#
def calc_e(digits):
        global total

        d = digits+1
        n_terms = int(terms(d))
        precision = math.ceil(d * math.log2(10)) + 4;
        print("d = ", d, ", n = ", n_terms, ", precision = ", precision)

        print("gmpy2 version:", gmpy2.version())
        print("MP version:", gmpy2.mp_version())
        print("MPFR version:", gmpy2.mpfr_version())
        max_precision = gmpy2.get_max_precision()
        if (max_precision < precision):
                print("Error! max precision", max_precision,
                        "is smaller than required precision", precision)
                return
        gmpy2.get_context().precision = precision
        print("Real precision = ", gmpy2.get_context().precision)

        emax = gmpy2.get_emax_max()
        print("gmpy2 emax max =", emax)
        gmpy2.get_context().emax = emax;

        total = n_terms * 2 - 1         # Max progress

        p, q = s(0, n_terms)

        pf = mpfr(p);
        qf = mpfr(q);
        ef = gmpy2.div(p, q)
        estr, exp, prec = mpfr.digits(ef)
        estr = estr[0:d]
        write_string(estr);

#
#  main program
#
argc = len(sys.argv)
if (argc >= 2):
        digits = int(sys.argv[1])
else:
        digits = 100000

calc_e(digits)

# End of e.py