Un problema de ciencia computacional.

El día de hoy me crucé con un problema bastante interesante que dice lo siguiente:

Dada una función cuya salida son números aleatorios entre 0 y 1, hallar el valor de pi.

Antes de suponer que ambas proposiciones son totalmente antagónicas, es bueno retomar unos conceptos básicos acerca de pi. Sabemos que el número pi expresa una relación entre el diametro y la longitud de una circunferenci tal como dice la siguiente expresión.

De igual manera, sabemos que entre más lados tenga un poligono éste se parecerá cada vez mas a una circunferencia, por lo que para un cuadrado de lado 1, nuestra fórmula nos dice que

De igual manera usando la apotema como el radio, para un hexagono de lado 1 nos da un valor de 3.45000000001 y para un octógono 3.32.

Sabiendo esto, reformulemos el problema de la siguiente manera:

Teniendo un circulo de radio 1 inscrito en un cuadrado de lado 2, calcula pi.

Tenemos las siguientes dos fórmulas para el área del cuadrado y del circulo

 

 Y el área del cuadrado puede ser reescrita como . Dividiendo ambas áreas

Regresando al problema original, la función generadora de números aleatorios puede ser llamado dos veces, pudiendo así simular puntos en un plano cartesiano y los puntos dentro del circulo pueden ser entendidos como el área del circulo. Con esta interpretación quiero decir que

Dónde Pc es la cantidad de puntos dentro del circulo y Ps la cantidad de puntos dentro del cuadrado, luego                   

Y ésta es la ecuación que usaremos para calcular pi.

Ahora, para determinar la si un punto está dentro del circulo basta con verificar que la distancia del origen al punto sea menor que 1.

 El programa queda de la siguiente manera.

import random

def rnd():

    return random.random()

def solve():

    points = []

    i = 0

    count = 0

    while i < 1000 :

        points.append([rnd(),rnd()])

        i += 1

    f = lambda x, y: (x**2 + y**2)**0.5

    for point in points:

        y = f(point[0], point[1])

        if y <= 1:

            count += 1

    return 4*count/len(points)

print(solve()) 

[Running] python -u ".\Documents\solve.py"

3.196

[Done] exited with code=0 in 0.8 seconds

Entre más grande sea el número de iteraciones, más acertado será el resultado.

Este es el problema original