Logo Questions Linux Laravel Mysql Ubuntu Git Menu
 

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

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?

like image 996
lindelof Avatar asked Mar 02 '12 04:03

lindelof


People also ask

What is a better way to compare floating point values?

To compare two floating point or double values, we have to consider the precision in to the comparison. For example, if two numbers are 3.1428 and 3.1415, then they are same up to the precision 0.01, but after that, like 0.001 they are not same.

Why should we never compare floating point values by using exact equality comparison?

Comparing for equality Simple values like 0.1 cannot be precisely represented using binary floating point numbers, and the limited precision of floating point numbers means that slight changes in the order of operations or the precision of intermediates can change the result.

How do you compare two float values?

The compare() method of Float Class is a built-in method in Java that compares the two specified float values. The sign of the integer value returned is the same as that of the integer that would be returned by the function call. Parameters: The function accepts two parameters: f1: The first float value to be compared.

How accurate are floating point numbers?

The data type float has 24 bits of precision. This is equivalent to only about 7 decimal places. (The rest of the 32 bits are used for the sign and size of the number.) The number of places of precision for float is the same no matter what the size of the number.


2 Answers

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.

like image 176
yabba dabba doo Avatar answered Oct 21 '22 22:10

yabba dabba doo


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.

like image 39
fang Avatar answered Oct 21 '22 23:10

fang