Help, best way to compare doubles
Author Message Help, best way to compare doubles

[ Followup-To: comp.lang.c.moderated ]

Quote:

>: Having a hard time to figure out a function on how to compare doubles.
>: Any ideas?
>No function necessary.. C has this really cool boolean equality operator
>that does this sort of thing for you.  Have a look:

[...]

It is unwise to compare floating point numbers for equality without
thinking hard first.  Roundoff error means that mathematically
equivalent expressions may yield different results for different
orderings; for example, a+b+c may very well compare unequal to
a+c+b (try this out with a=1.0, b=-1.0 and c=1.E-20, for example).
Roundoff error can also mean that algorithms which should converge
don't, so beware.  (I even know one system where a/2*2 != a for one
quarter of floating point numbers - no points for guessing it's
a /370 :-)

The correct way to compare two floating point numbers for "equality" to
use fabs(a-b) < eps.  The problem arises in determining what to use for
eps, of course :-).  What's reasonable depends very much on your problem
and on the floating point characteristics of your system (for example,
do your values span many orders of magnitude?  Are they close to unity?
Could they have opposing signs?), so it's really impossible to tell
without a fairly thorough analysis.
--

The joy of engineering is to find a straight line on a double
logarithmic diagram.

Wed, 01 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

>[ Followup-To: comp.lang.c.moderated ]

>>: Having a hard time to figure out a function on how to compare doubles.
>>: Any ideas?
...
>The correct way to compare two floating point numbers for "equality" to
>use fabs(a-b) < eps..

Or to use _normalized_ differences.

First version (has a division and a square root):
fabs( a - b )/sqrt( a*a + b*b ) < eps

better (has no division):
fabs( a - b ) < eps * sqrt( a*a + b*b )

better yet (neither division nor square root):
( t = a - b )*t < eps * ( a*a + b*b )

Thu, 02 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

>>: Having a hard time to figure out a function on how to compare doubles.
>>: Any ideas?
>The correct way to compare two floating point numbers for "equality" to
>use fabs(a-b) < eps.  The problem arises in determining what to use for
>eps, of course :-).

Sometimes I do this:

if(tmp1= a+b, tmp2= (a-b)/64, tmp2+= tmp1, tmp1==tmp2) {
} else {
/* different */;
}

The idea is to test whether (a-b) falls off the mantissa when scaled to
the same exponent as (a+b).  The 64 is for slop.  Hopefully this code
won't be placed in the middle of a tight loop...

I believe the C compiler is forbidden to optimize the code into a
mathematically equivalent but numerically incorrect test.  Can it be
done without the temporaries ?  I don't read the standard anymore :-)
--

Thu, 02 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

>>[ Followup-To: comp.lang.c.moderated ]

>>>: Having a hard time to figure out a function on how to compare doubles.
>>>: Any ideas?
>...
>>The correct way to compare two floating point numbers for "equality" to
>>use fabs(a-b) < eps..
>Or to use _normalized_ differences.
>First version (has a division and a square root):
>    fabs( a - b )/sqrt( a*a + b*b ) < eps
>better (has no division):
>    fabs( a - b ) < eps * sqrt( a*a + b*b )
>better yet (neither division nor square root):
>   ( t = a - b )*t < eps * ( a*a + b*b )

Nice.  Is there a general method you are using here?  Details
would be nice.  If you couch it in terms of C, the moderator may
whistle and look the other way <G>.
Free lessons on whistling and averting ones eyes available here.

Sincerely,

Gene Wirchenko

C Pronunciation Guide:
y=x++;     "wye equals ex plus plus semicolon"
x=x++;     "ex equals ex doublecross semicolon"

Fri, 03 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

> Or to use _normalized_ differences.
> First version (has a division and a square root):
>     fabs( a - b )/sqrt( a*a + b*b ) < eps
> better (has no division):
>     fabs( a - b ) < eps * sqrt( a*a + b*b )
> better yet (neither division nor square root):
>    ( t = a - b )*t < eps * ( a*a + b*b )

^

Does any rule in C++ tell you that the latter t will be evaluated later
than the expression "t = a - b"?
[C++?  Oh, well.  I don't think there is one in C, anyway. -mod]

--
/-----------------------------------v------------------------------------\
|Isaac To, kar-keung                |  'y.a1j                            |
|MPhil, Computer Science            |  -p:b>w>G  -u>G:S\$h=R5{            |

|URL: http://www.cs.hku.hk/~kkto    |  :t-6!Ghttp://www.cs.hku.hk/~kkto  |
\-----------------------------------^------------------------------------/

Fri, 03 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

>>The correct way to compare two floating point numbers for "equality" to
>>use fabs(a-b) < eps..

>Or to use _normalized_ differences.
>First version (has a division and a square root):
>    fabs( a - b )/sqrt( a*a + b*b ) < eps

This only works correctly if a*a + b*b does not overflow your machine
number, and also that a*a+b*b does not underflow.  Still some thought
required...

A cleaner overall solution would be to use interval arithmetic instead
of conventional floating point variables.  But by then, we'd probably
be outside of what is normally considered the scope of C.
--

The joy of engineering is to find a straight line on a double
logarithmic diagram.

Fri, 03 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

> The correct way to compare two floating point numbers for "equality" to
> use fabs(a-b) < eps.  The problem arises in determining what to use for
> eps, of course :-).  What's reasonable depends very much on your problem
> and on the floating point characteristics of your system (for example,
> do your values span many orders of magnitude?  Are they close to unity?

I often use macros like:

#define TOLERENCE 1.0e6
#define EQUAL( a, b ) (fabs(a-b)<=(fabs(a)/TOLERENCE))

By relating epsilon  to the values used, the macro becomes more general.
Note that the "<=" is important for the case a==b==0.

--
#####################################################################

Emmenjay Consulting
Applications Support              PC Installation/Maintenance
Training                                    Technical Writing
Windows/Unix               General Consulting/Troubleshooting

PO Box 909                                             Ph 018 240 704
Kensington 2033
AUSTRALIA
#####################################################################

Fri, 03 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

|>
|> > Or to use _normalized_ differences.
|>
|> > First version (has a division and a square root):
|> >     fabs( a - b )/sqrt( a*a + b*b ) < eps
|>
|> > better (has no division):
|> >     fabs( a - b ) < eps * sqrt( a*a + b*b )
|>
|> > better yet (neither division nor square root):
|> >    ( t = a - b )*t < eps * ( a*a + b*b )
|>                    ^
|>
|> Does any rule in C++ tell you that the latter t will be evaluated later
|> than the expression "t = a - b"?
|> [C++?  Oh, well.  I don't think there is one in C, anyway. -mod]
|>
While we're at it, note that the final version (even after *that* bug
is fixed) is _not_ equivalent to the ones which precede it. It *should*
be something like:

t = a - b;

if ( t*t < eps*eps*(a*a + b*b) ) {
/* 'a' and 'b' are "morally equal" */
} else {
/* they're not */
}

For small values of 'eps' (which are surely the most interesting), the
original (flawed) approach is far too tolerant ...

|> --
|> /-----------------------------------v------------------------------------\
|> |Isaac To, kar-keung                |  'y.a1j                            |
|> |MPhil, Computer Science            |  -p:b>w>G  -u>G:S\$h=R5{            |

|> |URL: http://www.cs.hku.hk/~kkto    |  :t-6!Ghttp://www.cs.hku.hk/~kkto  |
|> \-----------------------------------^------------------------------------/

--
Ed Hook                              |       Coppula eam, se non posit
Computer Sciences Corporation        |         acceptera jocularum.
NASA Langley Research Center         | Me? Speak for my employer?...<*snort*>

Sun, 05 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

>>> The correct way to compare two floating point numbers for "equality" to
>>> use fabs(a-b) < eps..
>>  Or to use _normalized_ differences.

>> First version (has a division and a square root): fabs( a - b )/sqrt(
>> a*a + b*b ) < eps

Thomas> This only works correctly if a*a + b*b does not overflow your
Thomas> machine number, and also that a*a+b*b does not underflow.  Still
Thomas> some thought required...

It also falls over if 'a' is less then 'eps', but b is exactly zero; you
end up with

a/a < eps ,   ==>    1.0 < eps.

This is false, even though 'a' was within 'eps' of 'b'.

I ended up with a big clunky macro for robustly dealing with all this,
which uses two little macros for determining if a number is zero,
(i.e. absolute error, (a<eps), and one for relative_error (as explained
above).

#define real_equal(a,b) \
( ((a)==0.0) ?                    /*   if (a==0.0) {                     */ \
( ((b)==0.0) ?                /*       if (b==0.0) {                 */ \
TRUE                      /*           return TRUE;              */ \
:                           /*       } else {                      */ \
is_zero (b) )             /*           return is_zero(b);        */ \
:                                 /*   } else {                          */ \
( ((b)==0.0) ?                /*       if (b==0.0) {                 */ \
is_zero(a)                /*           return is_zero(a);        */ \
:                            /*       } else {                      */ \
eq_rel_error(a,b,EPSILON) /*           return eq_rel_error;      */ \
)                             /*       }                             */ \
)                                 /*   }                                 */

Does anybody have a better solution ? It seems a bit horrific !

--
=~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~=

Computer Science Department, The University, Bristol, UK

Sun, 05 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

>> The correct way to compare two floating point numbers for "equality" to
>> use fabs(a-b) < eps.  The problem arises in determining what to use for

>#define TOLERENCE 1.0e6
>#define EQUAL( a, b ) (fabs(a-b)<=(fabs(a)/TOLERENCE))

I use a slightly more "expensive", and slightly less standard way.  The
advantage is that it works well for pairs of numbers no matter how large or
small they are.  I don't have the code I use here at home, but it's something
like this.  One difference between the real function I use and this one is
that the real function has code to handle comparisons between a real zero and
a very small positive or negative double, and where the exponents are off by
one because one number is slightly less than a power of two and the other is
slightly more.

The function "frexp()" is a POSIX fuction that breaks a floating
point number into a fraction and exponent.

Boolean dequal(double x, double y)
{
int      x_exp, y_exp;
double   x_mant, y_mant;

x_mant = frexp(x, &x_exp);
y_mant = frexp(y, &y_exp);

if (x_exp != y_exp)
return False;

if (fabs(x_mant - y_mant) < 1.0e6)
return True;

return False;

Quote:
}

--

"While you're waiting, read the free novel we sent you.  It's a Spanish story
about a guy named `Manual'" - Dilbert

Sun, 05 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

>[ Followup-To: comp.lang.c.moderated ]

>>: Having a hard time to figure out a function on how to compare doubles.
>>: Any ideas?

>>No function necessary.. C has this really cool boolean equality operator
>>that does this sort of thing for you.  Have a look:

>It is unwise to compare floating point numbers for equality without
>thinking hard first.

...as evidenced by the fact that many of the proposed "solutions" have
potential overflow or underflow problems.

Of course, it's often necessary to first make sure you've defined the
algorithm for equality that you want, then make sure it works for all
possible input.  The algorithm you want is highly dependent on the
application (whose context was either never posted or I missed it).

For example, if you're doing a homework assignment and are waiting for a
numerical calculation to stabilize, and you know more or less the range
of input values for the comparison, then define an epsilon, use a simple
test with fabs, and be done with it.

If you're writing a database server though, you're likely to need a
tunable comparison (one user wants 10 digits of precision, another wants
14), that's well documented, and works correctly for all input values.

I can think of two broad categories of equality algorithms.  I'll call
them interval and difference, for want of knowing if they have official
names.

An interval method would partition the real numbers into contiguous
subsets or intervals within which all real numbers would be defined as
"equal" for the purpose of the aplication.  For example, the interval
called "1.65" might include all real numbers > 1.645 and <= 1.655.

This has an advantage in doing some kinds of applications because it
corresponds to how we often think of currency.  In the US, we round
monetary figures to the nearest .01 dollars, for example.  The
disadvantage of this method is that at all places in the application
where conversion is done between these intervals and real numbers,
there's a potential for unexpected behavior.  1.646 and 1.653 are
considered "equal" but 1.646 and 1.644 are unequal.  More specifically,
a calculation that we expect to result in 1.645 could be labelled 1.64
or 1.65 depending on the order of the operands, etc.

This method is often most appropriate for systems that convert input
values once into a decimal format and performa calculations with them in
that format.

A difference method relies on the relative difference between two
numbers.  Instead of translating between real numbers and these
intervals, the results are kept in the natural machine format.  (Note
however that the natural machine format already is an interval of the
real numbers, but hopefully at a finer grain than what the application
cares about.)  When two of these numbers are compared, they are equal if
their difference is less than some other number.  This other number
could be a fixed amount (like .01) or an amount relative to the numbers
being compared (say, .01 times the number with the lowest magnitude).

One advantage of this method is that there are no discontinuities like
there is with the interval method.  A disadvantage is that if A=B and
B=C, A and C are not necessarily equal.

One can even combine fixed and relative comparisons of the difference
method for difference equality.  It might be desireable to have a
relative difference of .00001 times the number with the lowest magnitude
for equality (essentially comparing to 5 digits), but for values close to
0 use a fixed difference like .00001.  This ensures that results that
you expect to be 0 compare as equal, even if one is twice the other.
(For example, one result might be the lowest machine representable
value, and anotehr result might be twice that value, yet you still
probably want those to compare equal in many applications.)

Once you've decided on the algorithm that's most appropriate for your
application, you still have to implement it so that there's no chance of
over- or underflow.  This can be non-trivial, especially if you also
want to do it quickly in terms of run-time performance.

Also remember that if you have a special function for equality of real
numbers, you need to be careful if you're using comparison operators:

if (a < b)

might instead need to be written as:

if ((!real_equal(a,b)) && (a<b))

to prevent unexpected results!

As Thomas Koenig said, one should think carefully about what you need
before answering this question.  There is no single correct method of
doing the comparison.

Michael J. Zehr
Sr. Software Engineer
Kenan Systems Corporation

Sun, 05 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

>>: Having a hard time to figure out a function on how to compare doubles.
>>: Any ideas?
>The correct way to compare two floating point numbers for "equality" to
>use fabs(a-b) < eps.  The problem arises in determining what to use for
>eps, of course :-).

Back when I still did numerical programming (in fortran) I used
the equilvalent of

eps1 >= (fabs(a-b) / fmax(fabs(a), fabs(b), eps2)

If you are converging on a large number, then you want a relative
equality test (agree to X decimal places), however, if the sequence
happens to converge on zero than a relative test won't work, so
below "eps2" you switch over to an absolute "fuzz factor"

Sun, 05 Jul 1998 03:00:00 GMT  Help, best way to compare doubles
Quote:

>>: Having a hard time to figure out a function on how to compare doubles.
>>: Any ideas?

..>The correct way to compare two floating point numbers for "equality" to

Quote:
>use fabs(a-b) < eps..

You'll be in trouble if you apply the transitive law to this:
a approximately equal to b, b approximately equal to c
will not imply
a approximately equal to c.

Any epsilon you choose in the absence of further information
is likely to be arbitrary.  If it's important that two floats are or are
not equal, there's usually some other logical evidence (which you
may have foolishly thrown away) elsewhere in the program telling
you so.

If there isn't any evidence, then more often than not it's not terribly
important to tell the difference between the case of the floats being
equal, and the case of the floats being closer together than epsilon,
so use (fabs(a-b) < eps)

You should only mean what you say and say what you mean.
If you can't rely on the statement (a == b) then don't say something
else like (fabs(a-b) < eps) instead, because you won't mean it,
and you should only say what you mean, you see?

Julian Todd.

Sun, 05 Jul 1998 03:00:00 GMT  Help, best way to compare doubles

Quote:

> >>: Having a hard time to figure out a function on how to compare doubles.
> >>: Any ideas?
> ...
> >The correct way to compare two floating point numbers for "equality" to
> >use fabs(a-b) < eps..

> Or to use _normalized_ differences.
Better

> First version (has a division and a square root):
>     fabs( a - b )/sqrt( a*a + b*b ) < eps

A usually faster variant (assuming absolute values are taken faster then
multiplications:
fabs( a - b )/( fabs(a) + fabs(b) ) < eps
Avoiding the square also significantly reduces the chances of overflow or
underflow during evaluation of the expressions. Try to compare the
doubles 1e200 and 2e200 or 1e-200 and 2e-200 with both routines on a machine
that uses IEEE double precision arithmetics. Squaring them will give overflow
or underflow.
In the file <float.h> some usefull constants are defined: DBL_MAX, DBL_MIN
and especially DBL_EPSILON, the minimum positive floating point number x
such that 1.0 + x != 1.0  In most cases you can chose eps a few times
DBL_EPSILON when using normalised differences.
Quote:

> better (has no division):
>     fabs( a - b ) < eps * sqrt( a*a + b*b )

for completeness:
fabs( a - b ) < eps * ( fabs(a) + fabs(b) )
Quote:

> better yet (neither division nor square root):
>    ( t = a - b )*t < eps * ( a*a + b*b )

^
YUCK the order of evaluation is UNDEFINED does this t mean the NEW or the
ORIGINAL value of t?
(squaring argument remains also valid over here)

--
===============================================================================

There's no computer like an ARMed computer!
===============================================================================
... Sorry... my mind has a few bad sectors.

Sun, 05 Jul 1998 03:00:00 GMT  Help, best way to compare doubles
Slightly edited,

Quote:

>>>The correct way to compare two floating point numbers for "equality" to
>>>use fabs(a-b) < eps..
>>Or to use _normalized_ differences.
>>    fabs( a - b )/sqrt( a*a + b*b ) < eps
>>    fabs( a - b ) < eps * sqrt( a*a + b*b )
>>   ( t = a - b )*t < eps * ( a*a + b*b )

^^^^^^^^^^^^^^^ Others have pointed out that the assignment and
reference to t results in undefined behavior, but I will anyway,
in the hope that I can con the moderator into passing the rest
of the article.

Quote:
>     Nice.  Is there a general method you are using here?  Details
>would be nice.  If you couch it in terms of C, the moderator may
>whistle and look the other way <G>.

To _normalize_ a vector (or a quantity computed from the components of
a vector, as in this case) is to divide its components by the norm (or
length) of the vector.  If we treat the values (a,b) as a vector, its
length by the usual "Euclidian" definition is sqrt(a^2 + b^2) (where ^
is used here to mean "raised to the power", a common notation in groups
that are not comp.lang.c.*).

The point is this: If you are interested in computing results to (say)
six digits of accuracy, then you would want a = 100000 and b =
100000.5 to compare equal, but if you compare using fabs(a - b) <
1.0e-6, they will not, whereas a = 0.1 and b = 0.1000005 will.  The
process of dividing by the norm results in the comparisons being made
in a way that takes the number of digits of agreement into account
while factoring out the influence of the magnitudes of the values.

Some comments: (1) The Euclidian norm is not the only one that can be
used.  You can get the same effect (though with slightly different
results for some comparisons) by using max(fabs(a),fabs(b)) or
(fabs(a)+fabs(b)) as the norm.  This saves a sqrt() call (at the cost
of fabs() calls), and the "max norm" does not risk overflow.

(2) You may not want to normalize comparisons of very small numbers.
For example, if six digits is your goal, then you may want to accept
0.000001 and 0.0000015 as equal, but a normalized comparison will
not do so.  This situation is often dealt with by using
max(1.0, fabs(a), fabs(b)) or (1.0 + fabs(a) + fabs(b)) as the "norm".
I'm not sure if anything other than max(1.0, sqrt(a^2+b^2)) is needed
to deal with the Euclidian version.
--
Matthew Saltzman
Clemson University Math Sciences

Sun, 05 Jul 1998 03:00:00 GMT

 Page 1 of 3 [ 31 post ] Go to page:   

Relevant Pages