Comparing doubles

**Michael DOUBEZ** · Jul 9 '07, 11:15 AM

Re: Comparing doubles

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_li mits<double>::e psilon() 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
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 IEEErepresentat ion {
long long man : 52;
int exp : 11;
int sign : 1;
};

Here, you must mean struct.

>
union IEEEdouble {
double d;
IEEErepresentat ion 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

**ralpe** · Jul 9 '07, 11:15 AM

Re: Comparing doubles

On Jul 9, 12:07 pm, Thomas Kowalski <t...@gmx.dewro te:

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

Why is this a union?

[...]

**osmium** · Jul 9 '07, 12:45 PM

Re: Comparing doubles

"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
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** · Jul 9 '07, 01:15 PM

Re: Comparing doubles

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

Regards,
Tk

**Thomas Kowalski** · Jul 9 '07, 02:05 PM

Re: Comparing doubles

On Jul 9, 2:43 pm, "osmium" <r124c4u...@com cast.netwrote:

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_li mits<double>::e psilon() );
}

**Thomas Kowalski** · Jul 9 '07, 02:15 PM

Re: Comparing doubles

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_li mits<double>::e psilon() 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 (12345678901234 56, 123456789012345 7 ) );
ASSERT ( isEqual (1.234567890123 456, 1.2345678901234 57 ) );
ASSERT ( isEqual (0.000001234567 890123456,
0.0000012345678 90123457 ) );

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** · Jul 9 '07, 02:15 PM

Re: Comparing doubles

Thomas Kowalski a écrit :

On Jul 9, 2:43 pm, "osmium" <r124c4u...@com cast.netwrote:
>

>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_li mits<double>::e psilon() );
}
>

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_li mits<double>::e psilon()

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

Michael

**Michael DOUBEZ** · Jul 9 '07, 02:25 PM

Re: Comparing doubles

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_l imits<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 (12345678901234 56, 123456789012345 7 ) );
ASSERT ( isEqual (1.234567890123 456, 1.2345678901234 57 ) );
ASSERT ( isEqual (0.000001234567 890123456,
0.0000012345678 90123457 ) );

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.0000012345 67890123456;
dl d2=0.0000012345 67890123457;

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** · Jul 9 '07, 03:05 PM

Re: Comparing doubles

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 (12345678901234 56, 123456789012345 7 ) );
>ASSERT ( isEqual (1.234567890123 456, 1.2345678901234 57 ) );
>ASSERT ( isEqual (0.000001234567 890123456,
>0.000001234567 890123457 ) );

>
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
>>architectur e 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.0000012345 67890123456;
dl d2=0.0000012345 67890123457;
>
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
--
Please remove capital 'A's when replying by e-mail
I do not respond to top-posted replies, please don't ask

**Victor Bazarov** · Jul 9 '07, 03:05 PM

Re: Comparing doubles

Victor Bazarov wrote:

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 (12345678901234 56, 123456789012345 7 ) );
>>ASSERT ( isEqual (1.234567890123 456, 1.2345678901234 57 ) );
>>ASSERT ( isEqual (0.000001234567 890123456,
>>0.00000123456 7890123457 ) );

>>
>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
>>>architectu re 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.0000012345 67890123456;
>dl d2=0.0000012345 67890123457;
>>
>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>.

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
--
Please remove capital 'A's when replying by e-mail
I do not respond to top-posted replies, please don't ask

**Victor Bazarov** · Jul 9 '07, 03:05 PM

Re: Comparing doubles

Victor Bazarov wrote:

Victor Bazarov wrote:

>Michael DOUBEZ wrote:

>>Thomas Kowalski a écrit :
>>>>>But this a approach usually fails if comparing doubles of
>>>>>differen t magnitude since it's hard or not possible to find a
>>>>>suitable threshold value.
>>>>std::numeri c_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 (12345678901234 56, 123456789012345 7 ) );
>>>ASSERT ( isEqual (1.234567890123 456, 1.2345678901234 57 ) );
>>>ASSERT ( isEqual (0.000001234567 890123456,
>>>0.0000012345 67890123457 ) );
>>>
>>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.0000012345 67890123456;
>>dl d2=0.0000012345 67890123457;
>>>
>>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>.

>
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
--
Please remove capital 'A's when replying by e-mail
I do not respond to top-posted replies, please don't ask

**Juha Nieminen** · Jul 9 '07, 06:45 PM

Re: Comparing doubles

Thomas Kowalski wrote:

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(fab s(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(fab s(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)<epsilo n
comparison.

**Victor Bazarov** · Jul 9 '07, 07:15 PM

Re: Comparing doubles

Juha Nieminen wrote:

Thomas Kowalski wrote:

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

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(fab s(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(fab s(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)<epsilo n
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_TRHESH OLD * max(fabs(a),fab s(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
--
Please remove capital 'A's when replying by e-mail
I do not respond to top-posted replies, please don't ask

**kwikius** · Jul 9 '07, 10:15 PM

Re: Comparing doubles

On 9 Jul, 12:12, Michael DOUBEZ <michael.dou... @free.frwrote:

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_li mits<double>::e psilon() 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(epsilo n);
}

//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::pow(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

Comparing doubles

Comparing doubles

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment

Comment