#!/usr/bin/env python
"""Simple demonstration of Python's arbitrary-precision integers."""

# We need exact division between integers as the default, without manual
# conversion to float b/c we'll be dividing numbers too big to be represented
# in floating point.
from __future__ import division

def pi(n):
    """Compute pi using n terms of Wallis' product.

    Wallis' formula approximates pi as

    pi(n) = 2 \prod_{i=1}^{n}\frac{4i^2}{4i^2-1}."""
    
    num = 1
    den = 1
    for i in xrange(1,n+1):
	tmp = 4*i*i
	num *= tmp
	den *= tmp-1
    return 2.0*(num/den)

if __name__ == '__main__':
    # Simple convergence demo.

    import pylab as P
    import numpy as N

    nrng = N.linspace(10,2000,20).astype(N.int_)
    wpi = N.array(map(pi,nrng))
    diff = N.absolute(wpi-N.pi)
    # plot with lines
    P.semilogy(nrng,diff)
    # and superimpose red circles
    P.semilogy(nrng,diff,'ro')
    P.title(r"Convergence of Wallis' product formula for $\pi$")
    P.xlabel('Number of terms')
    P.ylabel(r'$|\rm{Error}|')
    P.grid()
    P.show()