sqrt(x+a) - sqrt(x) 的数值稳定评估

Question

对于完整的参数范围 x,a >= 0，是否有一种优雅的数值稳定评估以下表达式的方法？

f(x,a) = sqrt(x+a) - sqrt(x)

还有任何编程语言或库可以提供这种功能吗？如果是，以什么名字命名？我现在使用上面的表达式没有什么具体问题，但是以前遇到过很多次，一直以为这个问题之前一定已经解决了！

Answer 1

是的，有！如果x 和a 至少有一个是肯定的，则可以使用：

f(x, a) = a / (sqrt(x + a) + sqrt(x))

它在数值上是完全稳定的，但它本身并不值得一个库函数。当然x = a = 0时，结果应该是0。

解释：sqrt(x + a) - sqrt(x) 等于 (sqrt(x + a) - sqrt(x)) * (sqrt(x + a) + sqrt(x)) / (sqrt(x + a) + sqrt(x))。现在将前两项相乘得到sqrt(x+a)^2 - sqrt(x)^2，简化为a。

这是一个展示稳定性的示例：原始表达式的麻烦情况是x + a 和x 的值非常接近（或者等效地，当a 的大小比x 小得多时）。例如，如果x = 1 和a 很小，我们从1 周围的泰勒展开式知道sqrt(1 + a) 应该是1 + a/2 - a^2/8 + O(a^3)，所以sqrt(1 + a) - sqrt(1) 应该接近a/2 - a^2/8。让我们尝试选择小的a。这是原始函数（在本例中是用 Python 编写的，但您可以将其视为伪代码）：

def f(x, a):
    return sqrt(x + a) - sqrt(x)

这是稳定版：

def g(x, a):
    if a == 0:
        return 0.0
    else:
        return a / ((sqrt(x + a) + sqrt(x))

现在让我们看看x = 1 和a = 2e-10 得到了什么：

>>> a = 2e-10
>>> f(1, a)
1.000000082740371e-10
>>> g(1, a)
9.999999999500001e-11

我们应该得到的值是（取决于机器精度）：a/2 - a^2/8 - 对于这个特定的a，三次和更高阶项在 IEEE 754 双精度浮点数的上下文中是微不足道的，它只提供大约 16 位小数的精度。让我们计算该值以进行比较：

>>> a/2 - a**2/8
9.999999999500001e-11