Ways to calculate Pi (in Python)

Enjoy (requires gmpy)!

UPDATE: More Taylor series approximations (I love them), polygon method now works.

"""
    Calculate pi.
"""

import math, cmath, gmpy
import time

# Euler
machin_euler_set = ((5, 7), (2, gmpy.mpf(79, 1024) / 3))
# F. C. W. Stormer
machin_1896_set = ((44, 57), (7, 239), (-12, 682), (24, 12943))
# Kikuo Takano
machin_1982_set = ((12, 49), (32, 57), (-5, 239), (12, 110443))
# Hwang Chien-Lih
machin_2003_set = ((183, 239), (32, 1023), (-68, 5832), (12, 113021), (-100, 6826318), (-12, 33366019650), (12, 43599522992503626068))

def factorial(n):
    if (n == 0):
        return 1
    
    if (n <= 1):
        return n
    
    return n * factorial(n - 1)

def double_factorial(n):
    if (n == 0) or (n == -1) or (n == 1):
        return 1
    
    return n * double_factorial(n - 2)

def taylor_sin(x, precision, n = 100):
    """ Return the taylor series approximation of sine. """
    gmpy.set_minprec(precision)
    sum = gmpy.mpf(0)
    
    for i in xrange(n):
        sum += pow(-1, i) * (gmpy.mpf(x) ** ((2 * i) + 1)) / factorial((2 * i) + 1)
    
    return sum

def taylor_cos(x, precision, n = 100):
    """ Return the taylor series approximation of cosine. """
    gmpy.set_minprec(precision)
    sum = gmpy.mpf(0)
    
    for i in xrange(n):
        sum += pow(-1, i) * (gmpy.mpf(x) ** (2 * i)) / factorial(2 * i)
    
    return sum

def taylor_atan(x, precision, n = 100):
    """ Return the taylor series approximation of arctan. """
    gmpy.set_minprec(precision)
    sum = gmpy.mpf(0)
    
    for i in xrange(n):
        sum += pow(-1, i) * (gmpy.mpf(x) ** ((2 * i) + 1)) / ((2 * i) + 1)
    
    return sum

def calcpi_polygon(sides, precision):
    """ The classic polygon method. """
    gmpy.set_minprec(precision)
    tmp = gmpy.mpf(2.0) - (2.0 * taylor_cos(((gmpy.mpf(360) / sides) * (gmpy.pi(precision) / 180.0)), precision, precision / 4))
    
    return gmpy.mpf(sides) / (gmpy.mpf(2) * gmpy.fsqrt(gmpy.mpf(1) / tmp))

def calcpi_leibniz(n, precision):
    """ The Leibniz method. Converges very slowly. """
    gmpy.set_minprec(precision)
    sum = gmpy.mpf(0)
    
    for i in xrange(n):
        sum += pow(-1, i) * gmpy.mpf(1) / gmpy.mpf((2 * i) + 1)
    
    return sum * 4

def calcpi_newton(n, precision):
    """ The Newton method. """
    gmpy.set_minprec(precision)
    sum = gmpy.mpf(0)
    
    for i in xrange(n):
        sum += gmpy.mpf(factorial(i)) / gmpy.mpf(double_factorial((2 * i) + 1))
    
    return sum * 2

def calcpi_madhava(n, precision):
    """ The Madhava method. """
    gmpy.set_minprec(precision)
    sum = gmpy.mpf(0)
    
    for i in xrange(n):
        sum += pow(-1, i) * (gmpy.mpf(1) / (((2 * i) + 1) * pow(3, i)))
    
    return sum * gmpy.fsqrt(12)

def calcpi_machin(n, precision, set):
    """ The Machin method. Can use multiple sets """
    gmpy.set_minprec(precision)
    sum = gmpy.mpf(0)
    
    for a, b in set:
        sum += gmpy.mpf(a) * taylor_atan(gmpy.mpf(1) / gmpy.mpf(b), precision, precision / 4)
    
    return sum * 4

def calcpi_chudnovsky(n, precision):
    """ The Chudnovsky method. Converges so fast you can't believe it. """
    gmpy.set_minprec(precision)
    sum = gmpy.mpf(0)
    
    for i in xrange(n):
        # (-1)i
        a0 = pow(-1, i)
        # (6i)!
        a1 = factorial(6 * i)
        # 13591409 + 545140134i
        a2 = (545140134 * i) + 13591409
        # product of As.
        a = a0 * a1 * a2
        # (3i)!
        b0 = factorial(3 * i)
        # (i!) ^ 3
        b1 = pow(factorial(i), 3)
        # 640320 ^ (3i + 3/2)
        c = gmpy.mpf(640320)
        b2 = gmpy.fsqrt(c ** ((6 * i) + 3))
        # product of Bs
        b = gmpy.mpf(b0) * gmpy.mpf(b1) * b2
        
        div = gmpy.mpf(a) / b
        sum += div
    
    sum *= 12
    
    return gmpy.mpf(1) / sum

def test_calcpi(func, machin_set = None):
    c = 2
    ls = ""
    t = time.clock()
    final_time = 0
    
    print "Bits:"
    n = sys.stdin.readline()
    
    while True:
        ts = time.clock()
        
        if machin_set is None:
            s = str(func(c, int(n)))
        else:
            s = str(func(c, int(n), machin_set))
        
        final_time = time.clock() - ts
        
        if ls == s:
            break
        
        ls = s
        
        sys.stdout.write("\ri: " + str(c) + (" " * (8 - len(str(c)))) + " pi: " + s)
        c += 1
    
    print "\ntime: " + str(time.clock() - t) + " seconds"
    print "iterations (until no more change): " + str(c - 1)
    print "final time: " + str(final_time) + " seconds"

if __name__ == "__main__":
    import sys
    print "  ----  ----\nTesting Polygon method..."
    print "  ---- HISTORY ----\nI used to use math.cos() for this but it sucks so I wrote my own Taylor series based version of cosine. And it finally works (it's very slow)!"
    test_calcpi(calcpi_polygon)
    print "  ----  ----\nTesting Newton method..."
    print "  ---- HISTORY ----\nThis the first method I implemented. Works pretty slowly (but it's way faster than the Polygon method)."
    test_calcpi(calcpi_newton)
    print "  ----  ----\nTesting Madhava method..."
    print "  ---- HISTORY ----\nThis method is so fast (and old)!"
    test_calcpi(calcpi_madhava)
    print "  ----  ----\nTesting Machin (Euler) method..."
    print "  ---- HISTORY ----\nThese Machin methods' number sets were initially hardcoded, but then I decided to use one function and multiple \"sets\"."
    test_calcpi(calcpi_machin, machin_euler_set)
    print "  ----  ----\nTesting Machin (1896) method..."
    test_calcpi(calcpi_machin, machin_1896_set)
    print "  ----  ----\nTesting Machin (1982) method..."
    test_calcpi(calcpi_machin, machin_1982_set)
    print "  ----  ----\nTesting Machin (2003) method..."
    test_calcpi(calcpi_machin, machin_2003_set)
    print "  ----  ----\nTesting Chudnovsky method..."
    print "  ---- HISTORY ----\nI wanted to try this out because it's used by a pi program on the GMP website. It is SO FAST!"
    test_calcpi(calcpi_chudnovsky)
    print "  ----  ----\nTesting Leibniz method..."
    print "  ---- HISTORY ----\nHit Control+C now, this will take forever..."
    test_calcpi(calcpi_leibniz)