Floating Point; Division vs Multiplication

Question

You have 2 functions;

f(x)= x(((x+1)^(1/2))-(x^(1/2)))    
g(x)= x/(((x+1)^(1/2))+(x^(1/2)))

Which one is more accurate?

Side note: If you could explain why, that would really help me out, I have a feeling it's f(x) as there isn't a denominator but I'm not 100% certain.

Answer 1

Floating-point arithmetic is quite different from real mathematics. In particular, floating-point arithmetic is not associative or distributive. Therefore mathematically equivalent expressions are not necessarily equivalent when evaluated with finite-precision floating-point arithmetic. This is the case even when individual operations, such as addition, subtraction, multiplication, division, and square root produce correctly rounded results, as required by the IEEE-754 floating-point standard.

In this case g() will be much more accurate than f() on average, and also for the specific case of x = 500. The reason is that f() suffers from subtractive cancellation. This occurs during effective subtraction of two floating-point numbers that are almost identical in magnitude. The leading digits cancel during the subtraction, leaving only a few remaining trailing digits that enter into subsequent computation. Also, any rounding error accumulated in the trailing digits of the original operands being subtracted may be magnified by subsequent computation. An expanded explanation with example can be found in this Wikipedia article.

In this case, sqrt(x+1) and sqrt(x) are almost the same magnitude, in particular as the magnitude of x increases. Using the example of x = 500 and employing IEEE-754 single-precision arithmetic, we find:

x = 500   f(x) = 0x1.659ae0p+3 (11.175156)   reference = 0x1.659798p+3 (11.174755)
x = 500   g(x) = 0x1.659798p+3 (11.174755)   reference = 0x1.659798p+3 (11.174755)

The error in f(500) is 420 ulps, while g(500) delivers the correctly rounded single-precision result.