14 de marzo de 2015

Pensando en Pi

Esta entrada del blog de EduPython es para conmemorar el Día de Pi del 2015. Como sabemos, π (pi) es una constante matemática que representa la relación que tiene la circunferencia de un círculo con respecto a su diámetro. Esta constante es un número irracional, lo cual significa que posee infinitas cifras decimales sin patrón de repetición alguno. Sus primeras 20 cifras son: 3.1415926535897932384


A partir del formato de fechas que se utiliza en los Estados Unidos (mes/día), el Día de Pi se celebra el 3/14 (14 de marzo) de cada año, dado que 3, 1 y 4 son los tres primeros dígitos de π. Como feliz coincidencia, el 14 de marzo también fue el día que nació Albert Einstein.

Albert Einstein nació el 14 de marzo
de 1879 en Ulm, Alemania.

Cada año el Día de Pi§ es celebrado por estudiantes, profesores, matemáticos y computólogos alrededor de todo el mundo. Las celebraciones incluyen actividades relacionadas con π, por ejemplo: competencias que consisten en recitar de memoria la mayor cantidad de dígitos de π, hornear y comer pays, o incluso hasta escribir blogs acerca de π.

Un rico pay de moras para celebrar el Día de Pi.

El año 2015 es especial, ya que una sola vez en cada siglo ocurre una combinación de fecha y hora que incluye los primeros 10 dígitos de π: 3/14/15 a las 9:26:53 hrs.

Poniéndonos más técnicos, veremos ahora algunas formas de obtener valores aproximados de π usando Python. Todos los ejemplos de código que se muestran a continuación fueron probados en Python versión 3.4.

π en Python

Cuando un programa en Python requiere del valor de π, lo más recomendable es importarlo directamente del módulo math. Podemos demostrar esto usando el shell de Python:
>>> from math import pi
>>> pi
3.141592653589793
Podemos observar que el valor de pi tiene 16 dígitos de precisión debido a que ésa es la cantidad máxima de dígitos decimales que puede tener un valor de tipo float en una implementación típica de Python. Dicha aproximación es adecuada para la mayoría de los cálculos que comúnmente se requieren en las diversas disciplinas científicas e ingenieriles.

A continuación exploraremos diversas formas programáticas de calcular el valor de π.

π como un número racional

En tiempos pasados era común usar fracciones como 22/7 y 355/113 para aproximar π.  Evaluando dichas fracciones:
>>> 22/7
3.142857142857143
>>> 355/113
3.1415929203539825
Podemos notar que el resultado de la primera división tiene tres dígitos correctos, mientras que la segunda división acierta siete.

Serie de Gregory–Leibniz

En los siglos XVII y XVIII, James Gregory y Gottfried Leibniz descubrieron una serie infinita que sirve para calcular π: $$ \begin{align*} \pi &= 4 \left ( \sum_{k=1}^{\infty} \frac{(-1)^{(k + 1)}}{2k-1} \right ) \\ &= 4 \left ( 1 - \frac{1}{3} + \frac{1}{5} - \frac{1}{7} + \frac{1}{9} - \frac{1}{11} \cdots \right ) \end{align*} $$ Resulta relativamente fácil traducir esta fórmula a una función en Python que nos permita aproximar el valor de π tomando los n primeros términos de la serie:
def gregory_leibniz(n):
    """Calcula y devuelve el valor de pi usando
    los primeros n términos de la serie de
    Gregory–Leibniz.
    
    π = 4(1 - 1/3 + 1/5 - 1/7 + ...)
    """
    s = 0
    for k in range(1, n + 1):
        s += (-1)**(k + 1) / (2 * k - 1)
    return 4 * s
Podemos probar desde el shell nuestra función con diferentes valores de n:
>>> gregory_leibniz(1)
4.0
>>> gregory_leibniz(10)
3.0418396189294032
>>> gregory_leibniz(100)
3.1315929035585537
>>> gregory_leibniz(1000)
3.140592653839794
>>> gregory_leibniz(10000)
3.1414926535900345
>>> gregory_leibniz(100000)
3.1415826535897198
>>> gregory_leibniz(1000000)
3.1415916535897743
>>> gregory_leibniz(10000000)
3.1415925535897915
Se puede observar que se obtiene un dígito más de precisión cada vez que multiplicamos por 10 el argumento de la función gregory_leibniz(). Sin embargo, tal como es de esperarse, entre más grande sea el valor de n más tiempo tarda la función en completar.

Método de Montecarlo

Supongamos que tenemos un tablero para jugar a los dardos. Dicho tablero es un cuadrado de un metro de cada lado. El tablero contiene un cuarto de círculo tal como se muestra en la siguiente imagen:


Ahora lanzamos t dardos al tablero. Si todos los dardos caen de manera aleatoria dentro del tablero con una distribución uniforme, algunos dardos caerán en la región gris y otros en la región blanca. Los puntos en la siguiente imagen representan los lugares donde pudieron haber caído los t dardos:


Sabemos que el área del tablero es: 1 metro × 1 metro = 1 metro2. El área de la región gris es un cuarto del área de un círculo, es decir: \((\pi \times r^2) \div 4\). Dado que el radio r mide 1 metro (lo que mide un lado del tablero), entonces el área de la región gris es: $$ \frac{\pi \times (1\;\textrm{m}^2)}{4} = \frac{\pi}{4} \textrm{m}^2 $$ Si t es el total de dardos lanzados y g es la cantidad de esos dardos que cayeron dentro de la región gris, entonces podemos considerar que \(g/t\) debe ser aproximadamente igual a la división del área del cuarto de círculo (π/4 metros2) entre el área de todo el tablero (1 metro2): $$ \frac{\frac{\pi}{4} \textrm{m}^2}{1 \; \textrm{m}^2} \approx \frac{g}{t} $$ $$ \frac{\pi}{4} \approx \frac{g}{t} $$ $$ \pi \approx \frac{4 g}{t} $$ Haciendo los despejos necesarios, obtenemos otra forma de aproximar el valor de π. El método de Montecarlo para calcular π quedaría así:
  • Inicializar g en 0.
  • Repetir t veces lo siguiente:
    • Generar de manera aleatoria un punto con coordenadas (x, y) dentro del área del tablero. Dado que que cada lado del tablero mide un metro, los valores de x y y deben ser números reales entre 0 y 1.
    • Calcular la distancia d entre el centro del círculo (la esquina inferior izquierda del tablero) y el punto (x, y). Para ello utilizamos el teorema de Pitágoras: \(d = \sqrt{x^2 + y^2}\)
    • Si d es menor a un metro (el radio del círculo) entonces el punto (x, y) está dentro del área del círculo. En ese caso incrementamos en uno el valor de g.
  • El valor aproximado de π es: \( \frac{4 g}{t} \).
Se le llama Montecarlo a este método en referencia al casino que se encuentra en Mónaco, el cual es considerado por muchos como la capital mundial de los juegos de azar.

Casino de Montecarlo en el Principado
de Mónaco.

La implementación del método de Montecarlo en Python es bastante directa. Para generar los números aleatorios usamos la función random() del módulo random. Dicha función devuelve un número de punto flotante al azar en el intervalo [0, 1), es decir, dicho número es mayor o igual a cero pero menor a uno. Y esto es justo lo que nuestro algoritmo necesita. El código quedaría así:
from random import random
from math import sqrt

def montecarlo(t):
    """Calcula y devuelve el valor aproximado
    de pi usando el método de Montecarlo a
    partir de t puntos.
    """
    g = 0
    for i in range(t):
        x = random()
        y = random()
        d = sqrt(x ** 2 + y ** 2)
        if d < 1:
            g += 1
    return 4 * g / t
Llamando la función con t = 100,000,000 (cien millones) podemos obtener casi cinco dígitos de precisión:
>>> montecarlo(100000000)
3.14166808
>>> montecarlo(100000000)
3.14167544
>>> montecarlo(100000000)
3.1415684
Debido al uso de números aleatorios, cada invocación a la función montecarlo() produce resultados (ligeramente) diferentes.

Calculando muchos dígitos de π

Un problema que tienen todos los esquemas discutidos anteriormente para el cálculo de π es que están limitados a la representación interna del tipo float de Python. La mayoría de las implementaciones de Python usan números de punto flotante de precisión doble de 64 bits tal como lo describe el estándar IEEE 754 (también conocido como IEC 60559). Para fines prácticos, este tipo de dato permite representar valores numéricos de 15 a 17 dígitos decimales significativos. Pero, ¿qué pasa si queremos calcular más dígitos de π? Para eso tenemos los algoritmos de espita (spigot en inglés) que únicamente requieren números enteros. Así como gotea el agua de una espita (grifo, válvula o llave), a los algoritmos de espita se les llama así debido a que producen dígitos individuales de π que no se reusan después de que son calculados. Esto contrasta con las series infinitas o algoritmos iterativos, en los que se retienen y utilizan todos los dígitos intermedios hasta que se produce el resultado final.

Una espita
(grifo, válvula o llave)

El primer algoritmo de espita se le atribuye a A.H.J. Sale, que en 1968 presentó una manera de calcular muchos dígitos de e. En 1995 Stan Wagon y Stanley Rabinowitz publicaron un algoritmo de espita que permite calcular una cantidad arbitraria de dígitos de π.

Boris Gourévitch hace un buen trabajo explicando el algoritmo de espita para calcular π. Dado que el tema es algo complicado, no entraré en más detalles aquí. Solo me limitaré a mostrar un programa escrito por John Zelle en el 2006 que implementa el algoritmo de espita publicado un año antes por Jeremy Gibbons que a su vez mejoró el trabajo original de Wagon y Rabinowitz.
def espita(d):
    """Regresa una lista con los primeros d dígitos
       de pi utilizando el algoritmo de espita
       diseñado por Jeremy Gibbons. Implementación
       de John Zelle con ligeras alteraciones por
       Ariel Ortiz.
    """
    x = []
    q,r,t,k,n,l = 1,0,1,1,3,3
    while len(x) < d:
        if 4*q+r-t < n*t:
            x.append(n)
            q,r,t,k,n,l = (
                10*q,10*(r-n*t),t,k,
                (10*(3*q+r))//t-10*n,l)
        else:
            q,r,t,k,n,l = (
                q*k,(2*q+r)*l,t*l,k+1,
                (q*(7*k+2)+r*l)//(t*l),l+2)
    return x
Probando la función:
>>> espita(1)
[3]
>>> espita(3)
[3, 1, 4]
>>> espita(5)
[3, 1, 4, 1, 5]
>>> espita(10)
[3, 1, 4, 1, 5, 9, 2, 6, 5, 3]
>>> espita(20)
[3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 
 3, 2, 3, 8, 4]
>>> espita(100)
[3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 
 3, 2, 3, 8, 4, 6, 2, 6, 4, 3, 3, 8, 3, 2, 7, 
 9, 5, 0, 2, 8, 8, 4, 1, 9, 7, 1, 6, 9, 3, 9, 
 9, 3, 7, 5, 1, 0, 5, 8, 2, 0, 9, 7, 4, 9, 4, 
 4, 5, 9, 2, 3, 0, 7, 8, 1, 6, 4, 0, 6, 2, 8, 
 6, 2, 0, 8, 9, 9, 8, 6, 2, 8, 0, 3, 4, 8, 2, 
 5, 3, 4, 2, 1, 1, 7, 0, 6, 7]
>>> espita(1000)
[3, 1, 4, 1, 5, 9, 2, 6, 5, 3, 5, 8, 9, 7, 9, 
 3, 2, 3, 8, 4, 6, 2, 6, 4, 3, 3, 8, 3, 2, 7, 
 9, 5, 0, 2, 8, 8, 4, 1, 9, 7, 1, 6, 9, 3, 9, 
 9, 3, 7, 5, 1, 0, 5, 8, 2, 0, 9, 7, 4, 9, 4, 
 4, 5, 9, 2, 3, 0, 7, 8, 1, 6, 4, 0, 6, 2, 8, 
 6, 2, 0, 8, 9, 9, 8, 6, 2, 8, 0, 3, 4, 8, 2, 
 5, 3, 4, 2, 1, 1, 7, 0, 6, 7, 9, 8, 2, 1, 4, 
 8, 0, 8, 6, 5, 1, 3, 2, 8, 2, 3, 0, 6, 6, 4, 
 7, 0, 9, 3, 8, 4, 4, 6, 0, 9, 5, 5, 0, 5, 8, 
 2, 2, 3, 1, 7, 2, 5, 3, 5, 9, 4, 0, 8, 1, 2, 
 8, 4, 8, 1, 1, 1, 7, 4, 5, 0, 2, 8, 4, 1, 0, 
 2, 7, 0, 1, 9, 3, 8, 5, 2, 1, 1, 0, 5, 5, 5, 
 9, 6, 4, 4, 6, 2, 2, 9, 4, 8, 9, 5, 4, 9, 3, 
 0, 3, 8, 1, 9, 6, 4, 4, 2, 8, 8, 1, 0, 9, 7, 
 5, 6, 6, 5, 9, 3, 3, 4, 4, 6, 1, 2, 8, 4, 7, 
 5, 6, 4, 8, 2, 3, 3, 7, 8, 6, 7, 8, 3, 1, 6, 
 5, 2, 7, 1, 2, 0, 1, 9, 0, 9, 1, 4, 5, 6, 4, 
 8, 5, 6, 6, 9, 2, 3, 4, 6, 0, 3, 4, 8, 6, 1, 
 0, 4, 5, 4, 3, 2, 6, 6, 4, 8, 2, 1, 3, 3, 9,
 3, 6, 0, 7, 2, 6, 0, 2, 4, 9, 1, 4, 1, 2, 7,
 3, 7, 2, 4, 5, 8, 7, 0, 0, 6, 6, 0, 6, 3, 1,
 5, 5, 8, 8, 1, 7, 4, 8, 8, 1, 5, 2, 0, 9, 2,
 0, 9, 6, 2, 8, 2, 9, 2, 5, 4, 0, 9, 1, 7, 1,
 5, 3, 6, 4, 3, 6, 7, 8, 9, 2, 5, 9, 0, 3, 6,
 0, 0, 1, 1, 3, 3, 0, 5, 3, 0, 5, 4, 8, 8, 2,
 0, 4, 6, 6, 5, 2, 1, 3, 8, 4, 1, 4, 6, 9, 5,
 1, 9, 4, 1, 5, 1, 1, 6, 0, 9, 4, 3, 3, 0, 5,
 7, 2, 7, 0, 3, 6, 5, 7, 5, 9, 5, 9, 1, 9, 5,
 3, 0, 9, 2, 1, 8, 6, 1, 1, 7, 3, 8, 1, 9, 3,
 2, 6, 1, 1, 7, 9, 3, 1, 0, 5, 1, 1, 8, 5, 4,
 8, 0, 7, 4, 4, 6, 2, 3, 7, 9, 9, 6, 2, 7, 4,
 9, 5, 6, 7, 3, 5, 1, 8, 8, 5, 7, 5, 2, 7, 2,
 4, 8, 9, 1, 2, 2, 7, 9, 3, 8, 1, 8, 3, 0, 1,
 1, 9, 4, 9, 1, 2, 9, 8, 3, 3, 6, 7, 3, 3, 6,
 2, 4, 4, 0, 6, 5, 6, 6, 4, 3, 0, 8, 6, 0, 2,
 1, 3, 9, 4, 9, 4, 6, 3, 9, 5, 2, 2, 4, 7, 3,
 7, 1, 9, 0, 7, 0, 2, 1, 7, 9, 8, 6, 0, 9, 4,
 3, 7, 0, 2, 7, 7, 0, 5, 3, 9, 2, 1, 7, 1, 7,
 6, 2, 9, 3, 1, 7, 6, 7, 5, 2, 3, 8, 4, 6, 7,
 4, 8, 1, 8, 4, 6, 7, 6, 6, 9, 4, 0, 5, 1, 3,
 2, 0, 0, 0, 5, 6, 8, 1, 2, 7, 1, 4, 5, 2, 6,
 3, 5, 6, 0, 8, 2, 7, 7, 8, 5, 7, 7, 1, 3, 4,
 2, 7, 5, 7, 7, 8, 9, 6, 0, 9, 1, 7, 3, 6, 3,
 7, 1, 7, 8, 7, 2, 1, 4, 6, 8, 4, 4, 0, 9, 0,
 1, 2, 2, 4, 9, 5, 3, 4, 3, 0, 1, 4, 6, 5, 4,
 9, 5, 8, 5, 3, 7, 1, 0, 5, 0, 7, 9, 2, 2, 7,
 9, 6, 8, 9, 2, 5, 8, 9, 2, 3, 5, 4, 2, 0, 1,
 9, 9, 5, 6, 1, 1, 2, 1, 2, 9, 0, 2, 1, 9, 6,
 0, 8, 6, 4, 0, 3, 4, 4, 1, 8, 1, 5, 9, 8, 1,
 3, 6, 2, 9, 7, 7, 4, 7, 7, 1, 3, 0, 9, 9, 6,
 0, 5, 1, 8, 7, 0, 7, 2, 1, 1, 3, 4, 9, 9, 9,
 9, 9, 9, 8, 3, 7, 2, 9, 7, 8, 0, 4, 9, 9, 5,
 1, 0, 5, 9, 7, 3, 1, 7, 3, 2, 8, 1, 6, 0, 9,
 6, 3, 1, 8, 5, 9, 5, 0, 2, 4, 4, 5, 9, 4, 5,
 5, 3, 4, 6, 9, 0, 8, 3, 0, 2, 6, 4, 2, 5, 2,
 2, 3, 0, 8, 2, 5, 3, 3, 4, 4, 6, 8, 5, 0, 3,
 5, 2, 6, 1, 9, 3, 1, 1, 8, 8, 1, 7, 1, 0, 1,
 0, 0, 0, 3, 1, 3, 7, 8, 3, 8, 7, 5, 2, 8, 8,
 6, 5, 8, 7, 5, 3, 3, 2, 0, 8, 3, 8, 1, 4, 2,
 0, 6, 1, 7, 1, 7, 7, 6, 6, 9, 1, 4, 7, 3, 0,
 3, 5, 9, 8, 2, 5, 3, 4, 9, 0, 4, 2, 8, 7, 5,
 5, 4, 6, 8, 7, 3, 1, 1, 5, 9, 5, 6, 2, 8, 6,
 3, 8, 8, 2, 3, 5, 3, 7, 8, 7, 5, 9, 3, 7, 5,
 1, 9, 5, 7, 7, 8, 1, 8, 5, 7, 7, 8, 0, 5, 3,
 2, 1, 7, 1, 2, 2, 6, 8, 0, 6, 6, 1, 3, 0, 0,
 1, 9, 2, 7, 8, 7, 6, 6, 1, 1, 1, 9, 5, 9, 0,
 9, 2, 1, 6, 4, 2, 0, 1, 9, 8]
Al momento de estar escribiendo esta entrada, el mayor número de dígitos que se han calculado de π son 12.1 billones (1.21 × 1013) por Alexander Yee y Shigeru Kondo el 28 de diciembre del 2013.  Aplicando el algoritmo de Chudnovsky tardaron 94 días en total usando un equipo de cómputo con dos procesadores Intel Xeon E5-2690 a 2.9 GHz con 16 núcleos en total, 128 GB de memoria RAM y más de 60 TB de espacio en varios discos duros. Eso sí es como para pensar en π.

¡Feliz día de Pi!


Notas

§ El Día de Pi fue fundado por Larry Shaw y se celebró por primera vez en 1988 en el Exploratorium de San Francisco, California.
Referencia: http://en.wikipedia.org/wiki/Pi_Day.

En México pay es la castellanización de la palabra inglesa pie. En España se le llama tarta, mientras que en otros países de América Latina se le conoce como pastel. A fin de cuentas la broma consiste en que en inglés pie y π se pronuncian igual.

Realmente la serie Gregory–Leibniz fue descubierta tres siglos antes por el matemático indio Madhava de Sangamagrama. Por eso también se le conoce a esta fórmula como la serie de Madhava-Leibniz.

10 comentarios: