Algoritmo Gauss-Legendre en python

algorithm pi

Question

Necesito ayuda para calcular Pi. Estoy tratando de escribir un programa de Python que calculará Pi a X dígitos. He intentado varios de la lista de correo de python, y es lento para mi uso. He leído sobre el algoritmo de Gauss-Legendre y he tratado de portarlo a Python sin éxito.

Estoy leyendo desde aquí , y agradecería cualquier opinión sobre dónde me estoy equivocando.

Salidas: 0.163991276262

from __future__ import division
import math
def square(x):return x*x
a = 1
b = 1/math.sqrt(2)
t = 1/4
x = 1
for i in range(1000):
y = a
a = (a+b)/2
b = math.sqrt(b*y)
t = t - x * square((y-a))
x = 2* x
pi = (square((a+b)))/4*t
print pi
raw_input()

Answer 1

1. Si desea calcular PI a 1000 dígitos, necesita usar un tipo de datos que admita 1000 dígitos de precisión (p. Ej., MxNumber )
2. Necesita calcular a, b, t y x hasta | ab | <10 ** - dígitos, no iterar dígitos veces.
3. Calcula cuadrado y pi como @JF sugiere.

Answer 2

pi = (square((a+b)))/4*t

debiera ser

pi = (square((a+b)))/(4*t)

Answer 3

1. Olvidó paréntesis alrededor de 4*t :

   pi = (a+b)**2 / (4*t)

2. Puede usar decimal para realizar cálculos con mayor precisión.

   #!/usr/bin/env python
   from __future__ import with_statement
   import decimal
   def pi_gauss_legendre():
   D = decimal.Decimal
   with decimal.localcontext() as ctx:
   ctx.prec += 2
   a, b, t, p = 1, 1/D(2).sqrt(), 1/D(4), 1
   pi = None
   while 1:
   an = (a + b) / 2
   b = (a * b).sqrt()
   t -= p * (a - an) * (a - an)
   a, p = an, 2*p
   piold = pi
   pi = (a + b) * (a + b) / (4 * t)
   if pi == piold: # equal within given precision
   break
   return +pi
   decimal.getcontext().prec = 100
   print pi_gauss_legendre()

Salida:

3.141592653589793238462643383279502884197169399375105820974944592307816406286208/
998628034825342117068