Gauss-Legendre Algorithm in python

Question

I need some help calculating Pi. I am trying to write a python program that will calculate Pi to X digits. I have tried several from the python mailing list, and it is to slow for my use. I have read about the Gauss-Legendre Algorithm, and I have tried porting it to Python with no success.

I am reading from Here, and I would appreciate any input as to where I am going wrong!

It outputs: 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()

Traitor

Answer 1

1. You forgot parentheses around 4*t:

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

2. You can use decimal to perform calculation with higher precision.

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

Borneol

Unless you change to use some other data type the best you can get is 24 or 53 digits of precision using either 32- or 64-bit floating point arithmetic. For more info see en.wikipedia.org/wiki/IEEE_754. – Disrobe

@tvanfosson: I've posted version that uses decimal. It allows arbitrary precision. – Borneol

+1 -- didn't know that Python had decimal and mxNumber was the first item that popped up in Google. – Disrobe

Answer 2

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

should be

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

Sweetie

Answer 3

1. If you want to calculate PI to 1000 digits you need to use a data type that supports 1000 digits of precision (e.g., mxNumber)
2. You need to calculate a,b,t, and x until |a-b| < 10**-digits, not iterate digits times.
3. Calculate square and pi as @J.F. suggests.

Disrobe

decimal module is sufficient for 1000 digits. – Borneol