Python 中的 BBP 算法 - 使用任意精度算术

Question

我正在写一篇关于 pi 计算的学期论文。虽然我已经完成了理论站点，但我现在正在努力用 Python 实现 BBP 算法。

您可以在这里找到 BBP 算法： http://en.wikipedia.org/wiki/Bailey%E2%80%93Borwein%E2%80%93Plouffe_formula

这是我在 Python 中的实现：

from sympy.mpmath import *
pi = mpf(0)
mp.dps = 30;
for n in range(0,500):
    pi = pi + mpf((1/16**n)*(4/(8*n+1)-  2/(8*n+4)-  1/(8*n+5)-  1/(8*n+6)) )


print(pi)

我的问题是，无论我将 k 设置多高，或者将 pi 的小数位数设置多高，我都无法获得比 16 位数字更高的精度。

我使用mpmath来获得更高的精度，因为我之前遇到了一些问题。

我如何改进我的代码，以便获得更多数字？

Answer 1

默认情况下，python 将使用 IEEE-754 定义的标准浮点数。它的精度约为 12 位数字，可以表示低至2 ^-1022mpf的数字，现在您可以通过在过程的早期调用运算符来解决此问题，因此更具可读性和更精确的版本是：

from sympy.mpmath import *
pi = mpf(0)
mp.dps = 30;
for n in range(0,500):
    u = 4.0/(8*n+1)-2.0/(8*n+4)-1.0/(8*n+5)-1.0/(8*n+6)
    u = mpf(u)
    d = mpf(16.0/1.0)**n
    pi += u/d

print(pi)

然而，当您计算u零件时，这仍然存在问题。要准确地做到这一点，您可以使用：

from sympy.mpmath import *
pi = mpf(0)
mp.dps = 50;
for n in range(0,50000):
    u = mpf(4)/mpf(8*n+1)-mpf(2)/mpf(8*n+4)-mpf(1)/mpf(8*n+5)-mpf(1)/mpf(8*n+6)
    d = mpf(16.0)**n
    pi += u/d

print(pi)

前 50 位数字正确的是：

3.1415926535 8979323846 2643383279 502884197 1693993751

（已添加空格）