数学稳定评估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