用Python计算圆周率的前1000位数字

Question

我一直在考虑这个问题，但是我无法解决它。也许你可以帮助我。问题是我的代码不能在Python编程语言中输出1000位的π。

以下是我的代码：

def make_pi():
    q, r, t, k, m, x = 1, 0, 1, 1, 3, 3
    while True:
        if 4 * q + r - t < m * t:
            yield m
            q, r, t, k, m, x = (10*q, 10*(r-m*t), t, k, (10*(3*q+r))//t - 10*m, x)
        else:
            q, r, t, k, m, x = (q*k, (2*q+r)*x, t*x, k+1, (q*(7*k+2)+r*x)//(t*x), x+2)

digits = make_pi()
pi_list = []
my_array = []
for i in range(1000):
    my_array.append(str("hello, I'm an element in an array \n" ))
big_string = "".join(my_array)

print "here is a big string:\n %s" % big_string 

我知道这段代码可以修复，但我不确定需要修复什么......print语句中的here是一个大字符串，my_array.append(str("hello, im an element in an array \n))目前只是一个占位符。我知道所有代码的工作原理，但就像我之前说的那样，我无法使它输出那段代码。

- bobimo

---

这看起来像是 pi 出水口算法的一个版本，它确实是吗？ - Dan D.

你能否再明确一下问题所在？行为与你预期的有何不同？ - Scott Hunter

这段代码看起来非常类似于这里的代码和这里的代码。 - Sven Marnach

这个问题看起来非常熟悉。它之前在这里被问过并被删除了吗？ - Kris Harper

相关：Python中的高斯-勒让德算法 - jfs

Answer 1

如果您不想实现自己的算法，可以使用mpmath。

try:
    # import version included with old SymPy
    from sympy.mpmath import mp
except ImportError:
    # import newer version
    from mpmath import mp

mp.dps = 1000  # set number of digits
print(mp.pi)   # print pi to a thousand places

更新：代码支持旧版和新版SymPy（请参见评论）。*

参考资料

- Garrett Hyde

---

当你在类似于Anaconda创建的环境中安装sympy时，它会单独安装mpmath，并且导入语句变为from mpmath import mp。这个答案中的代码在现代版本的sympy中不再起作用了。现在没有sympy.mpmath了。 - Zelphir Kaltstahl

Answer 2

运行这个

def make_pi():
    q, r, t, k, m, x = 1, 0, 1, 1, 3, 3
    for j in range(1000):
        if 4 * q + r - t < m * t:
            yield m
            q, r, t, k, m, x = 10*q, 10*(r-m*t), t, k, (10*(3*q+r))//t - 10*m, x
        else:
            q, r, t, k, m, x = q*k, (2*q+r)*x, t*x, k+1, (q*(7*k+2)+r*x)//(t*x), x+2


my_array = []

for i in make_pi():
    my_array.append(str(i))

my_array = my_array[:1] + ['.'] + my_array[1:]
big_string = "".join(my_array)
print "here is a big string:\n %s" % big_string 

请点击链接了解关于yield操作符的内容："yield"关键字是什么意思？

以下是答案：

3.1415926535897932384626433832795028841971693993751058209749445923078164062862089986280348253421170679821480865132823066470938446095505822317253594081284811174502841027019385211055596446229489549303819644288109756659334461284756482337

- batbaatar

---

我有一个疯狂的想法，让树莓派在24/7运行数年后计算出数百万位的圆周率。这个问题/答案中提供的信息似乎非常有帮助。感谢您发布这个。我用50000的范围运行了代码，它完成了，但我不知道它是否准确。没有生成错误。需要想办法将部分结果卸载到磁盘而不是保留在内存中。 - WyomingGeezer

"for j in range(1000)" 使其生成的数字远少于1000个。应该改回“while True”或其他条件。 - Don Hatch

我真的不喜欢这段代码，首先 - 你改变了“make_pi”函数并加入了一些限制，甚至没有尝试将其提取为参数，但这并不是必要的 - make_pi()返回一个生成器，你可以使用.next()从中读取多次。其次，为什么要调用两次make_pi？ - Jerzyk

除了尝试之外，有没有任何解释说明这段代码为什么应该给出PI？它非常难以阅读。它让我想起了用于在对数时间内生成斐波那契数列的变换。 - Zelphir Kaltstahl

你为什么声明了 pi_list 并赋值 digits 却没有使用它们？这只会让人感到困惑。 - Arthur Dent

通过执行 j in range(4319)，这个答案将会给出1000位数字。 - wazawoo

Answer 3

正如评论中所指出的那样，被接受的答案是不正确的。

原帖作者的代码似乎基于从这里复制的Spigot算法实现。

为了修复代码以满足原帖作者的问题（尽管我已经重命名变量和函数以匹配原始源代码），一个解决方案可能是：

#!/usr/bin/env python

DIGITS = 1000

def pi_digits(x):
    """Generate x digits of Pi."""
    q,r,t,k,n,l = 1,0,1,1,3,3
    while x >= 0:
        if 4*q+r-t < x*t:
            yield n
            x -= 1
            q,r,t,k,n,l = 10*q, 10*(r-n*t), t, k, (10*(3*q + r))/t-10*n, l
        else:
            q,r,t,k,n,l = q*k, (2*q+r)*l, t*l, k+1, (q*(7*k+2)+r*l)/(t*l), l+2

digits = [str(n) for n in list(pi_digits(DIGITS))]
print("%s.%s\n" % (digits.pop(0), "".join(digits)))

此外，这里有一个更快的实现，也显然基于Spigot算法：

#!/usr/bin/env python

DIGITS = 1000

def pi_digits(x):
    """Generate x digits of Pi."""
    k,a,b,a1,b1 = 2,4,1,12,4
    while x > 0:
        p,q,k = k * k, 2 * k + 1, k + 1
        a,b,a1,b1 = a1, b1, p*a + q*a1, p*b + q*b1
        d,d1 = a/b, a1/b1
        while d == d1 and x > 0:
            yield int(d)
            x -= 1
            a,a1 = 10*(a % b), 10*(a1 % b1)
            d,d1 = a/b, a1/b1

digits = [str(n) for n in list(pi_digits(DIGITS))]
print("%s.%s\n" % (digits.pop(0), "".join(digits)))

我在这个在线圆周率数字生成器上多次测试了两种方法。

所有功劳归deeplook的这个代码片段所有。

* 基于测试10,000个数字，其中我分别用了约7秒和1秒的时间。

- Alex Harvey

Answer 4

如果需要获取长达1百万位数的π值，请使用math_pi 模块（注意：我是该模块的作者）

通过pip安装：

pip install math-pi

在Python中：

>>> import math_pi
>>> print(math_pi.pi(b=1000))
3.1415926535...

- mathcat

Answer 5

最初的回答： 来自Fabrice Bellard网站的 Pi Computation 算法。非常抱歉，这个实现方式比较直接。计算1000位小数足够快（我只用了0.1秒），但是计算10000位就不太快了 - 要71秒 :-(

import time
from decimal import Decimal, getcontext

def compute(n):
    getcontext().prec = n
    res = Decimal(0)
    for i in range(n):
        a = Decimal(1)/(16**i)
        b = Decimal(4)/(8*i+1)
        c = Decimal(2)/(8*i+4)
        d = Decimal(1)/(8*i+5)
        e = Decimal(1)/(8*i+6)
        r = a*(b-c-d-e)
        res += r
    return res

if __name__ == "__main__":
    t1 = time.time()
    res = compute(1000)
    dt = time.time()-t1
    print(res)
    print(dt)

- Guest007

Answer 6

我是在5-6年前使用下面的公式解决的。 Machin-like formula 维基百科: https://en.wikipedia.org/wiki/Machin-like_formula 对于代码质量不佳我很抱歉，变量名可能没有意义。

#-*- coding: utf-8 -*-

# Author:    Fatih Mert Doğancan
# Date:      02.12.2014

def arccot(x, u):
    sum = ussu = u // x
    n = 3
    sign = -1
    while 1:
        ussu = ussu // (x*x)
        term = ussu // n
        if not term:
            break
        sum += sign * term
        sign = -sign
        n += 2
    return sum

def pi(basamak):
    u = 10**(basamak+10)
    pi = 4 * (4*arccot(5,u) - arccot(239,u))
    return pi // 10**10

if __name__ == "__main__":
    print pi(1000) # 1000 

- Fatih Mert Doğancan

Answer 7

在这里，您可以检查您的程序是否正确输出了1000个数字： http://spoj.com/CONSTANT

当然，您也可以使用diff或tc，但是您需要从其他地方复制这1000个数字，而在此处您只需提交您的程序并检查分数是否大于999即可。

您可以尝试在那里打印更多的数字，从而获得更多的分数。也许您会喜欢它。

- user1

Answer 8

我不熟悉你的算法。它是BBP的实现吗？

无论如何，你的make_pi是一个生成器。尝试在for循环中使用它：

for digit in make_pi():
    print digit

请注意，此循环是无限的：make_pi() 永远不会抛出 StopIteration。

- Managu

---

该算法基于Spigot算法：https://en.wikipedia.org/wiki/Spigot_algorithm - Avizipi

Answer 9

华莱士公式可以得到3.141592661439964，但需要一种更有效的方法来解决这个问题。

https://www.youtube.com/watch?v=EZSiQv_G9HM

现在是我的代码

x, y, summing = 2, 3, 4

for count in range (0,100000000):
    summing *= (x/y)
    x += 2
    summing *= (x/y)
    y += 2

print (summing)

- tanny

---

由于整数除法，这仅适用于Python3。 - lenik

@lenik，但是Python 2.x现在已经停止支持。 - TheTechRobo the Nerd

Answer 10

这个做你想要的吗？

i = 0;
pi_str = ""
for x in make_pi():
    pi_str += str(x)
    i += 1
    if i == 1001:
        break

print "pi= %s.%s" % (pi_str[0],pi_str[1:])

- Scott Hunter