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

【问题标题】：Why is the computing of the value of pi using the Machin Formula giving a wrong value?为什么使用 Machin 公式计算 pi 的值会给出错误的值？
【发布时间】：2016-09-27 23:47:59
【问题描述】：

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

我用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 次重复时

pi = 3.141592653589793408632493

第一个循环重复 100 次

pi = 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)))

【问题讨论】：

我没有检查你的代码，但是写 Decimal(1/5) 没有问题吗？它不提供 0.2 而是 0.2000....1110223024... 因为 1/5 首先转换为浮点数，而浮点数不能准确存储 0.2。写入 Decimal(1)/Decimal(5) 正好提供 0.2

标签： python math pi

【解决方案1】：

这么多术语足以让您获得超过 50 个小数位。问题是您将 Python 浮点数与小数混合在一起，因此您的计算被这些浮点数中的错误所污染，这些浮点数仅精确到 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) 复制它。

此外，您不需要对一切使用小数。对于count 和a 之类的东西，您可以只使用简单的Python 整数。例如，您对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位是垃圾，但其余的都很好。

【讨论】：