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

Question

对于我的学校项目,我试图计算使用不同方法的价值.我发现的公式之一是可以使用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)))    

Answer 1

这么多条款足以让你超过50个小数位.问题是你将Python浮点数与Decimals混合在一起,因此你的计算会被那些浮点数中的错误所污染,这些错误只能精确到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位数字是垃圾,但其余的都很好.