Accurate sqrt(1 + (x/2)^2) + x/2

Question

I need to compute sqrt(1 + (x/2)^2) + x/2 numerically, for positive x. Using this expression directly fails for very large values of x. How can I rewrite it to obtain a more accurate evaluation?

Answer 1

For very large x you can factor out an x/2:

sqrt(1 + (x/2)^2) + x/2
 = (x/2) * sqrt( 1/(x/2)^2 + (x/2)^2/(x/2)^2) + x/2
 = (x/2) * sqrt( (2/x)^2 + 1 ) + x/2

For x > 2/sqrt(eps) the square root will actually evaluate to 1 and your whole expression will simplify to just x. Assuming you need to cover the entire range [0, infinity], I would suggest just branching at that point and return x in this case and your original formula, otherwise:

if x > 2/sqrt(eps)  // eps is the machine epsilon of your float type
    return x
else
    return sqrt(1 + (x/2)^2) + x/2

Answer 2

Many programming languages offer a function hypot(x,y) that computes sqrt (x*x + y*y) while avoiding overflow and underflow in intermediate computation. Many implementations of hypot also compute the result more accurately than the naive expression. These advantages come at the expense of a moderate increase in run time.

With this function, the given expression can be written as hypot (1.0, 0.5*x) + 0.5*x. If your programming language of choice does not support hypot or an equivalent function, you may be able to adapt the implementation I provided in this answer.