#!/usr/bin/env python3
#
#  e-mpz.py - Calculate Eular's number e
#
import sys
import time
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
grad    = 0
step    = 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_init(max):
        global count, total, grad, step

        total = max
        count = 0
        step = int(total / 1000)
        grad = int(step / 2)

def progress():
        global count, total, grad, step

        if (count > grad):
                grad += step
                g = int(math.floor(72.5*count/total+0.5))
                p = int(math.floor(1000.5*count/total+0.5))
                msg = "H" * g + "-" * (72-g) + " " + str(p/10) + "%\r"
                if (grad > total):
                        msg += "\n"

                print(msg, sep="", end="", flush=True)

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

        fd.write("  e = ")
        fd.write(digit_string[0])
        fd.write(".")

        for c in range(1, len(digit_string)-1, 50):
                if (c != 1):
                        fd.write("\t")

                fd.write(digit_string[c:c+50])

                if ((c % 1000) == 951):
                        fd.write(" << ")
                        fd.write(str(c+49))
                        fd.write("\r\n")
                elif ((c % 500) == 451):
                        fd.write(" <\r\n")
                else:
                        fd.write("\r\n")

        # Final new-line
        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()
        print("max_precision =", max_precision)
        max_emax = gmpy2.get_emax_max()
        print("max_emax =", max_emax)

        if (max_precision < precision):
                print("Error! Max precision is too small! Program terminated.")
                return

        gmpy2.get_context().precision = precision
        gmpy2.get_context().emax = max_emax
        print("Real precision = ", gmpy2.get_context().precision)
        progress_init(n_terms * 2 - 1)          # Initialize progress bar

        start_time = time.monotonic_ns()
        p, q = s(0, n_terms)
        end_time = time.monotonic_ns()
        elapsed = (end_time - start_time) / 1000000000
        print("Recursion:", elapsed, "seconds.")

        start_time = time.monotonic_ns()
        pf = mpfr(p)
        qf = mpfr(q)
        ef = gmpy2.div(pf, qf)
        end_time = time.monotonic_ns()
        elapsed = (end_time - start_time) / 1000000000
        print("Grand division:", elapsed, "seconds.")

        start_time = time.monotonic_ns()
        estr, exp, prec = mpfr.digits(ef)
        estr = estr[0:d]
        end_time = time.monotonic_ns()
        elapsed = (end_time - start_time) / 1000000000
        print("Convert to decimal digits:", elapsed, "seconds.")

        start_time = time.monotonic_ns()
        write_string(estr)
        end_time = time.monotonic_ns()
        elapsed = (end_time - start_time) / 1000000000
        print("Write file:", elapsed, "seconds.")

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

        calc_e(digits)

# End of e-mpz.py