Python中的Gauss-Legendre算法

Z时代
2024-01-10
分类：问答

我需要一些帮助来计算Pi。我正在尝试编写一个将Pi转换为X位数的python程序。我已经从python邮件列表中尝试了几种，但使用起来很慢。我已经阅读了有关Gauss-Legendre算法的文章，并且尝试将其移植到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()

回答：

您忘记了以下括号4*t：
```
pi = (a+b)**2 / (4*t)
```

您可以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()