Numerically stable evaluation of sqrt(x+a) - sqrt(x)

Question

Is there an elegant way of numerically stable evaluating the following expression for the full parameter range x,a >= 0?

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

Also is there any programming language or library that does provide this kind of function? If yes, under what name? I have no specific problem using the above expression right now, but encountered it many times in the past and always thought that this problem must have been solved before!

Answer 1

Yes, there is! Provided that at least one of x and a is positive, you can use:

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

which is perfectly numerically stable, but hardly worth a library function in its own right. Of course, when x = a = 0, the result should be 0.

Explanation: sqrt(x + a) - sqrt(x) is equal to (sqrt(x + a) - sqrt(x)) * (sqrt(x + a) + sqrt(x)) / (sqrt(x + a) + sqrt(x)). Now multiply the first two terms to get sqrt(x+a)^2 - sqrt(x)^2, which simplifies to a.

Here's an example demonstrating the stability: the troublesome case for the original expression is where x + a and x are very close in value (or equivalently when a is much smaller in magnitude than x). For example, if x = 1 and a is small, we know from a Taylor expansion around 1 that sqrt(1 + a) should be 1 + a/2 - a^2/8 + O(a^3), so sqrt(1 + a) - sqrt(1) should be close to a/2 - a^2/8. Let's try that for a particular choice of small a. Here's the original function (written in Python, in this case, but you can treat it as pseudocode):

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

and here's the stable version:

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

Now let's see what we get with x = 1 and a = 2e-10:

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

The value we should have got is (up to machine accuracy): a/2 - a^2/8 - for this particular a, the cubic and higher order terms are insignificant in the context of IEEE 754 double-precision floats, which only provide around 16 decimal digits of precision. Let's compute that value for comparison:

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