Comparing doubles

[jacob navia]

Hi everyone,
To determine equality of two doubles a and b the following is often
done:

bool isEqual ( double a, double b ) {
return ( fabs (a-b) < THRESHOLD );
}

But this a approach usually fails if comparing doubles of different
magnitude since it's hard or not possible to find a suitable threshold
value. Since an double consists of mantissa, exponent and sign I
assume that the "best" way in this case would be to just "ignore" the
last few (e.g. 4-8 bits) of the mantissa to check for equality (at
least for the x86 architecture, which follow IEEE 754). Have someone
here ever tried to implement a similar approach? If yes, which
experiences have been made?

However my first approach to implement this failed. Can anyone tell
what I might have done wrong?

union IEEErepresentation {
long long man : 52;
int exp : 11;
int sign : 1;
};

union IEEEdouble {
double d;
IEEErepresentation c;
};

const int number_of_bits_to_ignore = 8;

bool isEqual (double a, double b) {
IEEEdouble aa, bb;
aa.d = a;
bb.d = b;
return ((aa.c.man << number_of_bits_to_ignore) == (bb.c.man <<
number_of_bits_to_ignore));
}

int main() {
string a = isEqual (1324.5678901234, 1324.5678901256 ) ? "true" :
"false";
string b = isEqual (134.5678901234, 134.5678901256 ) ? "true" :
"false";
string c = isEqual (.5678901234, .5678901256 ) ? "true" : "false";
string d = isEqual (1324.5678901234, 124.5678901256 ) ? "true" :
"false";
string e = isEqual (134.5678901234, 134.5178901256 ) ? "true" :
"false";
string f = isEqual (.5678901234, .3678901256 ) ? "true" : "false";

cout << "a = " << a << endl;
cout << "b = " << b << endl;
cout << "c = " << c << endl;
cout << "d = " << d << endl;
cout << "e = " << e << endl;
cout << "f = " << f << endl;
return 0;
}

Thanks in advance,
Thomas Kowalski

The lcc-win32 compiler system provides fcmp:
--------------------------------------------------------cut here
/*
fcmp
Copyright (c) 1998-2000 Theodore C. Belding
University of Michigan Center for the Study of Complex Systems
<mailto:[email protected]>
<http://www-personal.umich.edu/~streak/>

This file is part of the fcmp distribution. fcmp is free software;
you can redistribute and modify it under the terms of the GNU Library
General Public License (LGPL), version 2 or later. This software
comes with absolutely no warranty. See the file COPYING for details
and terms of copying.

File: fcmp.c

Description: see fcmp.h and README files.
Knuth's floating point comparison operators, from:

Knuth, D. E. (1998). The Art of Computer Programming.
Volume 2: Seminumerical Algorithms. 3rd ed. Addison-Wesley.
Section 4.2.2, p. 233. ISBN 0-201-89684-2.

Input parameters:
x1, x2: numbers to be compared
epsilon: determines tolerance

epsilon should be carefully chosen based on the machine's precision,
the observed magnitude of error, the desired precision, and the
magnitude of the numbers to be compared. See the fcmp README file for
more information.

This routine may be used for both single-precision (float) and
double-precision (double) floating-point numbers.
*/

#ifdef HAVE_CONFIG_H
#include <config.h>
#endif

#include "fcmp.h"
#include <math.h>

int EXPORT fcmp(double x1, double x2, double epsilon) {
int exponent;
double delta;
double difference;

/* Get exponent(max(fabs(x1), fabs(x2))) and store it in exponent. */

/* If neither x1 nor x2 is 0, */
/* this is equivalent to max(exponent(x1), exponent(x2)). */

/* If either x1 or x2 is 0, its exponent returned by frexp would be 0, */
/* which is much larger than the exponents of numbers close to 0 in */
/* magnitude. But the exponent of 0 should be less than any number */
/* whose magnitude is greater than 0. */

/* So we only want to set exponent to 0 if both x1 and */
/* x2 are 0. Hence, the following works for all x1 and x2. */

frexp(fabs(x1) > fabs(x2) ? x1 : x2, &exponent);

/* Do the comparison. */

/* delta = epsilon * pow(2, exponent) */

/* Form a neighborhood around x2 of size delta in either direction. */
/* If x1 is within this delta neighborhood of x2, x1 == x2. */
/* Otherwise x1 > x2 or x1 < x2, depending on which side of */
/* the neighborhood x1 is on. */

delta = ldexp(epsilon, exponent);

difference = x1 - x2;

if (difference > delta)
return 1; /* x1 > x2 */
else if (difference < -delta)
return -1; /* x1 < x2 */
else /* -delta <= difference <= delta */
return 0; /* x1 == x2 */
}

#ifdef TEST
#include <float.h>
#include <stdio.h>
#define ASSERT assert
#include <assert.h>
int main() {
int result;
ASSERT ( fcmp (1234567890123456, 1234567890123457 ,DBL_EPSILON) );
ASSERT ( fcmp (1.234567890123456, 1.234567890123457 ,DBL_EPSILON) );
ASSERT ( fcmp (0.000001234567890123456,
0.000001234567890123457,DBL_EPSILON ) );
return 0;
}
#endif
---------------------------------------------------------end of fcmp.c

 

[Lionel B]

Suggested function names:

in_the_right_field(x,y); //courtesy of an old Boatbuilder. It was his
favourite saying.

in_the_ballpark(x,y);

so_far_off_the_line_you_cant_even_see_the_line(x,y) // Joey in "Friends"

on_another_planet(x,y); // usually used by others to describe myself ;-)

void nail_him_Tony(x,y);
// Favourite saying of a roofer I knew when a batten (or whatever) was
in the right spot, according to him anyway;-). Actually was "Nail 'im
Tone' ! ". Tony was his long suffering Boss

I have a series of functions nad_eq, nad_gt, nad_lt, ... as in "near-as-
dammit". Also has (possibly quite appropriate) echoes of "nads" as in "a
kick in the 'nads" (sorry for that).

 

[kwikius]

Hello everyone,
thanks to the input I got in this thread I figured out some possible
solution for the problem might be something like this:

const double THRESHOLD =
4096.0*std::numeric_limits<double>::epsilon(); // or any other
multiple of epsilon.

could do :

template <unsigned int N>
double threshold()
{
enum { mux = (1 << N) };
return mux * std::numeric_limits<double>::epsilon();
}

regards
Andy Little

 

[James Kanze]

To be more a bit more precise. I have two numbers of similar magnitude
(nearly the same exponent, different mantissa). I perform some
calculations like e.g. additions and multiplications. If I want to
check whether the number is in some range (let's say 0.0 to 1.0 ... or
any other arbitrary numbers) it's clear that thanks to rounding errors
I might get the wrong result from time to time (with a few billion
numbers it's nearly always we case).

If you want to check whether the number is in some range, then
it's clear that you don't want to compare for equality, but for
inequality. And that all that has been said before is totally
irrelevant.

Therefore I need a way to just
ignore the last few digits of the mantissa to eliminate rounding
errors while I want to preserve the accuracy of the _isEqual_
operation ( the operation isEqual is expected to tell me whether the
numbers I am comparing would be the same in case the calculation
(additions etc.) would be precise).

Well, it's pretty easy to ignore the last few digits of the
mantissa:

union U
{
double d ;
uint64_t u ;
} ;
U tmp ;
tmp.d = value ;
u &= ~(uint64_t(0)) << 4 ;
value = tmp.d ;

(Technically, this is undefined behavior, and in addition, it
uses a type which won't be present until the next version of the
standard. But pratically, it will do what you want on just
about any machine I can think of.)

Of course, it's also pretty useless in most cases. All you're
doing is creating new equivalence classes, but two numbers that
were originally only one bit off might still compare not equal.

If you're interested in a "fuzzy" equality, there are a number
of classical solutions, along the lines of "abs( (a-b)/(a + b) )
< someEpsilon", but be aware that they do NOT define equivalence
classes, and do not have the same properties of equality.
(They're not transitive, for example.) They also typically have
some weaknesses, which may require special casing. (Try
comparing 0.0 and 0.0---which I think we agree should be
equal---with the above, for example.)

Of course you are right. That's why "==" is a binary operator. Let's
just assume that same magnitude means same exponent.

frexp() ?

But again, I fail to see how this could be of any use.

 

[Pete Becker]

To be more a bit more precise. I have two numbers of similar magnitude
(nearly the same exponent, different mantissa). I perform some
calculations like e.g. additions and multiplications. If I want to
check whether the number is in some range (let's say 0.0 to 1.0 ... or
any other arbitrary numbers) it's clear that thanks to rounding errors
I might get the wrong result from time to time (with a few billion
numbers it's nearly always we case). Therefore I need a way to just
ignore the last few digits of the mantissa to eliminate rounding
errors while I want to preserve the accuracy of the _isEqual_
operation ( the operation isEqual is expected to tell me whether the
numbers I am comparing would be the same in case the calculation
(additions etc.) would be precise).

Now all that's left is to figure out which of the "last few digits"
reflect rounding errors and which ones are actually part of the
"precise" result. If you don't have a rigorous procedure for doing
that, the results you get from jamming zeros (or ones, for that matter)
into the "last few digits" are no better than the original values.
Doing that requires analyzing your algorithm to see where roundoff
errors are occurring. Having done that, you no longer need this sort of
sledge hammer.

 

[Pete Becker]

That's exactly what I mean. In my case I can't tell the magnitude of a
number. (Lets just define the magnitude as the exponent of the double,
since clearly 2 doubles are not equal if the exponent is not equal).

But there's no such requirement in the C++ standard, and in general,
it's not true. 1.0x10^3 is equal to 0.1x10^4, despite their different
exponents. Granted, most floating-point representations these days use
normalized values, so representable values with different magnitudes
have different values. But that's hardware specific, not guaranteed.

 

[Juha Nieminen]

Why divide *the other* value with *the same divisor* and not the other
way around?

Why not? You are just comparing the equality/inequality of two
values. If they have completely different orders of magnitude
(ie. exponent of 10) then it doesn't matter which one is divided
with the magnitude of the other value. They will just be vastly
different and the fabs(svalue1-svalue2)<epsilon will give, correctly,
the result that they are different.
It's only when the two values are really close to each other that
the way you compare becomes significant. When they are very close to
each other they have the same order of magnitude (give or take), in
which case dividing them by the same value scales them both properly.

So to answer your question: It doesn't matter which one of the
values you use to calculate the divisor for both values.

 

[James Kanze]

That's exactly what I mean. In my case I can't tell the magnitude of a
number. (Lets just define the magnitude as the exponent of the double,
since clearly 2 doubles are not equal if the exponent is not equal).

Clearly they're not equal if the mantissas are not equal,
either. I thought you were looking for some "almost equal"
function. I'd say that 1.9999999999999998 and 2 are about
equal---there's only one LSB difference in them (in an IEEE
double, at least). But they have different exponents.

 

[kwikius]

It is indeed true that the classic comparison fabs(a-b)<epsilon
is not very good if a and b are very large or very small. When they
are very large the comparison becomes more restrictive because
epsilon becomes much smaller in relation to a and b, and thus less
discrepancy is allowed. When a and b are very small then epsilon
effectively becomes larger with respect to them and the inaccuracy
of the comparison increases.

BTW One of the strengths of my Quan library is that you can deal with
very large or small values accurately. Its all down to the SI system,
so you can work in (e.g) nM or Mm if you wish:

http://tinyurl.com/22xhne

The example shows use of non si values, which reduces accuracy but
using only SI units will give much better results. The accuracy
attained in the above example is not possible with doubles. I won't go
into the deatils here though...

http://sourceforge.net/projects/quan

(n.b quan is soon to be superceded by my 'super' quan library.)

regards
Andy Little