TL;DR
Tu supuesto "valor exacto de pi" contra el que te comparabas no era exacto. Sólo lo era hasta el decimal 15. De hecho, el valor de pi que computaba tu bucle era más correcto (y por eso se apartaba del valor de math.pi
a partir del decimal 15).
De todas formas, el algoritmo que usas para computar pi tiene serios problemas de estabilidad a partir de un punto. [Ver ampliación final para una posible solución]
Comprender Decimal
Lo primero, el tipo Decimal()
no es una fórmula milagrosa. Cuando conviertes un float
a Decimal
no vas a obtener instantáneamente más precisión de la que tenía el float
de entrada. Por ejemplo, si tú le pides representar con 51 decimales la cantidad 3.14, no vas a obtener más decimales de pi, es más, ni siquiera vas a obtener 3.1400000...00000 (49 ceros) sino
>>> print(decimal.Decimal(3.14))
3.140000000000000124344978758017532527446746826171875
Como ves, a partir de un punto aparecen decimales "raros". Esos decimales ya estaban ahí en el número 3.14 porque su representación en binario tiene infinitos decimales, y ya que el guardarlo en float
eso se trunca, se está cometiendo un error. Dicho de otro modo, cuando haces a = 3.14
(tipo float), en realidad el número que queda almacenado en a
es una aproximación de 3.14 (y en concreto, la aproximación mostrada más arriba).
Para evitar estos errores de redondeo no deberías pasarle un float
al constructor de Decimal()
, sino una cadena. Si hacemos Decimal("3.14")
, entonces sí, el valor almacenado será exacto:
>>> print(decimal.Decimal("3.14"))
3.14
# Almacena en realidad
# 3.140000000000000000000000000000000000000000000000000
Aunque a la hora de imprimirlo no muestra los ceros no significativos, en realidad están correctamente almacenados dentro de la representación Decimal.
Digo todo esto porque comienzas con:
pi2 = Decimal(math.pi)
y consideras que esta variable es "el valor exacto de pi", pero no lo es ya que ha sido generado a partir de math.pi
que es una aproximación válida sólo hasta el decimal 15, debido a que es un número de punto flotante que tiene errores inherentes como ya fue explicado aqui
Compara el valor de pi2
con un valor más exacto generado por Wolfram Alpha hasta el decimal 65:
Exacto con 65 cifras: 3.141592653589793 23846264338327950288419716939937510582097494459231
Decimal(math.pi): 3.141592653589793 115997963468544185161590576171875
Algunos detalles adicionales sobre tu programa
En tu programa tienes varios decimal.Decimal()
que podrían eliminarse.
Por ejemplo sobra aqui:
x1 = decimal.Decimal(Decimal(2)+a).sqrt()
que puede abreviarse a:
x1 = (a+2).sqrt()
puesto que a
ya era de tipo Decimal
, por lo que a+2
será también Decimal
(python convertirá 2
en Decimal
, de forma exacta por tratarse de un entero, y hará la suma). Ya que el resultado es un Decimal()
y aplicamos sobre él el método .sqrt()
de los Decimal()
, el resultado de esa raiz será a su vez otro Decimal()
que podemos asignar a x1
directamente, sin necesidad de convertirlo de nuevo en Decimal()
Ejecutando tu programa corregido
Vamos a ejecutar una versión simplificada (quitando lo decimal.Decimal()
superfluos) en la que nos comparamos con el valor "verdadero" de pi, al menos hasta la cifra 65, que es el que hemos obtenido de Wolfram Alpha
Este es el programa:
# @title pi
import math
from decimal import Decimal
import decimal
decimal.getcontext().prec = 65
a = Decimal(0)
valorN = 0
i = 3
pi = Decimal(0)
pi2 = Decimal("3.14159265358979323846264338327950288419716939937510582097494459231")
valorN = int(input('Ingrese el valor de n que sea mayor a 2 y menor a infinito: '))
a = Decimal(0)
while i <= valorN:
x1= (2+a).sqrt()
x = (2-x1).sqrt()
y = Decimal(pow(2,i-1))
pi = x*y
a = (2+a).sqrt()
i = i + 1
print('El valor de pi es ')
print(pi)
# este es el valor real #
print (pi2)
Al ejecutarlo para N=28 obtenemos resultados correctos hasta el decimal 15
Ingrese el valor de n que sea mayor a 2 y menor a infinito: 28
El valor de pi es
3.141592653589793 1667462219700150778696165489637681660703093062077
3.141592653589793 23846264338327950288419716939937510582097494459231
Pero ¿qué esperabas? En 28 iteraciones del algoritmo esa es la precisión que se ha logrado. Si quieres más precisión, podemos seguir iterando. Probemos por ejemplo 35 iteraciones:
Ingrese el valor de n que sea mayor a 2 y menor a infinito: 35
El valor de pi es
3.1415926535897932384 582661602928534441638628573935649839361286949
3.1415926535897932384 6264338327950288419716939937510582097494459231
Como vemos ha mejorado la precisión que ahora tenemos 19 decimales correctos (que ya es más de lo que podríamos tener en un float
).
Ocurre que N iteraciones no te da N decimales correctos. El algoritmo no converge tan rápido. Para lograr 19 decimales debes iterar 35 veces. Pero a partir de un punto seguir iterando ya no mejora el resultado. De hecho lo empeora!
El problema con este enfoque
Este no es un buen mecanismo para encontrar decimales de pi, pues es inestable. A medida que iteras, x1
es cada vez más próximo a 2, por lo que 2-x1
es cada vez más próximo a cero, y por tanto x
también.
De hecho llega un momento en que los 65 decimales no son suficientes para guardar x
(que sería inferior a 10^-65) y su valor se trunca a cero. A partir de ese punto todo falla, pues si x
es cero, también lo será x**y
y por tanto de pronto tu aproximación de pi se convierte en cero.
Esto pasa relativamente pronto (para n=110). Pero bastante antes de que eso ocurra, aparecen errores importantes en el cómputo de x**y
, ya que x
tomará valores muy pequeños, mientras que y
va tomando valores muy grandes (potencias de 2). En los valores pequeños de x
hay un error debido a que estamos truncándolo en su decimal 65, y ese error se ve amplificado al elevarlo a y
. El resultado es que antes de que se produzca la "catástrofe" que hace que pi
pase a valer cero, ya las aproximaciones previas eran cada vez peores en vez de ser cada vez mejores.
Por ejemplo, para n=105
la aproximación de pi
que salía era:
Ingrese el valor de n que sea mayor a 2 y menor a infinito: 105
El valor de pi es
3.14 21373846360117824379164873499126756724026658776749425705467686
3.14 15926535897932384626433832795028841971693993751058209749445923
¡Toma! ¡Está mal ya a partir del tercer decimal!
Parece entonces que a medida que aumentamos N nos vamos aproximando a pi, pero sólo hasta un punto a partir del cual nos alejamos de nuevo debido a los errores causados porque x
es muy muy pequeño y no cabe en 65 decimales, hasta que finalmente el error es ya catastrófico pues nuestro pi vale cero.
Curiosidades finales
Como curiosidad he hecho una gráfica de cuál es el error absoluto cometido (diferencia entre la aproximación de pi y su valor "verdadero") a medida que aumentamos N (iteraciones del bucle):
He usado escala logarítmica en el eje Y. Vemos como hasta N aproximadamente igual a 53 la cosa iba bien. El error se iba reduciendo exponencialmente, pero a partir de ahi la cosa se torció. Los errores de redondeo debido a la inestabilidad x**y
hicieron que el error fuera creciendo de nuevo (de forma irregular, un tanto caótica), hasta llegar a un error igual a 3.1415... cuando pi
pasa a valer 0 (y a partir de ese punto el error ya se queda en 3.1415.. porque toda iteración siguiente del bucle produce pi = 0
.
Para N=53 se produce la mejor aproximación de todas, que es la siguiente:
Ingrese el valor de n que sea mayor a 2 y menor a infinito: 53
El valor de pi es
3.141592653589793238462643383279 4390490856382359047953080772238233
3.141592653589793238462643383279 5028841971693993751058209749445923
correcta hasta el decimal 30.
He pintado otra gráfica que muestra cuántos dígitos correctos tiene la aproximación a medida que aumentas valorN
. Es la siguiente:
Donde se ve de nuevo que el máximo se alcanza para N alrededor de 55 y también que hasta ese momento la progresión no es lineal. En ocasiones pasar de N a N+1 aumenta el número de dígitos correctos en 1, otras veces no lo aumenta. Y a partir de 55 el número de dígitos correctos disminuye, también "a saltos".
Parece por tanto que con una precisión de 65 llegas correctamente hasta el decimal 30. Ampliando la precisión podrías obtener más decimales correctos.
Ampliación
Aparentemente, según algunas pruebas que he estado haciendo, para obtener N dígitos correctos, tienes que trabajar con una precisión de al menos 2*(N+3) (es un heurístico para el que no tengo justificación). Así, para obtener correctamente 65 dígitos deberías trabajar con una precisión de
136.
En lugar de pedirle al usuario el número de iteraciones a realizar, podemos detectar automáticamente cuándo parar, y será cuando el valor que nos sale de pi
sea igual al de la iteración anterior (a partir de ese punto la aproximación comenzaría a empeorar).
Este sería por tanto el código:
import math
from decimal import Decimal
import decimal
valorN = int(input('Ingrese el valor de dígitos deseado: '))
# Trabajar con la precisión adecuada
decimal.getcontext().prec = 2*(valorN+3)
a = Decimal(0)
i = 3
pi = Decimal(0)
a = Decimal(0)
while True:
piold = pi
x1= (2+a).sqrt()
x = (2-x1).sqrt()
y = Decimal(pow(2,i-1))
pi = x*y
if piold==pi:
break
a = (2+a).sqrt()
i = i + 1
print(f"Detenido tras {i} iteraciones")
print('El valor aproximado de pi es ')
print(pi)
# este es el valor real con 65 digitos #
print("3.1415926535897932384626433832795028841971693993751058209749445923")
Resultado de la ejecución:
Ingrese el valor de dígitos deseado: 65
Detenido tras 116 iteraciones
El valor aproximado de pi es
3.141592653589793238462643383279502884197169399375105820974944592307844062894732236752971209989019846026985971418316643805175816308465281
3.1415926535897932384626433832795028841971693993751058209749445923
Naturalmente el resto de decimales que salen a partir de la posición 65 estarían mal. Podríamos truncar el resultado con:
print(str(pi)[:valorN+1])