Python计算圆周率?

20
我是一名Python初学者,想要计算圆周率。我尝试使用Chudnovsky算法,因为听说它比其他算法更快。
以下是我的代码:
from math import factorial
from decimal import Decimal, getcontext

getcontext().prec=100

def calc(n):
    t= Decimal(0)
    pi = Decimal(0)
    deno= Decimal(0)
    k = 0
    for k in range(n):
        t = ((-1)**k)*(factorial(6*k))*(13591409+545140134*k)
        deno = factorial(3*k)*(factorial(k)**3)*(640320**(3*k))
        pi += Decimal(t)/Decimal(deno)                                   
    pi = pi * Decimal(12)/Decimal(640320**(1.5))
    pi = 1/pi
    return pi

print calc(25)

由于某些原因,这段代码仅提供了15位小数的pi值,与可接受的值相比有所差异。我尝试通过增加精度值来解决此问题;这会增加数字的位数,但只有前15位是准确的。我尝试改变算法计算方式,但也没有起作用。因此,我的问题是,是否有办法使此代码更加准确,还是我必须使用另一种算法?我需要帮助,因为我不知道如何在Python中操作这么多位数。我希望能够控制程序确定和显示的(正确)数字数量 - 无论是10、100、1000等。

仅供比较,这是一些可行的Chudnovsky代码:http://www.craig-wood.com/nick/articles/pi-chudnovsky/ - runDOSrun
准确性受 Python 中 decimal 包的默认精度限制。正确配置即可解决问题。请参阅 https://docs.python.org/3/library/decimal.html。 - 4dummies
3个回答

21

看起来你在这一行代码中失去了精度:

pi = pi * Decimal(12)/Decimal(640320**(1.5))

尝试使用:

pi = pi * Decimal(12)/Decimal(640320**Decimal(1.5))

这是因为虽然Python可以处理任意大小的整数,但对浮点数的处理效果不佳。

奖励

使用另一个算法(BBP公式)的单行实现:

from decimal import Decimal, getcontext
getcontext().prec=100
print sum(1/Decimal(16)**k * 
          (Decimal(4)/(8*k+1) - 
           Decimal(2)/(8*k+4) - 
           Decimal(1)/(8*k+5) -
           Decimal(1)/(8*k+6)) for k in range(100))

非常感谢!它起作用了。对我来说,这样一个小改变的差异是不可思议的。 - Norrin Rad
@Juan Lopes 你好,能帮我写一下计算圆周率的BBP公式的伪代码吗? - Sumit Khanduri
还没有人指出来,但是BBP公式是一个十六进制的自动计算算法,而不是十进制的,对吧?我在维基百科上引用了这个。 - Jeff

12

对于那些只是想在Python中获得任意精度π的解决方案的人(源代码稍加修改):

import decimal

def pi():
    """
    Compute Pi to the current precision.

    Examples
    --------
    >>> print(pi())
    3.141592653589793238462643383

    Notes
    -----
    Taken from https://docs.python.org/3/library/decimal.html#recipes
    """
    decimal.getcontext().prec += 2  # extra digits for intermediate steps
    three = decimal.Decimal(3)      # substitute "three=3.0" for regular floats
    lasts, t, s, n, na, d, da = 0, three, 3, 1, 0, 0, 24
    while s != lasts:
        lasts = s
        n, na = n + na, na + 8
        d, da = d + da, da + 32
        t = (t * n) / d
        s += t
    decimal.getcontext().prec -= 2
    return +s               # unary plus applies the new precision

decimal.getcontext().prec = 1000
pi = pi()

1
这段代码出现在Python库文档中的decimal - ShreevatsaR

0
from decimal import *

#Sets decimal to 25 digits of precision
getcontext().prec = 25

def factorial(n):
    if n<1:
        return 1
    else:
        return n * factorial(n-1)

def plouffBig(n): #http://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(1)/(16**k))*((Decimal(4)/(8*k+1))-(Decimal(2)/(8*k+4))-(Decimal(1)/(8*k+5))-(Decimal(1)/(8*k+6)))
        k += 1
    return pi

def bellardBig(n): #http://en.wikipedia.org/wiki/Bellard%27s_formula
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(-1)**k/(1024**k))*( Decimal(256)/(10*k+1) + Decimal(1)/(10*k+9) - Decimal(64)/(10*k+3) - Decimal(32)/(4*k+1) - Decimal(4)/(10*k+5) - Decimal(4)/(10*k+7) -Decimal(1)/(4*k+3))
        k += 1
    pi = pi * 1/(2**6)
    return pi

def chudnovskyBig(n): #http://en.wikipedia.org/wiki/Chudnovsky_algorithm
    pi = Decimal(0)
    k = 0
    while k < n:
        pi += (Decimal(-1)**k)*(Decimal(factorial(6*k))/((factorial(k)**3)*(factorial(3*k)))* (13591409+545140134*k)/(640320**(3*k)))
        k += 1
    pi = pi * Decimal(10005).sqrt()/4270934400
    pi = pi**(-1)
    return pi
print "\t\t\t Plouff \t\t Bellard \t\t\t Chudnovsky"
for i in xrange(1,20):
    print "Iteration number ",i, " ", plouffBig(i), " " , bellardBig(i)," ", chudnovskyBig(i)

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接