Python pi calculation?

Question

I am a python beginner and I want to calculate pi. I tried using the Chudnovsky algorithm because I heard that it is faster than other algorithms.

This is my code:

from math import factorial
from decimal import Decimal, getcontext

getcontext().prec=100

def calc(n):
    t= Decimal(0)
    pi = Decimal(0)
    deno= Decimal(0)
    k = 0
    for k in range(n):
        t = ((-1)**k)*(factorial(6*k))*(13591409+545140134*k)
        deno = factorial(3*k)*(factorial(k)**3)*(640320**(3*k))
        pi += Decimal(t)/Decimal(deno)                                   
    pi = pi * Decimal(12)/Decimal(640320**(1.5))
    pi = 1/pi
    return pi

print calc(25)

For some reason this code yields the vakue of pi up to only 15 decimals as compared with the acceptable value. I tried to solve this by increasing the precision value; this increases the number of digits, but only the first 15 are still accurate. I tried changing the way it calculates the algorithm and it didn't work either. So my question is, is there something that can be done to this code to make it much more accurate or would I have to use another algorithm? I would appreciate help with this because I don't know how to operate with so many digits in python. I would like to be able to control the number of (correct) digits determined and displayed by the program -- whether 10, 100, 1000, etc.

Answer 1

from decimal import *

#Sets decimal to 25 digits of precision
getcontext().prec = 25

def 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 pi

def 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 pi

def 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 pi
print "\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)

Answer 2

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 decimal

def 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 precision

decimal.getcontext().prec = 1000
pi = pi()

Answer 3

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, getcontext
getcontext().prec=100
print 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))