8

En una pregunta reciente se plantea un algoritmo para aproximar decimales de pi. El algoritmo básicamente se reduce a iterar la siguiente relación de recurrencia:

num_iter = 10
a = 0
i = 1
while i <= num_iter:
  a = sqrt(2+a)
  x = sqrt(2-a)
  pi = x * 2**(i+1)
  i += 1

Esta relación funciona bien para valores bajos de num_iter, pero enseguida empiezan a producirse errores de redondeo, los cuales eran el objetivo de aquella pregunta, y que pueden subsanarse usando representaciones de coma fija (decimal), a costa de incrementar la memoria necesaria.

Mi pregunta ahora es ¿qué algoritmo es este que se está usando? No he podido encontrar referencias al mismo. Por otro lado no parece ser un algoritmo muy bueno (como se mostró en la respuesta a aquella pregunta) ¿Qué otros algoritmos hay mejores para computar pi? ¿En qué sentido son mejores?

3
  • muchas gracias por tu ayuda, en respuesta a tu pregunta, te comento que el método que usé para hallar PI fue el de aproximaciones planteada por Arquímedes, si bien hay mejores formas como acabo de enterarme, lo hice por que fue planteado por el profesor de razonamiento matemático, en una asamblea en el colegio por el día de PI, por cuestiones de aplicar la geometría, y como estoy aprendiendo Python, quise comprobarlo haciendo un programa para tal fin quedo agradecido por tu valioso aporte en mi aprendizaje
    – Marcus
    el 1 abr. 2021 a las 1:01
  • nota: el algoritmo que planteas tiene un error, pi = x * 2**(i+1), debe ser pi = x * 2**(i-1)
    – Marcus
    el 1 abr. 2021 a las 1:02
  • @Marcus No, la formula es correcta porque empiezo en i=1 en vez de 3, pues me pareció más natural.
    – abulafia
    el 1 abr. 2021 a las 6:29

1 respuesta 1

7

El algoritmo usado en aquella pregunta no es ninguno de los que Wikipedia muestra en su entrada sobre aproximaciones de pi.

Matemáticamente este algoritmo puede ser expresado mediante la siguiente relación de recurrencia:

Ecuacion

que puede ser visualizada más fácilmente si usamos la versión "expandida" para los casos n=2 y n=3 por ejemplo:

expansion

La única referencia que he encontrado a esta aproximación de pi proviene de una respuesta en Math.stackexchange, la cual a su vez era simplemente una copia de un post en Facebook, que incluía la siguiente imagen:

Foto papel

Por tanto parece una aproximación geométrica, diseñada por ser fácil de visualizar y no por ser precisamente eficiente. Hay aproximaciones mucho más eficientes de las cuales destacaré sólo un par (pues el tema de computar pi eficientemente es ciertamente un mundo).

Aproximación de Gauss-Legendre

Esta aproximación converge mucho más rápidamente que la usada en esta pregunta. En la usada en esta pregunta, el número de decimales correctos se incrementa tan solo en 1 (o a veces en 0) con cada nueva iteración. En cambio la aproximación de Gauss-Legendre el número de dígitos correctos se duplica en cada iteración, por lo que tras unas pocas ya se tiene un elevado número de decimales.

Los detalles matemáticos son complejos (se pueden consultar en Wikipedia), y aquí daré sólo una implementación en Python (adaptada de esta respuesta en SO):

from decimal import Decimal, getcontext

num_iter = int(input("iteraciones? "))
getcontext().prec = 65
a = 1
b = 1/Decimal(2).sqrt()
t = 1/Decimal(4)
p = 1                
i = 0
while i < num_iter:
  an = (a + b) / 2
  b  = (a * b).sqrt()
  t -= p * (a - an) * (a - an)
  a = an
  pi = (a + b) * (a + b) / (4 * t)
  p  = 2*p
  i = i + 1

# Aproximacion
print(pi)  
# Valor exacto (hasta decimal 65, obtenido de WolframAlpha)
print("3.1415926535897932384626433832795028841971693993751058209749445923")

Con solo tres iteraciones ya se obtienen muchos decimales correctos!

iteraciones? 3
3.141592653589793238 2795127748018639743812255048354469357873307022
3.141592653589793238 4626433832795028841971693993751058209749445923

Y con una iteración más se duplica el número de decimales correctos

iteraciones? 4
3.1415926535897932384626433832795028841971 146782836489215566171072
3.1415926535897932384626433832795028841971 693993751058209749445923

ya con cinco o más tenemos los 64 dígitos correctos (el último dígito no sale nunca correcto debido al error de redondeo causado por trabajar solo con 65 decimales)

De hecho, una pequeña variación del programa anterior nos puede permitir obtener cualquier número de decimales deseado. Basta trabajar con una precisión de un par de decimales más (por ejemplo si queremos 500 dígitos de pi, trabajamos con 502) para eliminar ese problema en el redondeo final. Y en lugar de iterar un número prefijado de veces, iteramos hasta que dos iteraciones sucesivas produzcan el mismo resultado, deteniéndonos entonces:

from decimal import Decimal, getcontext

n_digitos = int(input("Número de dígitos deseado? "))
getcontext().prec = n_digitos+2
a = 1
b = 1/Decimal(2).sqrt()
t = 1/Decimal(4)
p = 1                
i = 0
pi_anterior = None
while True:
  an = (a + b) / 2
  b  = (a * b).sqrt()
  t -= p * (a - an) * (a - an)
  a = an
  pi = (a + b) * (a + b) / (4 * t)
  if pi == pi_anterior:
    break
  p  = 2*p
  pi_anterior = pi
  i = i + 1

print("Obtenido resultado en {} iteraciones".format(i))
# Aproximacion
print(str(pi)[:-2])  # Quitar los dos últimos, que serían incorrectos  
# Valor exacto (hasta decimal 65, obtenido de WolframAlpha)
print("3.1415926535897932384626433832795028841971693993751058209749445923")

Número de dígitos deseado? 65
Obtenido resultado en 5 iteraciones
3.1415926535897932384626433832795028841971693993751058209749445923
3.1415926535897932384626433832795028841971693993751058209749445923

Vamos a lo grande!

Número de dígitos deseado? 1000
Obtenido resultado en 9 iteraciones
3.14159265358979323846264338327950288419716939937510582097494459230781640628620899
8628034825342117067982148086513282306647093844609550582231725359408128481117450284
1027019385211055596446229489549303819644288109756659334461284756482337867831652712
0190914564856692346034861045432664821339360726024914127372458700660631558817488152
0920962829254091715364367892590360011330530548820466521384146951941511609433057270
3657595919530921861173819326117931051185480744623799627495673518857527248912279381
8301194912983367336244065664308602139494639522473719070217986094370277053921717629
3176752384674818467669405132000568127145263560827785771342757789609173637178721468
4409012249534301465495853710507922796892589235420199561121290219608640344181598136
2977477130996051870721134999999837297804995105973173281609631859502445945534690830
2642522308253344685035261931188171010003137838752886587533208381420617177669147303
5982534904287554687311595628638823537875937519577818577805321712268066130019278766
11195909216420198
3.1415926535897932384626433832795028841971693993751058209749445923

Algoritmos basados en obtener cifras

Si el algoritmo anterior se basa en aproximaciones (cada iteración produce un número más cercano a pi), hay otros algoritmos que se basan en ir obteniendo cifras de pi. La principal ventaja de este tipo de algoritmos sería que no requieren tanto espacio de almacenamiento, pues una vez un dígito ha sido computado, se imprime y no es necesario mantener en la memoria todos los dígitos anteriores, a la hora de computar el siguiente. Esto es una ventaja frente a métodos como el antes visto de Gauss-Legendre, en el cual si queremos llegar a obtener una precisión de un millón de cifras, debemos trabajar con un millón de dígitos ya desde el principio.

Uno algoritmo muy notable es el basado en la fórmula de Bailey–Borwein–Plouffe, pues permite computar el dígito que aparece en posición N sin necesidad de computar los N-1 anteriores. Por desgracia ese algoritmo produce los dígitos en base hexadecimal, y no en la base 10 a la que estamos acostumbrados. Conocer el dígito N-ésimo en hexadecimal no nos sirve para conocer un dígito concreto de la expansión decimal. Lo que sí es posible es escribir pi con N dígitos (hexadecimales), lo cual sería una aproximación de pi, y convertir ese número a decimal, para obtener otra aproximación en base 10, que sería correcta con una precisión de 2-4*N

No obstante esto de guardar los N primeros dígitos hexadecimales para pasarlos a decimal estropea esa ventaja de "no tener que computar ni almacenar los dígitos previos".

Otros métodos para ir obteniendo pi "dígito a dígito" se conocen globalmente como algoritmos Spigot. Uno de ellos (que he adaptado de esta respuesta en SO), dice ser bastante rápido, aunque no es nada evidente entender qué hace (ni aún leyendo el artículo en que se basa), por lo que podría incluso ser incorrecto, aunque parece dar correctamente los primeros 65 decimales:

def pi_digits():
  k,a,b,a1,b1 = 2,4,1,12,4
  while True:
    p,q,k = k * k, 2 * k + 1, k + 1
    a,b,a1,b1 = a1, b1, p*a + q*a1, p*b + q*b1
    d,d1 = a/b, a1/b1
    while d == d1:
        yield int(d)
        a,a1 = 10*(a % b), 10*(a1 % b1)
        d,d1 = a//b, a1//b1

num_decimales = 65
for i, digit in enumerate(pi_digits()):
  print(digit, end="")
  if i==0: print(".", end="")  # El punto decimal
  if i==num_decimales-1: break

print()
# Valor exacto para comparación
print("3.1415926535897932384626433832795028841971693993751058209749445923")

Resultado de la ejecución:

3.1415926535897932384626433832795028841971693993751058209749445923
3.1415926535897932384626433832795028841971693993751058209749445923

Este método no converge tan deprisa como el antes visto, pues para obtener 65 decimales ha iterado 65 veces. La función pi_digits() es un generador infinito, cada vez que se le llama proporciona el siguiente dígito de pi. El for desde el cual la uso tiene un break para detenerse tras 65 llamadas (de lo contrario sería infinito).

La ventaja de este método es que no depende de la precisión ni requiere decimal, pues trabaja con enteros, y que es más eficiente en memoria al no necesitar retener todos los dígitos de pi que haya computado hasta el momento, sino solo unas pocas variables (p,q,k,a,b,a1,1,d,d1). A cambio es más lento a partir de un cierto decimal.

Si le ponemos a computar 1000 dígitos vemos que obtiene exactamente el mismo resultado que el algoritmo de Gauss-Legendre, y tarda el doble. Si intentamos obtener 100.000 dígitos de pi, Gauss-Legendre tarda 21s, mientras que el algoritmo Spigot tarda del orden de 3 minutos. Cada dígito adicional va tardando un poco más.

Tu Respuesta

By clicking “Publica tu respuesta”, you agree to our terms of service and acknowledge you have read our privacy policy.

¿No es la respuesta que buscas? Examina otras preguntas con la etiqueta o formula tu propia pregunta.