Python pi calculation?
It seems you are losing precision in this line:
pi = pi * Decimal(12)/Decimal(640320**(1.5))
Try using:
pi = pi * Decimal(12)/Decimal(640320**Decimal(1.5))
This happens because even though Python can handle arbitrary scale integers, it doesn't do so well with floats.
Bonus
A single line implementation using another algorithm (the BBP formula):
from decimal import Decimal, getcontextgetcontext().prec=100print sum(1/Decimal(16)**k * (Decimal(4)/(8*k+1) - Decimal(2)/(8*k+4) - Decimal(1)/(8*k+5) - Decimal(1)/(8*k+6)) for k in range(100))
For people who come here just to get a ready solution to get arbitrary precision of pi with Python (source with a couple of edits):
import decimaldef pi(): """ Compute Pi to the current precision. Examples -------- >>> print(pi()) 3.141592653589793238462643383 Notes ----- Taken from https://docs.python.org/3/library/decimal.html#recipes """ decimal.getcontext().prec += 2 # extra digits for intermediate steps three = decimal.Decimal(3) # substitute "three=3.0" for regular floats lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24 while s != lasts: lasts = s n, na = n + na, na + 8 d, da = d + da, da + 32 t = (t * n) / d s += t decimal.getcontext().prec -= 2 return +s # unary plus applies the new precisiondecimal.getcontext().prec = 1000pi = pi()
from decimal import *#Sets decimal to 25 digits of precisiongetcontext().prec = 25def factorial(n): if n<1: return 1 else: return n * factorial(n-1)def plouffBig(n): #http://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula pi = Decimal(0) k = 0 while k < n: pi += (Decimal(1)/(16**k))*((Decimal(4)/(8*k+1))-(Decimal(2)/(8*k+4))-(Decimal(1)/(8*k+5))-(Decimal(1)/(8*k+6))) k += 1 return pidef bellardBig(n): #http://en.wikipedia.org/wiki/Bellard%27s_formula pi = Decimal(0) k = 0 while k < n: pi += (Decimal(-1)**k/(1024**k))*( Decimal(256)/(10*k+1) + Decimal(1)/(10*k+9) - Decimal(64)/(10*k+3) - Decimal(32)/(4*k+1) - Decimal(4)/(10*k+5) - Decimal(4)/(10*k+7) -Decimal(1)/(4*k+3)) k += 1 pi = pi * 1/(2**6) return pidef chudnovskyBig(n): #http://en.wikipedia.org/wiki/Chudnovsky_algorithm pi = Decimal(0) k = 0 while k < n: pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k))) k += 1 pi = pi * Decimal(10005).sqrt()/4270934400 pi = pi**(-1) return piprint "\t\t\t Plouff \t\t Bellard \t\t\t Chudnovsky"for i in xrange(1,20): print "Iteration number ",i, " ", plouffBig(i), " " , bellardBig(i)," ", chudnovskyBig(i)