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ódulomath. Podemos demostrar esto usando el shell de Python:>>> from math import pi >>> pi 3.141592653589793Podemos 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.1415929203539825Podemos 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 * sPodemos 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.1415925535897915Se 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} \).
 |
| 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 / tLlamando 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.1415684Debido 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 tipofloat 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 xProbando 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.
Hola! muchas gracias, estoy aprendiendo a programar en python y esto me pareció fabuloso!
ResponderBorrarGracias mdasuaje, me alegra saber que te gustó.
Borrar¡Feliz día de Pi, profesor Ariel!
ResponderBorrarMuchas gracias, Allie. Te mando un fuerte abrazo.
BorrarExcelente aporte!
ResponderBorrar¡Gracias!
BorrarDe gran utilidad, Gracias!
ResponderBorrarGracias. Tu nota me resulto soberbia. Voy a afinar con tu ayuda mi programación en phyton.
ResponderBorrarCuando copeo el programa tal y cual al python tommy no corre falla en lan fila 1 ya comenzando con ese error. O como le hago paaa que corra el programa?
ResponderBorrargrandee
ResponderBorrar