from matplotlib.pyplot import *
from numpy import *

# MC integration 1
def compute_pi1(N):
    step = 0
    sum = 0

    while step < N:
        x = random.random()
        y = random.random()

        if x*x + y*y < 1.0:
            sum +=  1
        step += 1

    return (4.0 * sum) / N

# MC integration 2
def compute_pi2(N):
    step = 0
    sum = 0.0

    while step < N:
        x = random.random()
        sum += (1.0 - x*x) ** 0.5
        step += 1

    return 4.0 * sum / N

# Numerical integration 3
def compute_pi3(N):
    x = 0.0
    sum = 0.0

    while x < 1.0:
        sum += (1.0 - x*x) ** 0.5
        x += 1.0/N

    return 4.0 * sum / N


N = []
err1s = []
err2s = []
err3s = []

n_trial = 10

for i in range(21):
    n = 2**i
    N.append(n)
    print "i={} N={}".format(i, n)

    err1 = err2 = err3 = err4 = 0
    for j in range(n_trial):
        pi1 = compute_pi1(n)
        err1 += abs(math.pi - pi1)
        
        pi2 = compute_pi2(n)
        err2 += abs(math.pi - pi2)
        
        pi3 = compute_pi3(n)
        err3 += abs(math.pi - pi3)
        
    err1s.append(err1/n_trial)
    err2s.append(err2/n_trial)
    err3s.append(err3/n_trial)

    print "pi1={:e} pi2={:e} pi3={:e}".format(pi1, pi2, pi3)

loglog(N, err1s, 'o-', N, err2s, 'x-', N, err3s, '*-')

show()