sympy/ntheory/bbp_pi.py

   1 '''
   2 This implimentation is a heavily modified fixed point implimentation of BBP_formula
   3 for calculating the nth position of pi. The original hosted at:
   4 http://en.literateprograms.org/Pi_with_the_BBP_formula_(Python)
   5
   6 # Permission is hereby granted, free of charge, to any person obtaining
   7 # a copy of this software and associated documentation files (the
   8 # "Software"), to deal in the Software without restriction, including
   9 # without limitation the rights to use, copy, modify, merge, publish,
  10 # distribute, sublicense, and/or sell copies of the Software, and to
  11 # permit persons to whom the Software is furnished to do so, subject to
  12 # the following conditions:
  13 #
  14 # The above copyright notice and this permission notice shall be
  15 # included in all copies or substantial portions of the Software.
  16 #
  17 # THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
  18 # EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
  19 # MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT.
  20 # IN NO EVENT SHALL THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY
  21 # CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT,
  22 # TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
  23 # SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
  24
  25
  26 Modifications:
  27 1.Once the nth digit is selected the number of digits of working precision is
  28 caluclated to ensure that the 14 Hexidecimal representation of that region is
  29 accurate. This was found eimperically to be int((math.log10(n//1000))+18).
  30 This was found by searching for a value of working precision for the n = 0 and
  31 n = 1 then n was increased untill the result was less precise, therefore increased again
  32 this was repeated for increasing n and an effective fit was found between n and
  33 the working precision value.
  34
  35 2. The while loop to evaluate whether the series has convergered has be replaced
  36 with a fixed for loop, that option was selected because in a very large number of
  37 cases the loop convereged to a point where no difference can be detected in less than
  38 15 iterations. (done for more accurate memory and time banking).
  39
  40 3. output hex string constrained to 14 characters (accuarcy assured to n = 10**7)
  41
  42 4. pi_hex_digits(n) changed to have coeffiecient to the formula in an array (perhaps just
  43 a matter of preference).
  44
  45 '''
  46 import math
  47 def Series(j, n):
  48
  49     # Left sum from the bbp algorithm
  50     s = 0
  51     D = dn(n)
  52     for k in range(0, n+1):
  53         r = 8*k+j
  54         s = (s + (pow(16,n-k,r)<<4*(D))//r)
  55
  56     # Right sum. should iterate to infinty, but now just iterates to the point where
  57     # one iterations change is beyond the resolution of the data type used
  58
  59     t = 0
  60     for k in range(n+1, n+15):
  61         xp = int(16**(n-k) * (16**(D))  )
  62         t = t + xp // (8*k+j)
  63     total = s+t
  64
  65     return total
  66
  67 def pi_hex_digits(n):
  68
  69     # main of implimentation arrays holding formulae coefficients
  70     n -= 1
  71     a=  [4,2,1,1]
  72     j = [1,4,5,6]
  73
  74     #formulae
  75     x =  + (a[0]*Series(j[0], n)
  76          - a[1]*Series(j[1], n)
  77          - a[2]*Series(j[2], n)
  78          - a[3]*Series(j[3], n)) & (16**(dn(n)) -1)
  79
  80     s=("%014x" % x)
  81     #s is constrained between 0 and 14
  82     return s[0:14]
  83
  84 def dn(n):
  85     # controller for n dependence on precision
  86     if (n < 1000):
  87         f=16
  88     else:
  89         f = int((math.log10(n//1000))+18)
  90     return f
  91
  92
  93