Gauss-Legendre Algorithm em python

Question

Preciso de ajuda cálculo Pi. Estou tentando escrever um programa python que irá calcular Pi para um dígito X. Eu tentei vários da lista de discussão python, e é para retardar para meu uso. Eu li sobre o Gauss-Legendre Algoritmo , e eu tentei portá-la para Python com sem sucesso.

Estou lendo de Aqui , e eu gostaria de receber qualquer entrada como a onde estou indo errado!

Ele produz: 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. Você esqueceu parênteses em torno 4*t:

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

2. Você pode usar decimal para realizar cálculo com maior precisão.

   #!/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()
   

Output:

3.141592653589793238462643383279502884197169399375105820974944592307816406286208\
    998628034825342117068

Answer 2

1. Se você quiser calcular PI 1000 dígitos você precisa usar um tipo de dados que suporta 1000 dígitos de precisão (por exemplo, mxNumber )
2. Você precisa calcular a, b, t, e x até | a-b | <10 ** -. Dígitos, não iterate dígitos vezes
3. quadrado Calcular e pi como @ J.F. sugere.

Answer 3

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

deve ser

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