fork download 
  1. #!/usr/bin/env python3
  2. #
  3. #  e-mpz.py - Calculate Eular's number e
  4. #
  5. import sys
  6. import time
  7. import math
  8. import gmpy2
  9. from gmpy2 import mpfr
  10. from gmpy2 import mpz
  11.  
  12. #
  13. # Constants used in Stirling's approximation
  14. #
  15. E       = float(2.718281828459045235360287)
  16. pi      = float(3.141592653589793238462643)
  17. C       = math.log10(2*pi) / 2
  18.  
  19. #
  20. # Global Variables
  21. #
  22. count   = 0
  23. total   = 0
  24. grad    = 0
  25. step    = 0
  26.  
  27. #
  28. # Stirling's approximation
  29. #
  30. def logfactorial(n):
  31.         return (C + math.log10(n)/2 + n*(math.log10(n)-math.log10(E)))
  32.  
  33. #
  34. # Estimate how many terms in the serie sould be calculated.
  35. #
  36. def terms(digits):
  37.         upper = 2
  38.         lower = 1
  39.         while (logfactorial(upper)<digits):
  40.                 upper <<= 1
  41.         else:
  42.                 lower = upper/2
  43.  
  44.         while ((upper-lower) > 1):
  45.                 n = (upper+lower)/2
  46.                 if (logfactorial(n) > digits):
  47.                         upper = n
  48.                 else:
  49.                         lower = n
  50.  
  51.         return n
  52.  
  53. #
  54. # Show Progress
  55. #
  56. def     progress_init(max):
  57.         global count, total, grad, step
  58.  
  59.         total = max
  60.         count = 0
  61.         step = int(total / 1000)
  62.         grad = int(step / 2)
  63.  
  64. def progress():
  65.         global count, total, grad, step
  66.  
  67.         if (count > grad):
  68.                 grad += step
  69.                 g = int(math.floor(72.5*count/total+0.5))
  70.                 p = int(math.floor(1000.5*count/total+0.5))
  71.                 msg = "H" * g + "-" * (72-g) + " " + str(p/10) + "%\r"
  72.                 if (grad > total):
  73.                         msg += "\n"
  74.  
  75.                 print(msg, sep="", end="", flush=True)
  76.  
  77. #
  78. # Write digit string
  79. #
  80. def write_string(digit_string):
  81.         fd = open("e-py.txt", mode="w")
  82.  
  83.         fd.write("  e = ")
  84.         fd.write(digit_string[0])
  85.         fd.write(".")
  86.  
  87.         for c in range(1, len(digit_string)-1, 50):
  88.                 if (c != 1):
  89.                         fd.write("\t")
  90.  
  91.                 fd.write(digit_string[c:c+50])
  92.  
  93.                 if ((c % 1000) == 951):
  94.                         fd.write(" << ")
  95.                         fd.write(str(c+49))
  96.                         fd.write("\r\n")
  97.                 elif ((c % 500) == 451):
  98.                         fd.write(" <\r\n")
  99.                 else:
  100.                         fd.write("\r\n")
  101.  
  102.         # Final new-line
  103.         fd.write("\r\n")
  104.  
  105.         fd.close()
  106.  
  107. #
  108. # Recursive funcion.
  109. #
  110. def s(a, b):
  111.         global count
  112.  
  113.         m = math.ceil((a + b) / 2)
  114.  
  115.         if (a == b):
  116.                 q = mpz(1)
  117.                 if (a == 0):
  118.                         p = mpz(1)
  119.                 else:
  120.                         p = mpz(0)
  121.         elif (b - a == 1):
  122.                 if (a == 0):
  123.                         p = mpz(2)
  124.                         q = mpz(1)
  125.                 else:
  126.                         p = mpz(1)
  127.                         q = mpz(b)
  128.         else:
  129.                 p1, q1 = s(a, m)
  130.                 p2, q2 = s(m, b)
  131.  
  132.                 # Merge
  133.                 p = gmpy2.add(gmpy2.mul(p1, q2), p2)
  134.                 q = gmpy2.mul(q1, q2)
  135.  
  136.         count += 1
  137.         progress()
  138.  
  139.         return p, q
  140.  
  141. #
  142. # Calculate e
  143. #
  144. def calc_e(digits):
  145.         global total
  146.  
  147.         d = digits+1
  148.         n_terms = int(terms(d))
  149.         precision = math.ceil(d * math.log2(10)) + 4
  150.         print("d = ", d, ", n = ", n_terms, ", precision = ", precision)
  151.  
  152.         print("gmpy2 version:", gmpy2.version())
  153.         print("MP version:", gmpy2.mp_version())
  154.         print("MPFR version:", gmpy2.mpfr_version())
  155.  
  156.         max_precision = gmpy2.get_max_precision()
  157.         print("max_precision =", max_precision)
  158.         max_emax = gmpy2.get_emax_max()
  159.         print("max_emax =", max_emax)
  160.  
  161.         if (max_precision < precision):
  162.                 print("Error! Max precision is too small! Program terminated.")
  163.                 return
  164.  
  165.         gmpy2.get_context().precision = precision
  166.         gmpy2.get_context().emax = max_emax
  167.         print("Real precision = ", gmpy2.get_context().precision)
  168.         progress_init(n_terms * 2 - 1)          # Initialize progress bar
  169.  
  170.         start_time = time.monotonic_ns()
  171.         p, q = s(0, n_terms)
  172.         end_time = time.monotonic_ns()
  173.         elapsed = (end_time - start_time) / 1000000000
  174.         print("Recursion:", elapsed, "seconds.")
  175.  
  176.         start_time = time.monotonic_ns()
  177.         pf = mpfr(p)
  178.         qf = mpfr(q)
  179.         ef = gmpy2.div(pf, qf)
  180.         end_time = time.monotonic_ns()
  181.         elapsed = (end_time - start_time) / 1000000000
  182.         print("Grand division:", elapsed, "seconds.")
  183.  
  184.         start_time = time.monotonic_ns()
  185.         estr, exp, prec = mpfr.digits(ef)
  186.         estr = estr[0:d]
  187.         end_time = time.monotonic_ns()
  188.         elapsed = (end_time - start_time) / 1000000000
  189.         print("Convert to decimal digits:", elapsed, "seconds.")
  190.  
  191.         start_time = time.monotonic_ns()
  192.         write_string(estr)
  193.         end_time = time.monotonic_ns()
  194.         elapsed = (end_time - start_time) / 1000000000
  195.         print("Write file:", elapsed, "seconds.")
  196.  
  197. #
  198. #  main program
  199. #
  200. if __name__ == '__main__':
  201.         argc = len(sys.argv)
  202.         if (argc >= 2):
  203.                 digits = int(sys.argv[1])
  204.         else:
  205.                 digits = 100000
  206.  
  207.         calc_e(digits)
  208.  
  209. # End of e-mpz.py
  210.  
Runtime error #stdin #stdout #stderr 0.15s 23400KB
comments (?)
stdin 
Standard input is empty
stdout 
Standard output is empty
stderr 
Traceback (most recent call last):
  File "./prog.py", line 8, in <module>
ModuleNotFoundError: No module named 'gmpy2'
https://ideone.com/3B6wdb
language:
Python 3 (python  3.12)
created:
4 years ago
visibility:
public

Share or Embed source code

Discover > Sphere Engine API

The brand new service which powers Ideone!

Discover > IDE Widget

Widget for compiling and running the source code in a web browser!