Python 对大浮点数计算错误

Question

我的算法 returns 圆周率的近似值。它通过使用等式

来工作

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

可能不是 python 本身有这个问题，但是在 2^22 你运行遇到了 python 内部使用的类型的数值问题。

对于python使用的64位浮点数，你得到52位的尾数；将您的准确性限制在逗号后面的大约 16（十进制）数字。

编辑：我真的不知道回到 Matlab 有什么帮助——默认情况下，Matlab 内部通常也使用相同的 64 位浮点类型。

Answer 2

您正在使用的迭代公式有些问题：随着迭代的进行，数量 n 变得比 an 大得多。所以在表达式

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

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

短篇小说：要使用最初编写的公式获得 pi 的 d 位数，您将需要在中间处理更多的 d 位数计算。通过一些工作，您可以证明您需要在内部使用多于 2d 位的精度才能获得 pi 的 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 浮点数就可以获得最多 15 位的精度。事实上，运行你的代码在这个新的迭代中，我在 30 次迭代后得到了 3.1415926535897927 的值，这只是通过单个 ulp（最后一个单位）对 pi 的最佳双精度近似地方）。

要获得比这更多的数字，我们需要使用 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 的精度。

Python 对大浮点数计算错误

Python miscalculating for large floats

python

algorithm

math

floating-point