python - 为什么使用 Machin 公式计算 pi 的值会给出错误的值？

Question

对于我的学校项目，我试图计算使用不同方法的值(value)。我找到的公式之一是 Machin 公式，它可以使用 arctan(x) 的泰勒展开来计算。

我用python写了下面的代码:

import decimal

count = pi = a = b = c = d = val1 = val2 = decimal.Decimal(0) #Initializing the variables      
decimal.getcontext().prec = 25 #Setting percision

while (decimal.Decimal(count) <= decimal.Decimal(100)): 
    a = pow(decimal.Decimal(-1), decimal.Decimal(count))
    b = ((decimal.Decimal(2) * decimal.Decimal(count)) + decimal.Decimal(1))
    c = pow(decimal.Decimal(1/5), decimal.Decimal(b))
    d = (decimal.Decimal(a) / decimal.Decimal(b)) * decimal.Decimal(c)
    val1 = decimal.Decimal(val1) + decimal.Decimal(d)
    count = decimal.Decimal(count) + decimal.Decimal(1)
    #The series has been divided into multiple small parts to reduce confusion

count = a = b = c = d = decimal.Decimal(0) #Resetting the variables

while (decimal.Decimal(count) <= decimal.Decimal(10)):
    a = pow(decimal.Decimal(-1), decimal.Decimal(count))
    b = ((decimal.Decimal(2) * decimal.Decimal(count)) + decimal.Decimal(1))
    c = pow(decimal.Decimal(1/239), decimal.Decimal(b))
    d = (decimal.Decimal(a) / decimal.Decimal(b)) * decimal.Decimal(c)
    val2 = decimal.Decimal(val2) + decimal.Decimal(d)
    count = decimal.Decimal(count) + decimal.Decimal(1)
    #The series has been divided into multiple small parts to reduce confusion

pi = (decimal.Decimal(16) * decimal.Decimal(val1)) - (decimal.Decimal(4) * decimal.Decimal(val2))
print(pi)

问题是，无论循环重复多少次，我都只能得到 pi 的正确值，直到小数点后 15 位。

例如:

第一个循环重复 11 次

圆周率 = 3.141592653589793408632493

第一个循环重复 100 次

圆周率 = 3.141592653589793410703296

我没有增加第二个循环的重复次数，因为 arctan(1/239) 非常小，重复几次就达到了一个极小的值，因此应该不会影响仅小数点后 15 位的 pi 值。

额外信息:

Machin 公式指出:

   π = (16 * Summation of (((-1)^n) / 2n+1) * ((1/5)^(2n+1))) - (4 * Summation of (((-1)^n) / 2n+1) * ((1/239)^(2n+1)))    

Answer 1

那么多术语足以让您保留超过 50 位小数。问题是您将 Python float 与小数混合在一起，因此您的计算被这些 float 中的错误所污染，这些 float 仅精确到 53 位(大约 15 位十进制数字)。

你可以通过改变来解决这个问题

c = pow(decimal.Decimal(1/5), decimal.Decimal(b))

到

c = pow(1 / decimal.Decimal(5), decimal.Decimal(b))

或

c = pow(decimal.Decimal(5), decimal.Decimal(-b))

显然，需要对

进行类似的更改

c = pow(decimal.Decimal(1/239), decimal.Decimal(b))

---

您可以使您的代码更易读很多。对于初学者，您应该将计算 arctan 级数的东西放入一个函数中，而不是为 arctan(1/5) 和 arctan(1/239) 复制它。

此外，您不需要对所有 使用 Decimal。您可以只使用简单的 Python 整数来表示 count 和 a 之类的东西。例如，您对 a 的计算可以写成

a = (-1) ** count

或者您可以在循环外将 a 设置为 1 并在每次循环中取反。

这是您的代码的更紧凑版本。

import decimal

decimal.getcontext().prec = 60 #Setting precision

def arccot(n, terms):
    base = 1 / decimal.Decimal(n)
    result = 0
    sign = 1
    for b in range(1, 2*terms, 2):
        result += sign * (base ** b) / b
        sign = -sign
    return result

pi = 16 * arccot(5, 50) - 4 * arccot(239, 11)
print(pi)

输出

3.14159265358979323846264338327950288419716939937510582094048

最后 4 位是垃圾，但其余的都很好。