comparing doubles for equality

[John Smith]

This code for the comparison of fp types is taken from the C FAQ.
Any problems using it in a macro?

/* compare 2 doubles for equality */
#define DBL_ISEQUAL(a,b) (fabs((a)-(b))<=(DBL_EPSILON)*fabs((a)))

Do the same issues involved in comparing 2 fp types for equality
apply to comparing a float to zero? E.g. is if(x == 0.0)
considered harmful?

 

[Tim Prince]

This code for the comparison of fp types is taken from the C FAQ.
Any problems using it in a macro?

/* compare 2 doubles for equality */
#define DBL_ISEQUAL(a,b) (fabs((a)-(b))<=(DBL_EPSILON)*fabs((a)))

Do the same issues involved in comparing 2 fp types for equality
apply to comparing a float to zero? E.g. is if(x == 0.0) considered
harmful?

Depending on how you use this, you could make your application highly
dependent on external issues, such as whether you compile for an extra
precision mode, or disable gradual underflow.

 

[Dik T. Winter]

> This code for the comparison of fp types is taken from the C FAQ.
> Any problems using it in a macro?
>
> /* compare 2 doubles for equality */
> #define DBL_ISEQUAL(a,b) (fabs((a)-(b))<=(DBL_EPSILON)*fabs((a)))
>
> Do the same issues involved in comparing 2 fp types for equality
> apply to comparing a float to zero? E.g. is if(x == 0.0)
> considered harmful?

If you insert 0.0 for 'a' you will see that the two give precisely the
same result. On the other hand, I think it is possible to construct
two floating point numbers 'a' and 'b', where the result is asymmetric.
But let that not deter you from using the macro, when it is asymmetric
you are in the outskirts of floating-point arithmetic.

 

[Thad Smith]

This code for the comparison of fp types is taken from the C FAQ.
Any problems using it in a macro?

/* compare 2 doubles for equality */
#define DBL_ISEQUAL(a,b) (fabs((a)-(b))<=(DBL_EPSILON)*fabs((a)))

This construction is misleading and I would never use it, because the
implied function, determining whether two doubles are equal, is not an
accurate description of the returned value.

Do the same issues involved in comparing 2 fp types for equality
apply to comparing a float to zero? E.g. is if(x == 0.0) considered
harmful?

Which issues are those? The test "if (x == 0.0)", in contrast, does not
have the inaccurate description that the DBL_ISEQUAL macro does.

In both cases, you should employ good analysis for floating point
comparisons.

 

[CBFalconer]

Depending on how you use this, you could make your application
highly dependent on external issues, such as whether you compile
for an extra precision mode, or disable gradual underflow.

That's why gcc has the Wfloat-equal switch. As far as the macro is
concerned, the a argument better not have any side effects.

 

[CBFalconer]

This construction is misleading and I would never use it, because
the implied function, determining whether two doubles are equal,
is not an accurate description of the returned value.

How about:

#define DUNEQUAL(a, b) (fabs((a)-(b)) > (DBL_EPSILON)*fabs((a)))

with a caveat against passing an a with side effects.

 

[Thad Smith]

How about:

#define DUNEQUAL(a, b) (fabs((a)-(b)) > (DBL_EPSILON)*fabs((a)))

with a caveat against passing an a with side effects.

That misses the point. The only true equality test is given by a==b.
If you need an epsilon, fine, but you should make it explicit.
DBL_EPSILON isn't necessarily the correct choice for a particular
application. In fact, the two tests above may not provide any advantage
over a strict equality.

 

[CBFalconer]

That misses the point. The only true equality test is given by
a==b. If you need an epsilon, fine, but you should make it
explicit. DBL_EPSILON isn't necessarily the correct choice for a
particular application. In fact, the two tests above may not
provide any advantage over a strict equality.

No, I think you miss the point. Floating point is inherently an
approximation, and the above says that anything within the
resolution (i.e. the approximation) is to be considered equal. Any
time you see (a == b) for floating point operands, it is probably a
bug. For example:

for (a = 1.0; a != 0; a -= 0.1) dosomethingwith(a);

will probably not behave. While:

for (a = 1.0; !DUNEQUAL(a, 0.0); a -= 0.1) dosomethingwith(a);

will behave. Of course:

for (a = 1.0; a > 0; a -= 0.1) dosomethingwith(a);

may behave. But it is a crapshoot whether it does an extra cycle.

 

[SM Ryan]

# This code for the comparison of fp types is taken from the C FAQ.
# Any problems using it in a macro?
#
# /* compare 2 doubles for equality */
# #define DBL_ISEQUAL(a,b) (fabs((a)-(b))<=(DBL_EPSILON)*fabs((a)))
#
# Do the same issues involved in comparing 2 fp types for equality
# apply to comparing a float to zero? E.g. is if(x == 0.0)
# considered harmful?

Floats and doubles use a small number of bits for their values,
so exact equality to zero or any other representable value is
meaningful. Two problems are

(1) most machines do not use radix 10 representation; some decimal
values, such as 0.1 can only be approximated. That means if you
write (x==0.1) you're not really testing x against one tenth but
some other value that is close to one tenth. So the test is not
what you think it is.

(2) many computations with reals are iterative approximations to
an unknown value. Because of the precision limits, they often
cannot reach the exact correct value, so you have to accept when
it is close enough.

 

[Dik T. Winter]

....
> For example:
> for (a = 1.0; a != 0; a -= 0.1) dosomethingwith(a);
> will probably not behave. While:
> for (a = 1.0; !DUNEQUAL(a, 0.0); a -= 0.1) dosomethingwith(a);

you wish to omit the exclamation mark here.

> will behave. Of course:
> for (a = 1.0; a > 0; a -= 0.1) dosomethingwith(a);
> may behave. But it is a crapshoot whether it does an extra cycle.

And now compare the three with:
for (a = 1.0; DUNEQUAL(0.0, a); a -= 0.1) dosomethingwith(a);

A better test would be:

#define DUNEQUAL(a, b) (fabs((a) - (b)) > \
(fabs(a) > fabs(b) ? DBL_EPSILON * fabs(a) : DBL_EPSILON * fabs(b)))

in this way there is no asymmetry.

 

[Thad Smith]

No, I think you miss the point. Floating point is inherently an
approximation, and the above says that anything within the
resolution (i.e. the approximation) is to be considered equal. Any
time you see (a == b) for floating point operands, it is probably a
bug. For example:

for (a = 1.0; a != 0; a -= 0.1) dosomethingwith(a);

While I agree that this is likely to fail ...

for (a = 1.0; !DUNEQUAL(a, 0.0); a -= 0.1) dosomethingwith(a);

.... there is no guarantee that this second form is correct, for the
definition supplied above. Why? Because each subtraction may invoke a
roundoff. After several roundings, you maybe off by more than one lsb
for a number with the magnitude of 1, let alone, say, one lsb of 1e-7.

Note also that DUNEQUAL, is not commutative, which is fails the
principle of least astonishment.

for (a = 1.0; a > 0; a -= 0.1) dosomethingwith(a);

may behave. But it is a crapshoot whether it does an extra cycle.

And, to make it robust, you could do

for (a = 1.0; a > 0.1/2; a -= 0.1) dosomethingwith(a);

or, to use the best approximation of of each nominal value of a:

int i;
for (i = 10; i > 0; i--) { dosomethingwith (i*0.1); }

or, for the pragmatic types

for (a = 10; a > 0; a -= 1) dosomethingwith (a/10);

The last form is not guaranteed, but performs correctly on any floating
point system I am familiar with, except maybe logarithmic systems, due
to the limited number of bits required to exactly represent the integer
values.

 

[CBFalconer]

you wish to omit the exclamation mark here.


And now compare the three with:
for (a = 1.0; DUNEQUAL(0.0, a); a -= 0.1) dosomethingwith(a);

A better test would be:

#define DUNEQUAL(a, b) (fabs((a) - (b)) > \
(fabs(a) > fabs(b) ? DBL_EPSILON * fabs(a) : DBL_EPSILON * fabs(b)))

in this way there is no asymmetry.

Not needed. The use of epsilon is only needed when a and b are
nearly equal, thus either can be used to derive it. The problem of
a zero remains, and ties in with the use of extremely small
operands in general. To handle all these cases we probably should
use a function rather than a macro.

 

[pete]

No, I think you miss the point. Floating point is inherently an
approximation, and the above says that anything within the
resolution (i.e. the approximation) is to be considered equal. Any
time you see (a == b) for floating point operands, it is probably a
bug.

I disagree.
It really all depends
on what you're doing with the floating point values.
If you're using qsort to sort an array of doubles,
then there's really no place for DBL_EPSILON in your compar function.

int compar(const void *arg1, const void *arg2)
{
return *(double *)arg2 > *(double *)arg1 ? -1
: *(double *)arg2 != *(double *)arg1;
}