%pylab inline
import scipy

def monte_carlo_pi(N):
    """
    Percent of area of square from (0,0) to (1,1) 
    that is inside x^2 + y^2 < 1 is pi/4.
    """
    x = scipy.random.rand(N)
    y = scipy.random.rand(N)
    pi = 4.0*sum(x*x+y*y < 1)/N
    pi_err = abs(pi - scipy.pi)
    return pi, pi_err

def leibniz_pi(N):
    """ 
    pi/4 = 1 - 1/3 + 1/5 - 1/7 + 1/9 - ...
    Obtained by Taylor expansion of arctan  (arctan(x) = x - x^3/3 + x^5/5 - x^7/7 + x^9/9 + ... )
    Also, noting arctan(1) = pi/4
    """
    pi = 4*sum((-1)**i/(2*i+1.0) for i in range(N))
    pi_err = abs(pi - scipy.pi)
    return pi, pi_err

def machin_pi(N):
    """
    pi/4 = 4 arctan(1/5) - arctan(1/239)
    """
    def arctan(x, N):
        return sum((-1)**i *(x**(2*i+1.0))/(2*i+1.0) for i in range(N))
    pi = 16*arctan(1/5., N//2) - 4 * arctan(1/239.0, N//2)
    pi_err = abs(pi - scipy.pi)
    return pi, pi_err

z = scipy.arange(6,100)
z_exp = scipy.floor((1.2)**z)

data = {}
for f in (monte_carlo_pi, leibniz_pi, machin_pi):
    data[f] = []
    for z in z_exp:
        data[f].append(f(int(z)))
    

mc, mc_err = zip(*data[monte_carlo_pi])
leibniz, leibniz_err = zip(*data[leibniz_pi])
machin, machin_err = zip(*data[machin_pi])
plot(z_exp,mc, 'r', label="Monte Carlo")
plot(z_exp, leibniz, 'b', label="Leibniz")
plot(z_exp, machin, 'g', label="Machin")
xscale('log')
axis([1,10**7,0,5])
legend(loc='lower right')

plot(z_exp,mc_err, 'r', label = "MC Error")
plot(z_exp, leibniz_err, 'b', label="Leibniz Error")
plot(z_exp, machin_err, 'g', label="Machin Error")
xscale('log')
yscale('log')
legend(loc='lower right')

import timeit
print "Monte Carlo Time:", timeit.repeat("monte_carlo_pi(100000)", setup="""def monte_carlo_pi(N):
    import scipy
    x = scipy.random.rand(N)
    y = scipy.random.rand(N)
    pi = 4.0*sum(x*x+y*y < 1)/N
    pi_err = abs(pi - scipy.pi)
    return pi, pi_err""", number=10, repeat=3)
print "Leibniz Time", timeit.repeat("leibniz_pi(100000)", setup="""def leibniz_pi(N):
    pi = 4*sum((-1)**i/(2*i+1.0) for i in range(N))
    pi_err = abs(pi - 3.141592653589793238)
    return pi, pi_err""", number=10, repeat=3)
print "Leibniz (no -1^i by expanding out two terms) Time", timeit.repeat("leibniz_pi(100000)", setup="""def leibniz_pi(N):
    pi = 4*sum(1.0/(4*i+1.0) - 1.0/(4*i+3.0) for i in range(N//2))
    pi_err = abs(pi - 3.141592653589793238)
    return pi, pi_err""", number=10, repeat=3)
print "Machin Time", timeit.repeat("machin_pi(100000)", setup="""def machin_pi(N):
    def arctan(x, N):
        return sum((-1)**i *(x**(2*i+1.0))/(2*i+1.0) for i in range(N))
    pi = 16*arctan(1/5., N//2) - 4 * arctan(1/239.0, N//2)
    pi_err = abs(pi - 3.141592653589793238)
    return pi, pi_err""", number=10, repeat=3)
print "Machin (N=22) Time", timeit.repeat("machin_pi(22)", setup="""def machin_pi(N):
    def arctan(x, N):
        return sum((-1)**i *(x**(2*i+1.0))/(2*i+1.0) for i in range(N))
    pi = 16*arctan(1/5., N//2) - 4 * arctan(1/239.0, N//2)
    pi_err = abs(pi - 3.141592653589793238)
    return pi, pi_err""", number=10, repeat=3)

print "Machin(22):", machin_pi(22)
print "Leibniz(100000):", leibniz_pi(100000)

from mpmath import mp
def machin_mp_pi(N):
    """
    pi/4 = 4 arctan(1/5) - arctan(1/239)
    """
    mp.dps = N if N > 25 else 25
    def arctan(x, N):
        return sum((x**(4*i+1.0))/(4*i+1.0) - (x**(4*i+3.0))/(4*i+3.0) for i in range(N//2))
    pi = 16*arctan(mp.mpf("1/5"), N//2) - 4 * arctan(mp.mpf("1/239"), N//2)
    pi_err = float(abs(pi - mp.pi))
    return pi, pi_err

for d in [10,20,30,40,50,60,70,80,90,100,200,300,400]:
    print "N: %d, Error: %g" % (d, machin_mp_pi(d)[1])