bbp_pi.py 5.9 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190
  1. '''
  2. This implementation is a heavily modified fixed point implementation of
  3. BBP_formula for calculating the nth position of pi. The original hosted
  4. at: https://web.archive.org/web/20151116045029/http://en.literateprograms.org/Pi_with_the_BBP_formula_(Python)
  5. # Permission is hereby granted, free of charge, to any person obtaining
  6. # a copy of this software and associated documentation files (the
  7. # "Software"), to deal in the Software without restriction, including
  8. # without limitation the rights to use, copy, modify, merge, publish,
  9. # distribute, sub-license, and/or sell copies of the Software, and to
  10. # permit persons to whom the Software is furnished to do so, subject to
  11. # the following conditions:
  12. #
  13. # The above copyright notice and this permission notice shall be
  14. # included in all copies or substantial portions of the Software.
  15. #
  16. # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  17. # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  18. # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  19. # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  20. # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  21. # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  22. # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  23. Modifications:
  24. 1.Once the nth digit and desired number of digits is selected, the
  25. number of digits of working precision is calculated to ensure that
  26. the hexadecimal digits returned are accurate. This is calculated as
  27. int(math.log(start + prec)/math.log(16) + prec + 3)
  28. --------------------------------------- --------
  29. / /
  30. number of hex digits additional digits
  31. This was checked by the following code which completed without
  32. errors (and dig are the digits included in the test_bbp.py file):
  33. for i in range(0,1000):
  34. for j in range(1,1000):
  35. a, b = pi_hex_digits(i, j), dig[i:i+j]
  36. if a != b:
  37. print('%s\n%s'%(a,b))
  38. Deceasing the additional digits by 1 generated errors, so '3' is
  39. the smallest additional precision needed to calculate the above
  40. loop without errors. The following trailing 10 digits were also
  41. checked to be accurate (and the times were slightly faster with
  42. some of the constant modifications that were made):
  43. >> from time import time
  44. >> t=time();pi_hex_digits(10**2-10 + 1, 10), time()-t
  45. ('e90c6cc0ac', 0.0)
  46. >> t=time();pi_hex_digits(10**4-10 + 1, 10), time()-t
  47. ('26aab49ec6', 0.17100000381469727)
  48. >> t=time();pi_hex_digits(10**5-10 + 1, 10), time()-t
  49. ('a22673c1a5', 4.7109999656677246)
  50. >> t=time();pi_hex_digits(10**6-10 + 1, 10), time()-t
  51. ('9ffd342362', 59.985999822616577)
  52. >> t=time();pi_hex_digits(10**7-10 + 1, 10), time()-t
  53. ('c1a42e06a1', 689.51800012588501)
  54. 2. The while loop to evaluate whether the series has converged quits
  55. when the addition amount `dt` has dropped to zero.
  56. 3. the formatting string to convert the decimal to hexadecimal is
  57. calculated for the given precision.
  58. 4. pi_hex_digits(n) changed to have coefficient to the formula in an
  59. array (perhaps just a matter of preference).
  60. '''
  61. from sympy.utilities.misc import as_int
  62. def _series(j, n, prec=14):
  63. # Left sum from the bbp algorithm
  64. s = 0
  65. D = _dn(n, prec)
  66. D4 = 4 * D
  67. d = j
  68. for k in range(n + 1):
  69. s += (pow(16, n - k, d) << D4) // d
  70. d += 8
  71. # Right sum iterates to infinity for full precision, but we
  72. # stop at the point where one iteration is beyond the precision
  73. # specified.
  74. t = 0
  75. k = n + 1
  76. e = D4 - 4 # 4*(D + n - k)
  77. d = 8 * k + j
  78. while True:
  79. dt = (1 << e) // d
  80. if not dt:
  81. break
  82. t += dt
  83. # k += 1
  84. e -= 4
  85. d += 8
  86. total = s + t
  87. return total
  88. def pi_hex_digits(n, prec=14):
  89. """Returns a string containing ``prec`` (default 14) digits
  90. starting at the nth digit of pi in hex. Counting of digits
  91. starts at 0 and the decimal is not counted, so for n = 0 the
  92. returned value starts with 3; n = 1 corresponds to the first
  93. digit past the decimal point (which in hex is 2).
  94. Parameters
  95. ==========
  96. n : non-negative integer
  97. prec : non-negative integer. default = 14
  98. Returns
  99. =======
  100. str : Returns a string containing ``prec`` digits
  101. starting at the nth digit of pi in hex.
  102. If ``prec`` = 0, returns empty string.
  103. Raises
  104. ======
  105. ValueError
  106. If ``n`` < 0 or ``prec`` < 0.
  107. Or ``n`` or ``prec`` is not an integer.
  108. Examples
  109. ========
  110. >>> from sympy.ntheory.bbp_pi import pi_hex_digits
  111. >>> pi_hex_digits(0)
  112. '3243f6a8885a30'
  113. >>> pi_hex_digits(0, 3)
  114. '324'
  115. These are consistent with the following results
  116. >>> import math
  117. >>> hex(int(math.pi * 2**((14-1)*4)))
  118. '0x3243f6a8885a30'
  119. >>> hex(int(math.pi * 2**((3-1)*4)))
  120. '0x324'
  121. References
  122. ==========
  123. .. [1] http://www.numberworld.org/digits/Pi/
  124. """
  125. n, prec = as_int(n), as_int(prec)
  126. if n < 0:
  127. raise ValueError('n cannot be negative')
  128. if prec < 0:
  129. raise ValueError('prec cannot be negative')
  130. if prec == 0:
  131. return ''
  132. # main of implementation arrays holding formulae coefficients
  133. n -= 1
  134. a = [4, 2, 1, 1]
  135. j = [1, 4, 5, 6]
  136. #formulae
  137. D = _dn(n, prec)
  138. x = + (a[0]*_series(j[0], n, prec)
  139. - a[1]*_series(j[1], n, prec)
  140. - a[2]*_series(j[2], n, prec)
  141. - a[3]*_series(j[3], n, prec)) & (16**D - 1)
  142. s = ("%0" + "%ix" % prec) % (x // 16**(D - prec))
  143. return s
  144. def _dn(n, prec):
  145. # controller for n dependence on precision
  146. # n = starting digit index
  147. # prec = the number of total digits to compute
  148. n += 1 # because we subtract 1 for _series
  149. # assert int(math.log(n + prec)/math.log(16)) ==\
  150. # ((n + prec).bit_length() - 1) // 4
  151. return ((n + prec).bit_length() - 1) // 4 + prec + 3