交错级数的朴素求和陷阱Naively summing an alternating series

假设你偶然看到了指数函数的幂级数，并决定动手写代码实现它。这是个好主意：你很可能会学到一些东西，尽管可能不是你原本预期的那些。

也许你认为 10^-12 的容差已经足够好了，于是你不断将各项相加，直到下一个要添加的项小于该容差为止。

from math import factorial, exp

def naive_exp(x):
    tolerance = 1e-12
    s = 0
    n = 0
    while True:
        delta = x**n / factorial(n)
        s += delta
        if abs(delta) < tolerance:
            return s
        n += 1

你想测试一下程序，于是通过向函数传入参数 1 来计算自然常数 e。如果将你的结果与调用 exp(1) 得到的结果进行比较，你会发现算出的所有数字完全正确。

现在你尝试计算 exp(-20)。调用 naive_exp(-20) 得到的结果是

    5.47893091802112e-10

但调用 exp(-20) 得到的结果是

    2.061153622438558e-09

不要把这种情况当作偶然现象或者编译器 bug [1] 而敷衍过去。这正是你学习新知识的绝佳机会。

也许你会添加一条 print 语句，来查看存储在变量 s 中的求和中间值。如果这样做，你会观察到部分和在最终稳定之前会发生剧烈的震荡。

这可能看起来是个错误，但如果你更仔细地观察这个级数就会发现：第 n 项是 x^n/n!。由于 x 是负数，各项的符号会交替变化。而且各项的绝对值在变小之前会先变大。当 x = -20 时，每个分子都是前一项的 20 倍，而每个分母只是前一项的 n 倍。因此，这些项会一直增大，直到 n > 20。所以，这种剧烈的震荡是真实存在的，并不是 bug。

绝对值最大的部分和是 21822593.77927747。你知道 exp(-20) 是一个非常小的数字，因此在部分和稳定为一个很小的数之前，必然会发生大量的相消。也许你听说过，这种“相消”正是数值计算丢失精度的原因。如果以前没听说过，那么现在你知道了！

再来看看那个最大的部分和。它的小数点后有 8 位数字。代码打印出的结果已经达到了其所能提供的最高精度，因此此时的误差在 10^-8 数量级。而我们要计算的数字在 10^-9 数量级，如果在这种情况下的计算结果中还有哪怕一位数字是准确的，那纯粹是巧合。

如果你回头用 x = -22 测试你的代码，结果会更糟，竟然计算出了一个负数，而从理论上讲这个值绝不可能为负。但你也能看出原因：你正在让代码去计算一个比代码本身精度还要接近于 0 的数字。

计算机在内部并不是用十进制来表示数字的，但在本例中，上述分析已经足够说明问题了。如果你想深入探究，可以去研究一下浮点数的内部结构。

要解决上述问题，有一个简单的变通方法，但过早发现它会让你错失整个学习过程。你可以将 exp(-20) 计算为 1/exp(20)，这样就能避免所有的相消问题，因为 exp(20) 的级数并不存在符号交替的情况。

[1] 编译器偶尔确实会存在 bug，但更有可能（且概率高出好几个数量级）的是你的代码写错了。