Velocity Reviews > C++ > Comparing doubles

Comparing doubles

Thomas Kowalski
Guest
Posts: n/a

 07-09-2007
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

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

Thomas Kowalski

ralpe
Guest
Posts: n/a

 07-09-2007
On Jul 9, 12:07 pm, Thomas Kowalski <(E-Mail Removed)> wrote:

> [...]
>
> union IEEErepresentation {
> long long man : 52;
> int exp : 11;
> int sign : 1;
>
> };

Why is this a union?

> [...]

Michael DOUBEZ
Guest
Posts: n/a

 07-09-2007
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.

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

> 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

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
Guest
Posts: n/a

 07-09-2007
"Thomas Kowalski" wrote:

> 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

<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
Guest
Posts: n/a

 07-09-2007
> Why is this a union?>
Right, this should of course be a struct.

Regards,
Tk

Thomas Kowalski
Guest
Posts: n/a

 07-09-2007
On Jul 9, 2:43 pm, "osmium" <(E-Mail Removed)> wrote:

> 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
Guest
Posts: n/a

 07-09-2007
> > 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.

>
> 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
Guest
Posts: n/a

 07-09-2007
Thomas Kowalski a écrit :
> On Jul 9, 2:43 pm, "osmium" <(E-Mail Removed)> wrote:
>
>> 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() );
> }
>

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
Guest
Posts: n/a

 07-09-2007
Thomas Kowalski a écrit :
>>> 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.

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

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

>
>
>> 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?

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;

{
// Equal !!!
}

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

Michael

Victor Bazarov
Guest
Posts: n/a

 07-09-2007
Michael DOUBEZ wrote:
> Thomas Kowalski a écrit :
>>>> 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.
>>> 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 ) );

>
> On my system, epsilon is 2.22045e-16. It is normal the first two
> assert trigger a fault, the third not.
>
>>
>>
>>> 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?

>
> 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;
>
> {
> // 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
--