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