Algoritmos para el cálculo de Pi
En toda la historia se han descubierto/creado muchas maneras para calcular aproximaciones del número Pi. En el post del día de Pi ya vimos algunas aproximaciones numéricas. Hoy os traigo dos algoritmos para el cálculo de aproximaciones de Pi. Vamos con ellos:
- Algoritmo de Gauss-Legendre:
Partimos de los siguientes datos iniciales:

A partir de ellos realizamos la siguientes operaciones:

Entonces Pi se aproxima de la siguiente forma:

El método tiene convergencia de segundo orden. Es decir, en cada iteración duplicamos el número de dígitos exactos obtenidos en la iteración anterior.
- Algoritmo de Borwein:
Partimos de los siguientes datos iniciales:

Con ellos operamos de la siguiente forma:

Entonces se tiene lo siguiente:

La convergencia de este método es cuártica. Es decir, en cada iteración se consiguen el cuádruple de dígitos exactos que en la iteración anterior. Existen variaciones de este método que consiguen en cada iteración muchos más dígitos exactos. En este enlace de la Wikipedia inglesa podéis ver algunas de ellas.
Evidentemente existen muchísimos más. Y también muchas fórmulas que involucran a Pi mediante las cuales podemos calcular aproximaciones de este número. Más adelante irán apareciendo en este blog muchas más. Se aceptan sugerencias en los comentarios.


Joan | 11 de Junio de 2007 | 10:49
Particularmente me gusta la idea de la aproximación por el método de Montecarlo:
Construimos un cuadrado de lado igual a 2 y inscribimos una circumferencia. Entonces la circumferencia tendra un radio igual a 1.
La superfície de la circumferencia serà: pi*r^2 = pi.
Y la superfície del cuadrado: l^2 = 4.
Entonces, si empezásemos a lanzar dardos aleatoriamente dentro del cuadrado, la probabilidad de que este entrase en el círculo seria de pi/4. Solo deberiamos contar las veces que hemos lanzado dardos y cuantos de estos han entrado en el círculo, dividirlo y multiplicarlo por cuatro de forma que:
pi = (dardos en círculo / dardos totales) * 4.
No es un método muy exacto (se necesitarian muchísimos tiros para aproximarse a pi), pero es bonito y ingenioso.
Hice un pequeño programa de simulación en Delphi. Si a alguien le interesa:
http://gftpklient.sourceforge.net/l2p/?p=15
El artículo está en catalán, pero el código fuente es el código fuente, jeje y abajo podéis encontrar el .pas y un .zip para descargarlo y ejectuarlo. Saludos!
Uno de tantos | 11 de Junio de 2007 | 11:42
Mmmm… no me acaban de convencer estos métodos.
Con el primero tendríamos, de entrada, saber el valor de raíz de 2, y luego calcular unas cuantas raices cuadradas en cada valor. No es que me cree un conflicto en mi interior el usar raíces cuadradas, es que tenemos que tener en cuenta que hay que cortar en algún lugar el cálculo de la raíz y eso afecta (y mucho) a lo precisos que queramos ser con Pi
Con el segundo tenemos raíces de orden cuarto… Pues más de lo mismo.
Con lo del método de lanzar dardos sobre una diana, hace ya algunos años hice la simulación y la verdad que saqué 6 decimales de pi con 100.000 tiradas, pero el problema fundamental es la serie aleatoria que uses.
De los 15 generados de números pseudoaleatorios que usé, sólo con uno conseguí converger a pi con más de 3 decimales. Y en teoría eran buenos generadores sacados de papers (uno incluso ofrecía una recompensa de 10.000 $ a quien lograse demostrale una correlación en su generador). Estuve tentado (MUY tentado) de mandarle la simulación, pero dejé correr un tupido velo de silencio.
Con el único con el que conseguí la convergencia de 6 decimales fué con el Lineal Congruencial que usa el Mathemática, pero era especialmente costoso, porque trabaja con un módulo de (2^31)^48 - (2^31)^8 + 1, lo cual es una auténtica burrada computacionalmente hablando.
Y aún así con este generador pasan 2 cositas:
i) No es perfecto (A. M. Ferrenberg, D.P. Landau, and Y.J. Wong. Monte carlo simulations: hidden errors from “good” random number generators. Phys. Rev. Lett., 69:3382-3384, 1992.)
ii) Todo generador de números pseudoaleatorios tiene un periodo con el que acaba volviendo a repetir la serie. Así que tenemos un límite al número de decimales que podríamos conseguir. Podríamos decir que: Mayor número de decimales = Mayor número de tiradas, pero las tiradas están limitadas por el periodo de la serie pseudoaleatoria
La única forma de conseguir sacar Pi con muchísimos decimales y sin error es conectar un contador geiger a un ordenador y ponerlo a medir la radiación de una muestra de algún material radioactivo.
El problema de eso es sería muchiiiiiisimo más lento que cualquier generador de números pseudoaleatorios porque la CPU tendría que estar esperando a que el contador diese medida, y la CPU es mucho más rápida (de orden del millón de veces más) de lo que el contador devuelve datos.
La solución, capturar los datos de contador, ponerlos en un fichero y usar esos datos ya tomados de antemano.
Alfonso Jiménez | 11 de Junio de 2007 | 12:23
Otro puede ser la aguja de Buffon: http://www.genciencia.com/2006/08/28-el-problema-de-la-aguja-de-buffon
Saludos!
wallace | 11 de Junio de 2007 | 13:19
¿¿Cómo se demuestra matemáticamente que esos algoritmos conducen efectivamente a PI al iterar una infinidad de veces??
Trackback | 11 Jun, 2007
meneame.net
mimetist | 11 de Junio de 2007 | 13:47
Más información sobre cómo calcular PI por el “método de Buffon” aquí.
Ondia, eso lo he escrito yo!! qué autobombo más indiscreto xD
kalceto | 11 de Junio de 2007 | 14:38
Un dia descubrí uno muy curioso, ineficiente como él sólo, pero curioso como ninguno. Es usando el Conjunto de Mandelbrot.
Se trata de que en la unión del cardiode con la circunferencia de la izquierda, realmente no existe. La gracias es que mientras más te acercas al punto exacto de la unión (-3/4) en el eje y, más iteraciones se necesitan para demostrar que está fuera del Conjunto de Mandelbrot. Si coges un número muy cercano del estilo -3/4 + εi, donde ε es un número muy pequeño, 0.00000001, p.e. Calculas el número de iteraciones en ese punto necesarias para determinar que el punto no pertenece al conjunto, y lo multiplicas por ε, te da pi, mientras más pequeño sea ε, más te aproximarás a pi.
Como he dicho el método es ineficiente de narices, pero es curioso, usando únicamente multiplicaciones y sumas. Le calculo una complejidad aproximada de 10 elevado al número de decimales que quieras.
ReiVaX18 | 11 de Junio de 2007 | 15:36
En el último número de la Gaceta de la RSME (vol 10 núm 1) hay un artículo muy interesante que describe muchos algoritmos y fórmulas para calcular pi.
Fernando J. Pereda | 11 de Junio de 2007 | 17:29
Hace tiempo lo comenté, y me parece curiosísimo (y cuya aplicación en la vida real es nula claro):
Utilizando la función iteración que genera el conjunto de mandelbrot:
z0 = 0
z = z^2 + c
c = -3/4 + a i
Para a = 0 el punto pertenece al conjunto; sin embargo si n es el número de iteraciones que tardamos en conseguir un z cuyo módulo sea mayor que 2 (radio de escape):
n * a = pi
Y la precisión es la misma que la de a
Es matemáticamente irrelevante quizá, costoso de calcular, pero curioso sin duda
- ferdy
Fernando J. Pereda | 11 de Junio de 2007 | 20:05
Para muestra, un botón: http://dev.gentoo.org/~ferdy/stuff/pi_mandel_gmp.c
- ferdy
witilongi | 11 de Junio de 2007 | 20:52
Os recomiendo el programa gratuito que circula por ahí PiFast para calcular mega archivos con cifras de pi hasta millones y millones en poco rato.
Daniel | 11 de Junio de 2007 | 21:58
Alguna justificación de estos algoritmos?
Roberto | 11 de Junio de 2007 | 21:58
Es interesante la pregunta de wallace.
Las convergencia de las iteraciones tiene su propia demostración en cada caso. Habitualmente, se basa en que el punto límite de la sucesión cumple un papel determinado que sólo puede cumplir el número pi.
Ahora mismo no recuerdo ninguna en concreto, pero síme viene a la cabeza una demostración de convergencia que se limitaba a comprobar que la tangente (en radianes) de la cuarta parte del límite era 1. Para obtener la tangente empleaba una transformada que simplificaba las fórmulas extraordinariamente.
Rober | 12 de Junio de 2007 | 1:56
Acabo de rehacer un método basado en el que ya pensó Arquímedes: el polígono inscrito.
Sabiendo que media circunferencia es pi/2, partimos del lado AB de un cuadrado inscrito: A:(1,0) B:(x=0,y=1), de longitud raiz(2) y aproximamos pi=2*raiz(2). Doblamos los lados del polígono calculando el corte de la recta que pasa por (0,0) y el punto medio del anterior: (Xm=(1+0)/2,Ym=(0+1)/2). El número de lados es 2^i, así que aproximamos pi=2^i*L.
En ecuaciones (ojo con i+1 que es el siguiente término):
i=iteración
—
Li=raiz(Xi^2+Yi^2)
~pi=2^i*Li
Xmi=(1+Xi)/2
Ymi=Yi/2
— para el siguiente término:
Xi+1=raiz(1/(1+Ymi/Xmi))
Yi+1=(Ymi/Xmi)*Xi+1
Li+1=raiz(Xi+1^2+Yi+1^2)
~pi=2^i+1*Li+1
Estos resultados me da, con OpenOffice Calc:
i=1
X1=0
Y1=1
~pi=2,8284271247461900
…
i=12
X12=0,9999997058628820
Y12=0,0007669903187427
~pi=3,1415925765848700
…
i=24
X24=0,9999999999999820
Y24=0,0000001872535141
~pi=3,14159265358979
y ahí se acaba la precisión del OpenOffice, en 14 decimales.
Según mis cuentas empieza sacando 0,5 decimales de pi por cada iteración y va aumentado hasta casi 0,6 decimales, pero la convergenica tiene la pinta (graficándola) de tener un límite o, al menos, de aumentar muy despacio.
Harina de otro costal es lidiar con la inexactitud arrastrada entre iteraciones al calcular raices, como bien dice Unodetantos. No obstante, parece que no “desbarra” y tiene como límite los propios decimales de exactitud con que se realice el cálculo. Pero esto habría que demostrarlo.
Bueno, no está mal para haberlo hecho en un ratico: un método con convergencia de orden “casi 6/10″.
leandro | 12 de Junio de 2007 | 22:30
estoy de acuerdo con witilongi, ya que he usado ese programa y he calculado 1.000.000 de digitos del numero pi. se podrian calcular mas pero el tiempo seria mucho mayor. se los recomiendo ese programa. ademas se puede calcular la raiz de 2 y el numero e.
leandro.
Fermi0n | 14 de Junio de 2007 | 12:24
No entiendo el porqué de tanta pérdida de tiempo. Todo el mundo sabe que pi es exactamente igual a 3
http://buscandoasusy.wordpress.com/2007/06/14/%cf%80-igual-a-3/
Por cierto, estoy de acuerdo con Uno de Tantos.
telaro | 15 de Junio de 2007 | 4:14
1. Toma una calculadora
2. Presiona los tres primeros dígitos de tu numero telefónico
3. Multiplícalo por 80.
4. Súmale 1.
5. Multiplícalo por 250.
6. Súmale los últimos cuatro dígitos de tu numero telefónico.
7. Otra vez, súmale los últimos cuatro dígitos de tu numero telefónico.
8. Réstale 250.
9. Divídelo por 2
Reconoces el numero?????? ????????
daniel | 15 de Junio de 2007 | 20:02
¿Qué tan eficiente y exacto sería usar lo del problema de Basilea para calcular Pi?
(Pi^2)/6 = 1/(1^2) + 1/(2^2) + 1/(3^2)+…+ 1/(n^2)
ó
Pi = Sqrt(6[1/(1^2) + 1/(2^2) + 1/(3^2)+...+ 1/(n^2)])
Adder | 16 de Junio de 2007 | 22:19
daniel, no se de exactitud, yo creo que se dice que un método es mejor que otro en función de la velocidad de convergencia.
manji | 17 de Junio de 2007 | 1:57
y que tal buscar la solución a la ecuación
sen(x) = 0 en [3,4]
que se puede resolver por muchos métodos iterativos (bisección*, newton-raphson…)
o se puede buscar el punto fijo de la función f:
f(x) = sen(x) - x
se empieza, en 3, por ejemplo y seguimos por f(3), f(f(3)), f(f(f(3)))…
*bisección basicamente consiste en dividir [3,4] en dos subintervalos y quedarnos con aquel que cumpla bolzano (sen(a)*sen(b)
Rober | 17 de Junio de 2007 | 14:14
Daniel, he calculado lo del probelma de Basilea y me sale bastante ineficiente: al inicio consigue unas 0,2 cifras por iteración, pero al cabo de 200 baja hasta 0,002. Necesita 600 iteraciones para llegar a la aproximación 3,14 y por entonces sólo avanza 0,01 cifras por iteración. Avanza MUY lentamente, pero lo que más gracia me hace es que ¡¡¡ sabemos que, igual que el resto de métodos, en el infinito también sale PI !!!
Sable | 18 de Junio de 2007 | 13:45
Tal vez la serie que propone Daniel sea muy lenta, pero ¡qué decir de la serie de Ramanujan! que se aproxima al valor de pi a velocidad de vértigo. Cada nuevo sumando proporciona 8 decimales exactos de pi. La fórmula se puede ver en
http://www.epsilones.com/paginas/i-formulas.html#formulas-piramanujan
Poco más arriba de este enlace aparece otra fórmula que proporciona 14 decimales…
Contestando al interrogante planteado en epsilones, Ramanujan ideó la fórmula investigando las funciones modulares(alguien de aquí lo entenderá…) Pueden ver la fórmula de otra manera
http://suanzes.iespana.es/ramanuj.htm
Y para los avanzados e interesados en más fórmulas de Ramanujan para el cálculo de pi
http://personal.auna.com/jguillera/ramaslides.pdf
Un gran genio y personaje. Y algo más de la vida de Ramanujan como no en http://criptociencia.blogspot.com/2007/05/calculo-mental.html
Jonas Castillo | 24 de Junio de 2007 | 3:48
Al estilo de Ramanujan
(150 + 1/200)^2
——————- = casi un entero (231)
pi ^4
Sebastian | 19 de Julio de 2007 | 17:07
Daniel, yo tuve la misma idea que vos….
de hecho, hice un programita en VB que itera basado en el problema de basilea y desp de 2.000.000 de iteraciones el resultado es este:
3. 14 15 92 17 61 25 09
Son 6 decimales para 2.000.000 de iteraciones
El problema, creo yo, es que la raiz cuadrada hace que pierda mucha potencia,
El metodo no es bueno para aproximar pi
OTERO | 25 de Agosto de 2007 | 20:08
Este mensaje va dirigido al que puso el comentario sobre una formula para el numero pi
con el nombre de rober con fecha 12 junio 2007
he intentado hacer funcionar tu formula en un programa de ordenador y no he obtenido el resultado esperado supongo que es debido a que la interpretacion de tu escrito no es la misma que la que tu posees.Quisiera que me dieses una explicacion mas detallada y clara sobre tu formula puedes hacerlo a traves de estos comentarios o mandarme un mensaje a la direccion de E-mail MSN messenger. oteropera@hotmail.com.Estoy interesado en cualquier formula o demostracion de formula para el numero pi por lo tanto si me envias una respuesta te estare muy agradecido.Y por ultimo quisiera hacer un comentario sobre tu formula si la variable Xi parecen ser el coseno de un angulo y la variable Yi el seno de un angulo segun mi interpretacion de tu formula la variable Li parece ser unicamente el radio de la circunferencia.Esperando una pronta respuesta.Cordiales saludos candido otero
OTERO | 25 de Agosto de 2007 | 20:13
Este mensaje es un añadido al anterior debido a un error mio al redactar el mensaje la direccion correcta de E-mail MSN MESSENGER es oteropera@hotmail.com
Dodovrosky Medrano | 18 de Septiembre de 2007 | 2:31
Talvez si los programadores utilizarían las ya clasicas fórmulas de Leibniz y Euler, les podría parecer más útil, en todo caso, yo no se mucho sobre programación pero en algoritmos matemáticos que justifiquen el valor de pi, talvez pueda dar algunas sugerencias
jose | 7 de Octubre de 2007 | 10:44
hola, me gustaría saber como puedo conseguir dibujar pi en la recta real de forma exacta. Gracias.
KATYA | 4 de Marzo de 2008 | 5:29
disculpa quisiera saber como ordena la siguientes funciones de un programa en java q es oriendato a objetos, acerca de estadistica para calcular la media,moda,desviacion estandar ,media armonica,el rango,varianza ,para datos agrupados claro. porfa aclara mis dudas.
ARTURO | 3 de Abril de 2008 | 22:51
gracias me ayudaron en mi trabajo
Joaquín Salazar | 21 de Abril de 2008 | 16:26
Hola soy alumno de bioquimica, pero aun asi amante de la matematica, y desarrollé un método bastante “bruto” para calcular pi a gusto.
Es bastante simple, el unico problema es que incluye la función seno, bueno soy un amateur, por lo tanto les dejo el desafio de seguir mejorándola..
Bueno acá va con explicación:
Si tienes una circunferencia de radio 1, le dibujas 4 radios que esten a 90º uno del otro y luego trazas diagonales entre el extremo de cada radio, por supuesto obtendremos 4 diagonales con valor squr (2), llamemos a esto una secuencia P, en su primer valor o sea P(1), esta nos lanzará el valor 4 squr (2). Si analizamos esta figura con mucha atención, podremos ver 4 triángulos rectángulos, su hipotenusa será 2, su cateto opuesto squr (2). y su cateto adyacente también. El asunto es que si somos aún mas observadores notaremos que este perimétro está dado por 4 sen 45º, o sea P(1) = 4 sen 45º. Sigamos la secuencia y a cada diagonal le buscamos un punto medio, y de este trazamos dos líneas que irán a los extremos de cada diagonal, así obteniendo 8 trazos. Ahora se ve más complicado pero no lo es. Solo hacemos el análisis anterior y notaremos que hay el doble de triángulos, y el ángulo q usamos anteriormente de redujo a la mitad. o sea si aplicamos el mismo concepto, P(2) = 8 sen 22,5º. Entonces si generalizamos obtendremos la siguiente fórmula…
Lim (2^(n + 1)) * sen (45º / 2^(n-1)) = pi
n->inf
El problema es que el seno nos complica la tarea… pero igualmente es efectiva en máquinas y calculadoras, y nos dará el valor de pi con la exactitud que queremos. Bueno espero que les sirva …
Bruno | 22 de Abril de 2008 | 14:19
José: no se puede. Fíjate los artículos sobre regla y compás.
Fran Rojas | 31 de Agosto de 2008 | 3:05
Para una libreria que estoy programando necesitaba calcular pi, y decidí usar la fórmula de Euler, pero pude comprobar lo lento que convergía.
Buscando por internet, encontré un interesante programa de ejemplo en Java que calculaba pi con precisión configurable.
Utiliza la fórmula:
pi/4 = 4 * arctan( 1/5 ) - arctan( 1/239 )
y para calcular el arcotangente, utiliza la serie:
arctan(x) = x - (x^3)/3 + (x^5)/5 - (x^7)/7 + (x^9)/9 …
que, al ser el argumento pequeño, converge bastante rápidamente.
Para lo que yo lo necesito, me viene perfecto.
(en mi PC calcula 10.000 decimales en menos de 2 segundos).
Seguramente no es el método más rápido para calcular pi, pero es sencillo y aceptablemente rápido.
el programa, lo podéis encontrar en:
http://blog.taragana.com/index.php/archive/calculate-pi-to-arbitrary-precision-sample-java-code/