21

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

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

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?

Community
  • 1
  • 1
Michal Czardybon
  • 2,795
  • 4
  • 28
  • 40

2 Answers2

24

"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-point numbers, you just test for overflow of the result, which captures both the cases where the divisor is null and where the divisor is so small as to make the result not representable with any precision.

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
  • 79,187
  • 7
  • 161
  • 281
  • 5
    "it doesn't tell you what value of epsilon is right for what you are doing": Please emphasize this sentence more, it is really important. – Alexandre C. Aug 09 '11 at 09:43
  • *There is no reason to test for a null divisor with floating-point numbers, you just test for overflow of the result.* Floating-point division-by-zero results in undefined behaviour in C and C++. – Max Barraclough Mar 29 '21 at 10:06
  • @MaxBarraclough The C standard does not mandate a floating-point system with infinities or NaN, so of course in the C standard floating-point division by zero is UB. Floating-point overflow is also UB in this minimal setting per C17 6.5:5. Again there is no need to single out division by zero. Now, common sense dictates that when a C compiler documents that it implements IEEE 754 for `float` and `double`, overflowing operations produce inf/NaN according to the rules of IEEE 754, and the rules of IEEE 754 also apply to division by zero. This is what “implementing IEEE 754” means. – Pascal Cuoq Mar 29 '21 at 19:51
  • @MaxBarraclough A compiler cannot claim to implement IEEE 754 and leave `0./0.` UB. Words have meaning. – Pascal Cuoq Mar 29 '21 at 19:52
3

Epsilon is used to determine whether two numbers subject to rounding error are close enough to be considered "equal". Note that it is better to test fabs(b/c - 1) < EPS than fabs(b-c) < EPS, and even better — thanks to the design of IEEE floats — to test abs(*(int*)&b - *(int*)&c) < EPSI (where EPSI is some small integer).

Your problem is of a different nature, and probably warrants testing the result rather than the inputs.

Marcelo Cantos
  • 181,030
  • 38
  • 327
  • 365
  • 2
    The "strict aliasing rules" crowd would have something to say about your `*(int*)&` construct. I personally do not agree with the attitude of compilers which emphasize benchmarks timings over existing practices and programmers' widespread understanding of the meaning of the language, but they have a point: your code will be miscompiled by modern compilers. – Pascal Cuoq Apr 28 '10 at 13:33
  • 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 Apr 28 '10 at 13:44
  • Careful with the `*(int*)&` trick: apart from the type punning problem Pascal mentions, it also won't work near 0: if b and c are very close but have opposite signs (e.g. b = 2e-45f, c = -0.0f) then your test will fail. – Mark Dickinson Apr 28 '10 at 13:45
  • 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 Apr 28 '10 at 14:02
  • 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 Apr 28 '10 at 14:21
  • 2
    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 Apr 28 '10 at 21:57
  • @MarceloCantos A `union` won't help there; it's only in C that type-punning through a `union` is legal, but in C++ it's not and is undefined behaviour, and context indicates we're talking about the latter. – underscore_d Sep 20 '17 at 13:50