Comparing doubles

[Thomas Kowalski]

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

 

[ralpe]

[...]

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

};

Why is this a union?

[...]

 

[Michael DOUBEZ]

Thomas Kowalski a écrit :

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?

I did try to make an switchsign() fuction by reseting the SIGN bit. It
was less efficient than a=-a;

And you realize of course that you may succeed on your compiler and
architecture but it will likely fail on others'.

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;
};

Here, you must mean struct.

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));
}

Should be >>. Little endian inverts bytes, not bits

Michael

 

[osmium]

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?

<snip>

I think bit fiddling should be a last resort, it forces anyone encountering
the code to go through the same torturous process you went through when you
coded it.

Have you looked to see what you can do with the notion of "relative error"?

 

[Thomas Kowalski]

Why is this a union?>
Right, this should of course be a struct.

Regards,
Tk

 

[Thomas Kowalski]

Have you looked to see what you can do with the notion of "relative error"?

Something like the following don't work really well, either for small
numbers (e.g. 0.0000031234434).

bool isEqual ( double a, double b ) {
return ( fabs ((a-b)/b < std::numeric_limits<double>::epsilon() );
}

 

[Thomas Kowalski]

But this a approach usually fails if comparing doubles of different

std::numeric_limits<double>::epsilon() is a good starting point.

Does it work for this test set? Maybe I did something wrong, but for
me its not working well.
ASSERT ( isEqual (1234567890123456, 1234567890123457 ) );
ASSERT ( isEqual (1.234567890123456, 1.234567890123457 ) );
ASSERT ( isEqual (0.000001234567890123456,
0.000001234567890123457 ) );

I did try to make an switchsign() fuction by reseting the SIGN bit. It
was less efficient than a=-a;

And you realize of course that you may succeed on your compiler and
architecture but it will likely fail on others'.

I do. It's a pity that there is no better way to do this (or at least
I don't know it).

Here, you must mean struct.

Yes, true.

Should be >>. Little endian inverts bytes, not bits

With ">>" its not really working either. I assume there is a problem
with shifting a 52 bit number?

Regards,
Thomas

 

[Michael DOUBEZ]

Thomas Kowalski a écrit :

Something like the following don't work really well, either for small
numbers (e.g. 0.0000031234434).

bool isEqual ( double a, double b ) {
return ( fabs ((a-b)/b < std::numeric_limits<double>::epsilon() );
}

If a = b + 10 and a!=b and b very big, the above assumption will be
true but a!=b.

The equality can be:
abs(x - y) < std::numeric_limits<double>::epsilon()

And inequality less or equal can be of the same format without abs:
x - y < std::numeric_limits<double>::epsilon();

Michael

 

[Michael DOUBEZ]

Thomas Kowalski a écrit :

Does it work for this test set? Maybe I did something wrong, but for
me its not working well.
ASSERT ( isEqual (1234567890123456, 1234567890123457 ) );
ASSERT ( isEqual (1.234567890123456, 1.234567890123457 ) );
ASSERT ( isEqual (0.000001234567890123456,
0.000001234567890123457 ) );

On my system, epsilon is 2.22045e-16. It is normal the first two assert
trigger a fault, the third not.

I do. It's a pity that there is no better way to do this (or at least
I don't know it).


Yes, true.


With ">>" its not really working either. I assume there is a problem
with shifting a 52 bit number?

No reason for that. Perhaps you can try simply to mask the double:
union dl
{
double d;
unsigned long l;
};

dl d1=0.000001234567890123456;
dl d2=0.000001234567890123457;

unsigned long mask = (~0)<<nb_bit;

if ( (d1.l&mask) == (d2.l&mask) )
{
// Equal !!!
}

I didn't try but it is not far from the solution you seek.

Michael

 

[Victor Bazarov]

Thomas Kowalski a écrit :

On my system, epsilon is 2.22045e-16. It is normal the first two
assert trigger a fault, the third not.


No reason for that. Perhaps you can try simply to mask the double:
union dl
{
double d;
unsigned long l;
};

dl d1=0.000001234567890123456;
dl d2=0.000001234567890123457;

unsigned long mask = (~0)<<nb_bit;

if ( (d1.l&mask) == (d2.l&mask) )
{
// Equal !!!
}

I didn't try but it is not far from the solution you seek.

The problem with using 'union' like that is that it's undefined.

If portability is a requirement, it's better to extract the mantissa
from each value and compare those; see 'frexp' function in <cmath>.

V

 

[Victor Bazarov]

The problem with using 'union' like that is that it's undefined.

If portability is a requirement, it's better to extract the mantissa
from each value and compare those; see 'frexp' function in <cmath>.

Of course, I've forgotten to mention that exponents need to be the same,
IOW the mantissas obtained from 'frexp' need to be corrected for the
minimal exponent before comparison...

V

 

[Victor Bazarov]

Of course, I've forgotten to mention that exponents need to be the
same, IOW the mantissas obtained from 'frexp' need to be corrected
for the minimal exponent before comparison...

Damn... I haven't paid attention before posting, have I? There is
no need to correct the mantissas. Just compare the exponents and if
they are not the same, you have your answer, if they are the same,
compare the mantissas using the chosen epsilon.

Sorry about the fuss...

V

 

[Juha Nieminen]

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.

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.
I suppose that the fabs(a-b)<epsilon comparison is a compromise
between good-enough accuracy and speed, as all those operations are
quite fast to perform.

If what you want is a comparison which does not depend on the
magnitude of the values and you don't care if the comparison becomes
somewhat slower, then what you want is to compare the first n
significant digits of the values.
This can be done by dividing the values by their magnitude.

The magnitude of a value (ie. how many digits it has at the left
of the decimal point) can be calculated with floor(log10(fabs(value))).
Now if you divide the value by pow(10, magnitude) you get a value
which is scaled to the range 0-1.
In other words: svalue = value/pow(10, floor(log10(fabs(value))));
Now svalue is value scaled to 0-1, which makes it easy to compare the
n significant digits of the value. You can divide the other value
with that same divisor and then do the fabs(svalue1-svalue2)<epsilon
comparison.

 

[Victor Bazarov]

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.

That statement is too vague. And there is no way to be more concrete
without any additional information: it all depends on the problem
domain. What do those numbers designate? Example: such comparison
is perfectly fine if 'a' and 'b' are coordinates and the 'isEqual'
means that the gap between those coordinates is smaller than some
predefined value ('THRESHOLD').

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.

Not true. It's no more restrictive within the imposed requirements.
If the requirements are such that the absolute difference between
the two values is smaller than some value, it's coded correctly.
The problem is that we just don't know what the requirements are.

When a and b are very small then epsilon
effectively becomes larger with respect to them and the inaccuracy
of the comparison increases.
I suppose that the fabs(a-b)<epsilon comparison is a compromise
between good-enough accuracy and speed, as all those operations are
quite fast to perform.

As I showed above, it actually may be the only solution acceptable.

If what you want is a comparison which does not depend on the
magnitude of the values and you don't care if the comparison becomes
somewhat slower, then what you want is to compare the first n
significant digits of the values.
This can be done by dividing the values by their magnitude.

There is no single magnitude when *two* numbers are concerned. How
do you decide what that magnitude is?

The magnitude of a value (ie. how many digits it has at the left
of the decimal point) can be calculated with
floor(log10(fabs(value))).

Right. Of *a* value.

Now if you divide the value by pow(10,
magnitude) you get a value which is scaled to the range 0-1.

Actually, it's scaled to the range 0.1-1.

In other words: svalue = value/pow(10, floor(log10(fabs(value))));
Now svalue is value scaled to 0-1,
0.1-1

which makes it easy to compare the
n significant digits of the value.

Yes. Compare to what?

You can divide the other value
with that same divisor and then do the fabs(svalue1-svalue2)<epsilon
comparison.

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

What you're describing here is comparison of significant digits. To
do it quickly and without extra calculation you can just do

if (fabs(a-b) < ADJUSTED_TRHESHOLD * max(fabs(a),fabs(b)) ) ...

which basically is the calculation of the relative error (kind of).
But nothing actually says that in the OP's problem domain it is the
right thing to do.

V

 

[kwikius]

Thomas Kowalski a écrit :




std::numeric_limits<double>::epsilon() is a good starting point.

IMO its useful to be able to supply the epsilon . Being given a hard
coded one seems arbitrary. You can also get more elaborate:

#include <limits>

//compare against some acceptable error range
// epsilon can be computed
// and may therefore be negative
// hence use of abs
inline
bool eq_( double lhs, double rhs, double epsilon)
{
return std::abs(lhs-rhs) <= std::abs(epsilon);
}

//overload for function to supply epsilon
inline
bool eq_( double lhs, double rhs, double (*f)() )
{
return std::abs(lhs-rhs) <= std::abs(f());
}

//overload for function to supply epsilon based on inputs
inline
bool eq_( double lhs, double rhs, double (*f)(double, double) )
{
return std::abs(lhs-rhs) <= std::abs(f(lhs,rhs));
}

// get some percent based on x,y
template <int N>
double exp10range( double x, double y)
{
double avg =(std::abs(x) + std::abs(y))/2;
return avg * std:ow(10.,N);
}

#include <iostream>
int main()
{
double x = -7.;
double y = -7.007;

std::cout << eq_(x,y, exp10range<-1> ) << '\n';
std::cout << eq_(x,y, exp10range<-2> ) << '\n';
std::cout << eq_(x,y, exp10range<-3> ) << '\n';
std::cout << eq_(x,y, exp10range<-4> ) << '\n';
}

regards
Andy Little

 

[Michael DOUBEZ]

kwikius a écrit :

Thomas Kowalski a écrit :

std::numeric_limits<double>::epsilon() is a good starting point.

IMO its useful to be able to supply the epsilon . Being given a hard
coded one seems arbitrary. You can also get more elaborate:

[snip]

I see what you mean but in the case of a customizable threshold, I
wouldn't call it an equality but a proximity.

All boils down to "what do the OP mean by isEqual(double,double) ?" as
Victor Bazarov pointed out.


Michael

 

[Thomas Kowalski]

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.

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).

If what you want is a comparison which does not depend on the
magnitude of the values

Yes, that's what I need.

and you don't care if the comparison becomes
somewhat slower, then what you want is to compare the first n
significant digits of the values.

That's the problem. In my case speed is second just to correctness.
Readability or portability on the other hand have no importance at
all. That's why I wanted to resort to the bit fiddling.

Regards,
Thomas Kowalski

 

[Thomas Kowalski]

It is indeed true that the classic comparison fabs(a-b)<epsilon

That statement is too vague.

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).

It's no more restrictive within the imposed requirements.
If the requirements are such that the absolute difference between
the two values is smaller than some value, it's coded correctly.

Well, in my scenario it is incorrect.

There is no single magnitude when *two* numbers are concerned. How
do you decide what that magnitude is?

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


Regards,
Thomas Kowalski

 

[Thomas Kowalski]

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.

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

or in case that a might be zero:

bool isEqual(double a, double b) {
if (a = 0.0) {
return fabs(b) < THRESHOLD;
}
else {
return fabs(1.0-(b/a)) < THRESHOLD;
}
}


The question on hand is now. Does anyone have an idea whether and how
this (or something similiar) could be speeded up (using C, not
resorting to assembler) ?

Regards,
Thomas Kowalski

 

[kwikius]

All boils down to "what do the OP mean by isEqual(double,double) ?" as
Victor Bazarov pointed out.

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

regards
Andy Little