Python中的高斯-勒让德算法

Question

我需要一些帮助来计算圆周率（Pi）。我想写一个Python程序，可以计算出圆周率到X位数字。我尝试过几个从Python邮件列表上找到的方法，但对我来说速度太慢了。

我读过关于高斯-勒让德算法的资料，也尝试把它移植到Python上，但没有成功。

我在这个网站上阅读相关内容，如果你能告诉我哪里出错了，我会非常感激！

它输出的结果是：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

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

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

应该是

Answer 2

1. 如果你想计算出1000位的圆周率（PI），你需要使用一种可以支持1000位精度的数据类型（比如mxNumber）。
2. 你需要计算a、b、t和x，直到|a-b|小于10的负digits次方，而不是简单地循环计算digits次。
3. 按照@J.F.的建议来计算平方和圆周率。

Answer 3

1. 你忘了在 4*t 周围加上括号：

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

2. 你可以使用 decimal 来进行更精确的计算。

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

输出：

3.141592653589793238462643383279502884197169399375105820974944592307816406286208\
    998628034825342117068