fork(1) download 
  1. #!/usr/bin/env python3
  2. #
  3. #  pi.py - Calculate Pi
  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. # Global Variables
  14. #
  15. count   = 0
  16. total   = 0
  17. grad    = 0
  18. step    = 0
  19.  
  20. #
  21. # Show Progress
  22. #
  23. def progress_init(max):
  24.         global count, total, grad, step
  25.  
  26.         total = max
  27.         count = 0
  28.         step = int(total / 1000)
  29.         grad = int(step / 2)
  30.  
  31. def progress():
  32.         global count, total, grad, step
  33.  
  34.         if (count > grad):
  35.                 grad += step
  36.                 g = int(math.floor(72.5*count/total+0.5))
  37.                 p = int(math.floor(1000.5*count/total+0.5))
  38.                 msg = "H" * g + "-" * (72-g) + " " + str(p/10) + "%\r"
  39.                 if (grad > total):
  40.                         msg += "\n"
  41.  
  42.                 print(msg, sep="", end="", flush=True)
  43.  
  44. #
  45. # Write digit string
  46. #
  47. def write_string(digit_string):
  48.         fd = open("pi-py.txt", mode="w")
  49.  
  50.         fd.write(" pi = ")
  51.         fd.write(digit_string[0])
  52.         fd.write(".")
  53.  
  54.         for c in range(1, len(digit_string), 50):
  55.                 if (c != 1):
  56.                         fd.write("\t")
  57.  
  58.                 fd.write(digit_string[c:c+50])
  59.  
  60.                 if ((c % 1000) == 951):
  61.                         fd.write(" << ")
  62.                         fd.write(str(c+49))
  63.                         fd.write("\r\n")
  64.                 elif ((c % 500) == 451):
  65.                         fd.write(" <\r\n")
  66.                 else:
  67.                         fd.write("\r\n")
  68.  
  69.         # Final new-line
  70.         fd.write("\r\n")
  71.  
  72.         fd.close()
  73.  
  74. #
  75. # Recursive funcion.
  76. #
  77. def s(a, b, max):
  78.         global count
  79.  
  80.         m = math.ceil((a + b) / 2)
  81.  
  82.         if (b - a == 1):
  83.                 if (a == 0):
  84.                         r = 120         # 6!
  85.                         q = mpz(640320**3)
  86.                         p = gmpy2.sub( gmpy2.mul(q, 13591409),
  87.                                 gmpy2.mul(r, 13591409+545140134) )
  88.                 else:
  89.                         r = mpz(8 * (a*6+1) * (a*6+3) * (a*6+5))
  90.                         q = mpz((b*640320)**3)
  91.                         if ((b%2) == 0):
  92.                                 p = gmpy2.mul(mpz(13591409 + b*545140134), r)
  93.                         else:
  94.                                 p = gmpy2.mul(mpz(-13591409 - b*545140134), r)
  95.         else:
  96.                 p1, q1, r1 = s(a, m, max)
  97.                 p2, q2, r2 = s(m, b, max)
  98.  
  99.                 # Merge
  100.                 p = gmpy2.add( gmpy2.mul(p1, q2), gmpy2.mul(p2, r1) )
  101.                 q = gmpy2.mul(q1, q2)
  102.  
  103.                 if (b != max):
  104.                         r = gmpy2.mul(r1, r2)
  105.                 else:
  106.                         r = 0
  107.  
  108.         count += 1
  109.         progress()
  110.  
  111.         return p, q, r
  112.  
  113. #
  114. # Calculate e
  115. #
  116. def calc_pi(digits):
  117.         global total
  118.  
  119.         d = digits+1
  120.         n_terms = math.ceil(d*math.log(10)/(3*math.log(53360)))
  121.         precision = math.ceil(d * math.log2(10)) + 4
  122.         print("d = ", d, ", n = ", n_terms, ", precision = ", precision, sep="")
  123.  
  124.         print("gmpy2 version:", gmpy2.version())
  125.         print("MP version:", gmpy2.mp_version())
  126.         print("MPFR version:", gmpy2.mpfr_version())
  127.  
  128.         max_precision = gmpy2.get_max_precision()
  129.         print("max_precision =", max_precision)
  130.         max_emax = gmpy2.get_emax_max()
  131.         print("max_emax =", max_emax)
  132.  
  133.         if (max_precision < precision):
  134.                 print("Error! Max precision is too small! Program terminated.")
  135.                 return
  136.  
  137.         gmpy2.get_context().precision = precision
  138.         gmpy2.get_context().emax = max_emax
  139.         print("Real precision = ", gmpy2.get_context().precision)
  140.         progress_init(n_terms * 2 - 1)                          # Initialize progress bar
  141.  
  142.         # Recursion
  143.         start_time = time.monotonic_ns()
  144.         p, q, r = s(0, n_terms, n_terms)
  145.         end_time = time.monotonic_ns()
  146.         elapsed = (end_time - start_time) / 1000000000
  147.         print("Recursion:", elapsed, "seconds.")
  148.  
  149.         start_time = time.monotonic_ns()
  150.         q = gmpy2.mul(q, 426880)
  151.         end_time = time.monotonic_ns()
  152.         elapsed = (end_time - start_time) / 1000000000
  153.         print("Multiply by 426880:", elapsed, "seconds.")
  154.  
  155.         start_time = time.monotonic_ns()
  156.         pf = mpfr(p)
  157.         qf = mpfr(q)
  158.         ef = gmpy2.div(qf, pf)
  159.         end_time = time.monotonic_ns()
  160.         elapsed = (end_time - start_time) / 1000000000
  161.         print("Grand Division:", elapsed, "seconds.")
  162.  
  163.         start_time = time.monotonic_ns()
  164.         ef = gmpy2.mul(ef, gmpy2.sqrt(10005))
  165.         end_time = time.monotonic_ns()
  166.         elapsed = (end_time - start_time) / 1000000000
  167.         print("Multiply by sqrt(10005):", elapsed, "seconds.")
  168.  
  169.         # Convert to decimal digits
  170.         start_time = time.monotonic_ns()
  171.         estr, exp, prec = mpfr.digits(ef)
  172.         estr = estr[0:d]
  173.         end_time = time.monotonic_ns()
  174.         elapsed = (end_time - start_time) / 1000000000
  175.         print("Convert to decimal digits:", elapsed, "seconds.")
  176.  
  177.         # Write to file
  178.         start_time = time.monotonic_ns()
  179.         write_string(estr)
  180.         end_time = time.monotonic_ns()
  181.         elapsed = (end_time - start_time) / 1000000000
  182.         print("Write to file:", elapsed, "seconds.")
  183.  
  184. #
  185. #  main program
  186. #
  187. if __name__ == '__main__':
  188.         argc = len(sys.argv)
  189.         if (argc >= 2):
  190.                 digits = int(sys.argv[1])
  191.         else:
  192.                 digits = 100000
  193.  
  194.         calc_pi(digits)
  195.  
  196. # End of pi.py
  197.  
Runtime error #stdin #stdout #stderr 0.16s 23364KB
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/6YO1zU
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!