Varios métodos para calcular PI

Serie armónica

def forma4m1(n):
        'da +1 ó -1 según un número sea de la forma 4m+1 o 4m-1'
	if (n+1)%4 == 0:
		return +1
	if (n-1)%4 == 0:
		return -1
	return +1

# Aquí guardaremos a lista de valores de pi por cada término
api=[]
# Los dos primeros términos son 1 y 1/2
api.append([1,1])
api.append([2,1.5])
mipi=1+0.5	

terminos=10000

for i in range(terminos)[3:]:
	fa=list(factor(i))
	if is_prime(i):
		signo=forma4m1(i)
	else:
		signo=1
		for f in fa:
                        for f2 in range(f[1]):
			    signo*=forma4m1(f[0])
	mipi+=signo*1./i
	api.append([i,mipi])

print mipi
list_plot(api, plotjoined=True)

3.11042957448254

Montecarlo

2 dimensiones

npuntos=10000
dentro=0.
fuera=0.
api=[]
for i in range(npuntos):
    x=random()
    y=random()
    r=x**2+y**2
    if r < 1:
        dentro+=1
    else:
        fuera+=1
    pi=4*dentro/(fuera+dentro)
    api.append(pi)
print float(pi)
list_plot(api,plotjoined=True)

3.166

3 dimensiones

npuntos=10000
dentro=0.
fuera=0.
api=[]
for i in range(npuntos):
    x=random()
    y=random()
    z=random()
    r=x**2+y**2+z**2
    if r < 1:
        dentro+=1
    else:
        fuera+=1
    mipi=6*dentro/(fuera+dentro)
    api.append(mipi)
print mipi
list_plot(api,plotjoined=True)

3.11220000000000

Pseudo-montecarlo (grid)

npuntos=10000
pts=sqrt(npuntos)
dentro=0.
fuera=0.
api=[]
for x in range(pts):
    for y in range(pts):
        r=(x/pts)**2+(y/pts)**2
        if r < 1:
            dentro+=1
        else:
            fuera+=1
    mipi=4*dentro/(fuera+dentro)
    api.append(mipi)
print float(mipi)
list_plot(api,plotjoined=True)

3.1796