Python 计算错误的大 float

Question

我的算法返回 pi 的近似值。它通过使用等式来工作

a_2n= math.sqrt(-2*(math.sqrt((n**2)-4*(a_n**2))-n)*n)/2

其中 a_n 是当 n 是某个数字时的近似值，而 a_2n 是当 n 是它的两倍时的近似值。它以 n = 4 和 a_n 为 2 开始，并应用该公式直到 n 足够高。 n越大计算越准确，只是当n>=2^22时突然停止收敛到pi

完整代码如下:

import math
def pi(n):
    value = 2
    round = 4
    loop = 2
    while round != n:
        value= a2n(value,round)
        round = round*2
        loop +=1
        print("n=2^", loop, " pi = ", value)
    return value


def a2n(an,n):
    return  math.sqrt(-2*(math.sqrt((n**2)-4*(an**2))-n)*n)/2


print(pi(2**25))

我相信数学没问题，所以我认为 python 在处理较大的数字时遇到了问题。它从“3.141596553704”变为“3.14167426502175”，并从那里变得更糟。

Answer 1

您使用的迭代公式有点病态:随着迭代的进行，数量 n 变得比 an 大得多。所以在表达式中

math.sqrt(n**2-4*an**2)-n

平方根的结果将接近n，因此外减法是两个几乎相等(相对意义上)的量的减法。现在，如果您正在使用具有 16 位十进制精度的常规 Python float 进行计算，那么随着迭代的进行，该减法将为您提供一个仅精确到少数数字的结果。请参阅 loss of significance 上的维基百科页面以获得对此问题的更一般性讨论。

短篇小说:要使用您最初编写的公式得到 d 位 pi，您需要处理的远不止 d 中间计算中的数字。通过一些工作，您可以证明您将需要在内部处理多于 2d 位的精度才能获得 d 位准确数字圆周率。即便如此，您也必须小心，只使用您需要的迭代次数:迭代次数过多，无论您使用多少中间精度，精度都会再次下降。

但是这里有比中间精度加倍更好的替代方法，那就是重写您的公式，以便首先避免失去重要性。如果将 sqrt(n**2-4*an**2)-n 乘以共轭表达式 sqrt(n**2-4*an**2)+n，你得到 -4*an**2。所以原来的差值可以改写为-4*an**2/(sqrt(n**2-4*an**2)+n)。将其插入您的原始公式并稍微简化会导致如下所示的迭代步骤:

def a2n(an, n):
    return an*math.sqrt(2*n/(math.sqrt(n*n-4*an*an)+n))

从数学上讲，这与您的 a2n 函数进行的计算完全相同，但从数字的角度来看，它的表现要好得多。

如果您使用此迭代步骤代替原始步骤，您会看到更小的舍入误差，并且您应该能够仅使用 Python float 获得高达 15 位的精度。事实上，用这个新的迭代运行你的代码，我在 30 次迭代后得到 3.1415926535897927 的值，这只是通过单个 ulp(单位在最后一个地方)对 pi 的最佳 double 近似).

要获得比这更多的数字，我们需要使用 decimal模块。这是一个片段，基于您的代码，但使用我建议的迭代修改计算，使用 decimal 模块获取 pi 的值，该值精确到 51 位有效数字。它使用 55 位有效数字的内部精度，以允许累积舍入误差。

from decimal import getcontext

context = getcontext()
context.prec = 55    # Use 55 significant digits for all operations.
sqrt = context.sqrt  # Will work for both ints and Decimal objects,
                     # returning a Decimal result.

def step(an, n):
    return an*sqrt(2*n/(sqrt(n*n-4*an*an)+n)), n*2

def compute_pi(iterations):
    value, round = 2, 4
    for k in range(iterations):
        value, round = step(value, round)
    return value

pi_approx = compute_pi(100)
print("pi = {:.50f}".format(pi_approx))

结果如下:

pi = 3.14159265358979323846264338327950288419716939937511

使用原始公式，中间计算需要超过 100 的精度。