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