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

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

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

[ralpe]

On Jul 9, 12:07 pm, Thomas Kowalski <[e-mail]:

[...]

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

};

Why is this a union?

[...]

[osmium]

"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]

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

Regards,
Tk

[Thomas Kowalski]

On Jul 9, 2:43 pm, "osmium" <[e-mail]:

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

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]

Thomas Kowalski a écrit :

On Jul 9, 2:43 pm, "osmium" <[e-mail]:

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

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;

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]

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;

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]

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

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]

Victor Bazarov wrote:

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

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]

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(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]

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

[kwikius]

On 9 Jul, 12:12, Michael DOUBEZ <[e-mail]:

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.

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

[Michael DOUBEZ]

kwikius a écrit :

On 9 Jul, 12:12, Michael DOUBEZ <[e-mail]:

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

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

is not very good if a and b are very large or very small.

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]

On 10 Jul, 07:48, Michael DOUBEZ <[e-mail]:

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

[jacob navia]

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?

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:Te*********@umich.edu>
<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]

On Tue, 10 Jul 2007 01:47:00 -0700, kwikius wrote:

On 10 Jul, 07:48, Michael DOUBEZ <[e-mail]:

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

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

--
Lionel B

[kwikius]

On 10 Jul, 09:27, Thomas Kowalski <[e-mail]:

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]

On Jul 10, 9:36 am, Thomas Kowalski <[e-mail]:

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.

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

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.

frexp() ?

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

--
James Kanze (GABI Software) email:ja*********@gmail.com
Conseils en informatique orientée objet/
Beratung in objektorientierter Datenverarbeitung
9 place Sémard, 78210 St.-Cyr-l'École, France, [phone]

[Pete Becker]

On 2007-07-10 03:36:40 -0400, Thomas Kowalski <th***@gmx.desaid:

>
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
Roundhouse Consulting, Ltd. (www.versatilecoding.com) Author of "The
Standard C++ Library Extensions: a Tutorial and Reference
(www.petebecker.com/tr1book)

[Pete Becker]

On 2007-07-10 03:23:52 -0400, Thomas Kowalski <th***@gmx.desaid:

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

--
-- Pete
Roundhouse Consulting, Ltd. (www.versatilecoding.com) Author of "The
Standard C++ Library Extensions: a Tutorial and Reference
(www.petebecker.com/tr1book)

[James Kanze]

On Jul 10, 9:23 am, Thomas Kowalski <[e-mail]:

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

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.

--
James Kanze (GABI Software) email:ja*********@gmail.com
Conseils en informatique orientée objet/
Beratung in objektorientierter Datenverarbeitung
9 place Sémard, 78210 St.-Cyr-l'École, France, [phone]

[kwikius]

On 9 Jul, 19:42, Juha Nieminen <[e-mail]:

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