Which floating-point comparison is more accurate, and why?

Question

I'm experimenting with different implementations of the Newton method for calculating square roots. One important decision is when to terminate the algorithm.

Obviously it won't do to use the absolute difference between y*y and x where y is the current estimate of the square root of x, because for large values of x it may not be possible to represent its square root with enough precision.

So I'm supposed to use a relative criteria. Naively I would have used something like this:

static int sqrt_good_enough(float x, float y) {
  return fabsf(y*y - x) / x < EPS;
}

And this appears to work very well. But recently I have started reading Kernighan and Plauger's The Elements of Programming Style and they give a Fortran program for the same algorithm in chapter 1, whose termination criteria, translated in C, would be:

static int sqrt_good_enough(float x, float y) {
  return fabsf(x/y - y) < EPS * y;
}

Both are mathematically equivalent, but is there a reason for preferring one form over the other?

Answer 1

They're still not equivalent; the bottom one is mathematically equivalent to fabsf(y*y - x) / (y*y) < EPS. The problem I see with yours is that if y*y overflows (probably because x is FLT_MAX and y is chosen unluckily), then termination may never occur. The following interaction uses doubles.

>>> import math
>>> x = (2.0 - 2.0 ** -52) * 2.0 ** 1023
>>> y = x / math.sqrt(x)
>>> y * y - x
inf
>>> y == 0.5 * (y + x / y)
True

EDIT: as a comment (now deleted) pointed out, it's also good to share operations between the iteration and the termination test.

EDIT2: both probably have issues with subnormal x. The professionals normalize x to avoid the complications of both extremes.

Answer 2

The two are actually not exactly equivalent mathematically, unless you write fabsf(y*y - x) / (y*y) < EPS for the first one. (sorry for the typo in my original comment)

But I think the key point is to make the expression here match your formula for computing y in the Newton iteration. For example if your y formula is y = (y + x/y) / 2, you should use Kernighan and Plauger's style. If it is y = (y*y + x) / (2*y) you should use (y*y - x) / (y*y) < EPS.

Generally the termination criteria should be that abs(y(n+1) - y(n)) is small enough (i.e. smaller than y(n+1) * EPS). This is why the two expressions should match. If they don't match exactly, it is possible that the termination test decides that the residual is not small enough while the difference in y(n) is smaller than floating point error, due to different scaling. The result would be an infinite loop because y(n) has stopped changing and the termination criteria is never met.

For example the following Matlab code is exactly the same Newton solver as your first example, but it runs forever:

x = 6.800000000000002
yprev = 0
y = 2
while abs(y*y - x) > eps*abs(y*y)
    yprev = y;
    y = 0.5*(y + x/y);
end

The C/C++ version of it has the same problem.