floating-pointlanguage-agnosticnumericnumerical-computing

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

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?

Solution

  • 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