views:

235

answers:

2

To make the problem short let's say I want to compute expression: a / (b - c) on float's.

To make sure the result is meaningful, I can check if 'b' and 'c' are inequal:

float EPS = std::numeric_limits<float>::epsilon();
if ((b - c) > EPS || (c - b) > EPS)
{
    return a / (b - c);
}

but my tests show it is not enough to guarantee either meaningful results nor not failing to provide a result if it is possible.

Case 1: a = 1.0f; b = 0.00000003f; c = 0.00000002f;

Result: The if condition is NOT met, but the expression would produce a correct result 100000008 (as for the floats' precision).

Case 2: a = 1e33f; b = 0.000003; c = 0.000002;

Result: The if condition is met, but the expression produces not a meaningful result +1.#INF00.

I found it much more reliable to check the result, not the arguments:

const float INF = numeric_limits<float>::infinity();
float x = a / (b - c);
if (-INF < x && x < INF)
{
     return x;
}

But what for is the epsilon then and why is everyone saying epsilon is good to use?

+1  A: 
Marcelo Cantos
Pascal Cuoq
Thanks for the comment, Pascal. I could wrap it all in a union, but that would be a bit too much verbosity for such a tangential point.
Marcelo Cantos
Mark Dickinson
Mark: yes, there are edge cases, but yours isn't one of them. Positive and negative numbers are miles apart due to the sign bit (as are the zeros from each other and from non-zeros). I believe that denorms are a real edge case.
Marcelo Cantos
Marcelo: I'm confused. In my case, the two values are exactly 1 ulp apart: b is the smallest positive subnormal value, while c is -0.0, so by any reasonable measure they should be considered close. How is this not an edge case? The problem with the sign bit is exactly what I was referring to. The LHS of your comparison evaluates to 2147483647 (=2**31-1) in this case.And no, subnormals aren't a problem for this sort of comparison *unless* the signs differ.
Mark Dickinson
Ah, now I understand. You _want_ the two numbers to be considered nearly equal. It really depends in whether you are considering similarity in an arithmetic or geometric sense. For arithmetic similarity, you are quite right, but the ulp measure is only suitable for geometric similarity, in which case your two numbers are in entirely different universes, for all intents and purposes.
Marcelo Cantos
+7  A: 

"you must use an epsilon when dealing with floats" is a knee-jerk reaction of programmers with a superficial understanding of floating-point computations, for comparisons in general (not only to zero).

This is usually unhelpful because it doesn't tell you how to minimize the propagation of rounding errors, it doesn't tell you how to avoid cancellation or absorption problems, and even when your problem is indeed related to the comparison of two floats, it doesn't tell you what value of epsilon is right for what you are doing.

If you have not read What Every Computer Scientist Should Know About Floating-Point Arithmetic, it's a good starting point. Further than that, if you are interested in the precision of the result of the division in your example, you have to estimate how imprecise b-c was made by previous rounding errors, because indeed if b-c is small, a small absolute error corresponds to a large absolute error on the result. If your concern is only that the division should not overflow, then your test (on the result) is right. There is no reason to test for a null divisor with floating points, you just test for overflow of the result.

Regarding the propagation of rounding errors, there exists specialized analyzers that can help you estimate it, because it is a tedious thing to do by hand.

Pascal Cuoq