fork download
copy 
  1. #
  2. #  pi.py - Calculate Pi
  3. #
  4. import sys
  5. import math
  6. import gmpy2
  7. from gmpy2 import mpfr
  8. from gmpy2 import mpz
  9.  
  10. #
  11. # Global Variables
  12. #
  13. count = 0
  14. total = 0
  15. old_p = 0
  16.  
  17. #
  18. # Show Progress
  19. #
  20. def progress():
  21.         global count, old_p, total
  22.  
  23.         p = int(math.floor(1000*count/total+0.5))
  24.         if (p > old_p):
  25.                 old_p = p
  26.                 g = int(math.floor(72.5*count/total+0.5))
  27.                 for c in range(72):
  28.                         if (c<g):
  29.                                 print("H", sep="", end="")
  30.                         else:
  31.                                 print("-", sep="", end="")
  32.                 print(" ", p/10, "%\r", sep="", end="", flush=True)
  33.  
  34.         if (count == total):
  35.                 print("\n", sep="", end="")
  36.  
  37. #
  38. # Write digit string
  39. #
  40. def write_string(digit_string):
  41.         fd = open("pi-py.txt", mode="w")
  42.  
  43.         fd.write(" pi = ")
  44.         for c in range(len(digit_string)):
  45.                 if ((c != 1) and (c % 50 == 1)):
  46.                         fd.write("\t")
  47.                 fd.write(digit_string[c])
  48.                 if (c == 0):
  49.                         fd.write(".")
  50.                 elif ((c % 1000) == 0):
  51.                         fd.write(" << ")
  52.                         fd.write(str(c))
  53.                         fd.write("\r\n")
  54.                 elif ((c % 500) == 0):
  55.                         fd.write(" <\r\n")
  56.                 elif ((c % 50) == 0):
  57.                         fd.write("\r\n")
  58.                 elif ((c % 5) == 0):
  59.                         fd.write(" ")
  60.  
  61.         # Final new-line
  62.         if ((c%50) != 0):
  63.                 fd.write("\r\n")
  64.  
  65.         fd.close()
  66.  
  67. #
  68. # Recursive funcion.
  69. #
  70. def s(a, b, max):
  71.         global count
  72.  
  73.         m = math.ceil((a + b) / 2)
  74.  
  75.         if (b - a == 1):
  76.                 if (a == 0):
  77.                         r = 120         # 6!
  78.                         q = mpz(640320**3)
  79.                         p = gmpy2.sub( gmpy2.mul(q, 13591409),
  80.                                 gmpy2.mul(r, 13591409+545140134) )
  81.                 else:
  82.                         r = mpz(8 * (a*6+1) * (a*6+3) * (a*6+5))
  83.                         q = mpz((b*640320)**3)
  84.                         if ((b%2) == 0):
  85.                                 p = gmpy2.mul(mpz(13591409 + b*545140134), r)
  86.                         else:
  87.                                 p = gmpy2.mul(mpz(-13591409 - b*545140134), r)
  88.         else:
  89.                 p1, q1, r1 = s(a, m, max)
  90.                 p2, q2, r2 = s(m, b, max)
  91.  
  92.                 # Merge
  93.                 p = gmpy2.add( gmpy2.mul(p1, q2), gmpy2.mul(p2, r1) )
  94.                 q = gmpy2.mul(q1, q2)
  95.  
  96.                 if (b != max):
  97.                         r = gmpy2.mul(r1, r2)
  98.                 else:
  99.                         r = 0
  100.  
  101.         count += 1
  102.         progress()
  103.  
  104.         return p, q, r;
  105.  
  106. #
  107. # Calculate e
  108. #
  109. def calc_pi(digits):
  110.         global total
  111.  
  112.         d = digits+1
  113.         n_terms = math.ceil(d*math.log(10)/(3*math.log(53360)))
  114.         precision = math.ceil(d * math.log2(10)) + 4;
  115.         print("d = ", d, ", n = ", n_terms, ", precision = ", precision, sep="")
  116.  
  117.         gmpy2.get_context().precision = precision
  118.         gmpy2.get_context().emax = 1000000000000;
  119.         print("Real precision =", gmpy2.get_context().precision)
  120.         total = n_terms * 2 - 1         # Max progress
  121.  
  122.         p, q, r = s(0, n_terms, n_terms)
  123.  
  124.         q = gmpy2.mul(q, 426880)
  125.  
  126.         print("Recursion done. Processing division.")
  127.         ef = gmpy2.div(q, p)
  128.  
  129.         ef = gmpy2.mul(ef, gmpy2.sqrt(10005))
  130.  
  131.         print("Division done. Converting to digits.")
  132.         estr, exp, prec = mpfr.digits(ef)
  133.         estr = estr[0:d]
  134.         print("Writing to file.")
  135.         write_string(estr);
  136.  
  137. #
  138. #  main program
  139. #
  140. argc = len(sys.argv)
  141. if (argc >= 2):
  142.         digits = int(sys.argv[1])
  143. else:
  144.         digits = 100000
  145.  
  146. calc_pi(digits)
  147.  
  148. # End of pi.py
  149.  
Runtime error #stdin #stdout #stderr 0.13s 23540KB
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/jAiHdM
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!