Nth digit of PI

Peter Schneider-Kamp petersc at stud.ntnu.no
Tue Jun 6 03:09:23 EDT 2000


Kirby Urner wrote:
> 
> The original algorithm worked for PI expressed as
> a hexadecimal number, but I found a web page purporting
> to give the algorithm for a base 10 PI, but I couldn't
> figure out how to make it work in Python.[1][2]

Can you give the algorithm in some kind of pseudo code?
Pseudo python would be appreciated ...

If I understood it right the calculation has to be done
like this for every pair of primes in the denominator:

1) calculate k1 and k2 for 1/(a*b) = k1/a+k2/b
2) calculate m1, m2, m3 for (k1/a+k2/b)/c = m1/a+m2/b+m3/c

So that means one basically needs a continued fraction
algorithm and a factorization algorithm, right?

thinking-the-article-is-not-really-clear-about-it-ly
y'rs Peter

BTW: The algorithm for base 16 works ... very slowly:

def hexdcm(a,b,d,e):
  x = long(a % b)
  if d:
    for k in range(d):
      y = x << 4
      x = y % b
  else:
    y = a
  c = [y-x]
  for j in range(e):
    y = x << 4
    x = y % b
    c.append(y-x)
  sum = 0
  for j in range(e+1):
    sum = sum + (1L << ((e-j) << 2)) * (c[j] / b)
  return sum

def sum(d,e):
  sum = 0L
  for n in range(d+1):
    sum = sum + hexdcm(120L*n**2L+151L*n+47L,\
(8L*n+1L)*(2L*n+1L)*(8L*n+5L)*(4L*n+3L),d-n,e)
  return sum 

def pi(d,e):
  s = sum(d,e) % 16L**(e+1L)
  return int((s - s % 16L**e) / 16L**e)

def upto(n):
  for i in range(n+1):
    if not i % 50:
      if i:
        print s+"\n  ",
      else:
        print "3.",
      s = ""
    s = s + "%x" % pi(i+1,4)
  print "\n"
--
Peter Schneider-Kamp          ++47-7388-7331
Herman Krags veg 51-11        mailto:peter at schneider-kamp.de
N-7050 Trondheim              http://schneider-kamp.de




More information about the Python-list mailing list