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()

jfs · Accepted Answer · 2010-11-02 05:48:33Z

You forgot parentheses around 4*t:
```
pi = (a+b)**2 / (4*t) 
```

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

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.
@tvanfosson: I've posted version that uses decimal. It allows arbitrary precision.
+1 -- didn't know that Python had decimal and mxNumber was the first item that popped up in Google.

tvanfosson · Accepted Answer · 2008-12-07 18:22:15Z

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)
You need to calculate a,b,t, and x until |a-b| < 10**-digits, not iterate digits times.
Calculate square and pi as @J.F. suggests.

Bart Kiers · Accepted Answer · 2010-06-21 19:11:08Z

3

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

should be

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

edited Jun 21, 2010 at 19:11

Bart Kiers

171k38 gold badges307 silver badges297 bronze badges

answered Dec 7, 2008 at 16:37

Stefan Pantos

Collectives™ on Stack Overflow

Gauss-Legendre Algorithm in python

3 Answers 3

3 Comments

1 Comment

Comments

Linked

Hot Network Questions