Advertisement
jukaukor

BaileyBorweinPlouffePi_non_numba.py

Aug 5th, 2019
120
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.00 KB | None | 0 0
  1.  
  2. # Python program computes nth hex digit of Pi with BBP formula.
  3. # BBP = Bailey-Borwein-Plouffe_formula for Pi
  4. # (without numba compiler & njit directive)
  5. # Juhani Kaukoranta 5.8.2018
  6.  
  7. import math
  8. import time
  9.  
  10.  
  11. # Convert Pi to Hex-strings
  12. # (float number's decimals to hex digits)
  13.  
  14. def convHex(x) :
  15. hx = "0123456789ABCDEF"
  16. y = math.fabs(x)
  17. chx =''
  18. for i in range(0,16):
  19. y = 16.0 * (y - math.floor(y))
  20. chx += hx[int(y)]
  21. return chx
  22.  
  23. # Compute Modular Exponentiation
  24.  
  25. def compModExp(b, e, m):
  26. if (e == 0):
  27. return 1
  28. ans = compModExp(b, e // 2, m)
  29. ans = (ans * ans) % m
  30. if ((e % 2) == 1):
  31. ans = (ans * b) % m
  32. return ans
  33.  
  34. # Compute Sums
  35.  
  36. def Sums(j):
  37. s = 0.0
  38. # t // Each term of right summation
  39. # r // Denominator
  40. # k // Loop index
  41. # EPS = Loop-exit accuration of the right summation
  42.  
  43. EPS = 1.0e-17
  44. # Left Sum (0 ... n)
  45. k = 0
  46. for k in range(0,n+1):
  47. r = 8 * k + j
  48. t = compModExp(16, n - k, r)
  49. t /= r
  50. s += t - int(t)
  51. s -= int(s)
  52.  
  53. #Right sum (n + 1 ...)
  54. k = n+1
  55. while (1) :
  56. r = 8 * k + j
  57. t = 16**(n - k)
  58. t /= r
  59. if (t < EPS):
  60. break
  61. s += t
  62. s -= int(s)
  63. k += 1
  64.  
  65. return s # left + right
  66.  
  67. # paaohjelma kysyy halutun hex-muotoisen desimaalin järjestysnumeroa
  68. # tulostaa sen ja laskennassa käytetyn desimaalimuotoisen fraktion
  69.  
  70. n = int(input("kuinka mones Piin desimaali halutaan laskea?"))
  71. print("**** PI Computation ( hex digit number: ", n,", waiting..." )
  72. # Compute PI
  73. n = n - 1
  74. t0 = time.perf_counter() # time start
  75. pi = 4 * Sums(1) - 2 * Sums(4) - Sums(5) - Sums(6)
  76. pi = pi - int(pi) + 1
  77. pi_hex = convHex(pi)
  78. print("Hex digit: ", pi_hex[0])
  79. print("10 Hex digit: ",pi_hex[0:10])
  80. print("FRACTION : ", pi)
  81. t1 = time.perf_counter() # time end
  82. print("Aikaa kului ", t1 - t0," sekuntia")
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement