sanity check - floating point comparison

M

ma740988

template <class T> inline bool isEqual( const T& a, const T& b,
const T epsilon = std::numeric_limits<T>::epsilon() )
{
const T diff = a - b;
return ( diff <= epsilon ) && ( diff >= -epsilon );
}

int main()
{
std::deque<double> pt ;
pt.push_back ( 2.3123 );
pt.push_back ( 4.3445 );
pt.push_back ( 1.343 );
pt.push_back ( 4.3445 );
pt.push_back ( 4.3445 );
pt.push_back ( 2.3123 );


std::deque<double> jt ;
jt.push_back ( 4.3445 );
jt.push_back ( 2.3123 );

std::deque<double> results;
// results should have - 1, 3, 4, 0, 5

for ( int idx ( 0 ) ; idx < jt.size(); ++idx )
{
for ( int kdx ( 0 ); kdx < pt.size(); ++kdx )
{
if ( isEqual<double>( jt [ idx ], pt [ kdx ] ) )
{
results.push_back ( kdx );
}
}
}
std::copy ( results.begin(), results.end(),
std::eek:stream_iterator<int> ( std::cout, "\n" ) );
}

My intent. I'll search the pt container for the values in the jt
container. If found, I'll store - in the result container - the
index/location where the value was found in the pt container. For
example: Search pt container for first element in jt container. So I
found 4 at postions 1, 3 and 4. Similarily for 2. 2 was found within
the pt container at 0 and 5.

Result will print, 1, 3, 4, 0, 5. Works.

What makes me nervous though is the floating point comparsion. After
all numeric_limits is not defined on one platform (using gcc 2.96).
That said, I was opting to use iterators with - I think std::distance
but I'll end up doing floating point comparison anyways, in which case
my own version ( with my own comparator - isEqual) works best. Correct?
 
R

Richard Herring

template <class T> inline bool isEqual( const T& a, const T& b,
const T epsilon = std::numeric_limits<T>::epsilon() )
{
const T diff = a - b;
return ( diff <= epsilon ) && ( diff >= -epsilon );
}

int main()
{
std::deque<double> pt ;
pt.push_back ( 2.3123 );
pt.push_back ( 4.3445 );
pt.push_back ( 1.343 );
pt.push_back ( 4.3445 );
pt.push_back ( 4.3445 );
pt.push_back ( 2.3123 );


std::deque<double> jt ;
jt.push_back ( 4.3445 );
jt.push_back ( 2.3123 );

std::deque<double> results;
// results should have - 1, 3, 4, 0, 5

for ( int idx ( 0 ) ; idx < jt.size(); ++idx )
{
for ( int kdx ( 0 ); kdx < pt.size(); ++kdx )
{
if ( isEqual<double>( jt [ idx ], pt [ kdx ] ) )
{
results.push_back ( kdx );
}
}
}
std::copy ( results.begin(), results.end(),
std::eek:stream_iterator<int> ( std::cout, "\n" ) );
}

My intent. I'll search the pt container for the values in the jt
container. If found, I'll store - in the result container - the
index/location where the value was found in the pt container. For
example: Search pt container for first element in jt container. So I
found 4 at postions 1, 3 and 4. Similarily for 2. 2 was found within
the pt container at 0 and 5.

Result will print, 1, 3, 4, 0, 5. Works.

What makes me nervous though is the floating point comparsion. After
all numeric_limits is not defined on one platform (using gcc 2.96).
?

That said, I was opting to use iterators with - I think std::distance
Click to expand...

? Use of iterators versus indexing is a completely different question
from how to compare the stored values.
but I'll end up doing floating point comparison anyways, in which case
my own version ( with my own comparator - isEqual) works best. Correct?
Click to expand...
IMO no.

I can't imagine anyone deliberately setting out to produce a compiler on
which the same floating literal 2.3123 etc., wouldn't produce the same
double value internally in each place in the code where it's used, if
you use the same compiler flags. (If you're generating the numbers by
arithmetic it's a different matter, of course.)

Even if the compiler did produce different values for the same literal,
there's no reason to suppose that the results would pass isEqual().
epsilon() has its uses in numerical analysis, but I don't think this is
an appropriate one. In effect it's guaranteed that if X==1.0, (a) X and
X + epsilon() are distinct values and (b) X and X + epsilon()/2 are not,
but for other values of X one of those assertions may not be true.

To specify fuzzy floating-point comparisons correctly you generally need
some knowledge of the domain being modelled, but here you're trying to
second-guess a compiler problem that probably doesn't even exist.

Finally, if you really must use your own comparator, don't call it
isEqual.
 
M

ma740988

IMO no.

I can't imagine anyone deliberately setting out to produce a compiler on
which the same floating literal 2.3123 etc., wouldn't produce the same
double value internally in each place in the code where it's used, if
you use the same compiler flags. (If you're generating the numbers by
arithmetic it's a different matter, of course.)
Click to expand...
So if I understand you correctly. Simply doing
if ( jt [ idx ] == pt [ kdx ] ) )
would suffice?
 
B

Bernd Strieder

Hello,
template <class T> inline bool isEqual( const T& a, const T& b,
const T epsilon = std::numeric_limits<T>::epsilon() )
{
const T diff = a - b;
return ( diff <= epsilon ) && ( diff >= -epsilon );
}
Click to expand...


That epsilon is the difference of 1 and the least value strictly greater
than 1. If you take 1000000 and the least value strictly greater than
1000000, than you will see, that the difference is greater than
epsilon. I don't think you want those values to make isEqual return
false.

int main()
{
std::deque<double> pt ;
pt.push_back ( 2.3123 );
pt.push_back ( 4.3445 );
pt.push_back ( 1.343 );
pt.push_back ( 4.3445 );
pt.push_back ( 4.3445 );
pt.push_back ( 2.3123 );


std::deque<double> jt ;
jt.push_back ( 4.3445 );
jt.push_back ( 2.3123 );
Click to expand...
if ( isEqual<double>( jt [ idx ], pt [ kdx ] ) )
Click to expand...
my own version ( with my own comparator - isEqual) works best.
Correct?
Click to expand...

Almost, just don't abuse epsilon. Your input values seem to be rounded
to 4 digits after decimal point, so take 1e-4 instead of epsilon. You
will have to add epsilon as parameter to isEqual, either as template
argument via some bypass, or in the function argument, or hardcode
epsilon some other way.

Bernd Strieder
 
D

dan2online

ma740988 said:
What makes me nervous though is the floating point comparsion. After
all numeric_limits is not defined on one platform (using gcc 2.96).
Click to expand...

You need to upgrade your complier to gcc 3.x at least. Or you try to
find a patch for gcc 2.96 to fix the problem.
 
D

dan2online

Richard said:
To specify fuzzy floating-point comparisons correctly you generally need
some knowledge of the domain being modelled, but here you're trying to
second-guess a compiler problem that probably doesn't even exist.
Click to expand...
why not ? But it is a fact that gcc 2.96 has a problem in handling
"epsilon" .
 
D

dan2online

ma740988 said:
IMO no.

I can't imagine anyone deliberately setting out to produce a compiler on
which the same floating literal 2.3123 etc., wouldn't produce the same
double value internally in each place in the code where it's used, if
you use the same compiler flags. (If you're generating the numbers by
arithmetic it's a different matter, of course.)
Click to expand...
So if I understand you correctly. Simply doing
if ( jt [ idx ] == pt [ kdx ] ) )
would suffice?
Click to expand...

It doesn't work for float point number.
 
D

dan2online

Bernd said:
That epsilon is the difference of 1 and the least value strictly greater
than 1. If you take 1000000 and the least value strictly greater than
1000000, than you will see, that the difference is greater than
epsilon. I don't think you want those values to make isEqual return
false.
Click to expand...

epsilon depends on data type.

Almost, just don't abuse epsilon. Your input values seem to be rounded
to 4 digits after decimal point, so take 1e-4 instead of epsilon. You
will have to add epsilon as parameter to isEqual, either as template
argument via some bypass, or in the function argument, or hardcode
epsilon some other way.
Click to expand...

the hard code should work here. But it is not the best way that uses
epsilon.
If the data type is not double but other kind of struct, how do you
deal with it?

struct new_type
{
char x;
int y;
long z;
double w;
}
 
R

Richard Herring

IMO no.

I can't imagine anyone deliberately setting out to produce a compiler on
which the same floating literal 2.3123 etc., wouldn't produce the same
double value internally in each place in the code where it's used, if
you use the same compiler flags. (If you're generating the numbers by
arithmetic it's a different matter, of course.)
Click to expand...
So if I understand you correctly. Simply doing
if ( jt [ idx ] == pt [ kdx ] ) )
would suffice?
Click to expand...

It doesn't work for float point number.
Click to expand...

What do you mean "it doesn't work"? Equality is perfectly well defined
for floating types. It just isn't always what you want to test.
 
R

Richard Herring

epsilon depends on data type.
Click to expand...

Regardless, it's still the wrong thing to use, because it only
approximates "least detectable difference" for values near to 1.0.
the hard code should work here. But it is not the best way that uses
epsilon.
If the data type is not double but other kind of struct, how do you
deal with it?

struct new_type
{
char x;
int y;
long z;
double w;
}
Click to expand...

Then you understand the problem domain it's supposed to be modelling,
and define your comparison accordingly.
 
R

Richard Herring

IMO no.

I can't imagine anyone deliberately setting out to produce a compiler on
which the same floating literal 2.3123 etc., wouldn't produce the same
double value internally in each place in the code where it's used, if
you use the same compiler flags. (If you're generating the numbers by
arithmetic it's a different matter, of course.)
Click to expand...
So if I understand you correctly. Simply doing
if ( jt [ idx ] == pt [ kdx ] ) )
would suffice?
Click to expand...

Yes, *if* the numbers were floating literals inserted into jt and pt as
your previous code suggests.

If they are computed, you need a comparator named according to what it
does e.g. IsNearlyEqual() which compares with a tolerance you determine,
based on the nature of the domain you're trying to model.
 
R

Richard Herring

My frame of reference was the C++ FAQ. In that regard I tried to model
something akin to:

http://www.parashift.com/c++-faq-li....html#faq-29.17
In that regard I was just trying to get 'close' as the FAQ puts it.
Click to expand...
Fair enough, but note that their "epsilon" is not the function from
numeric_limits, but "some small number such as 1e-5" and the onus is on
you to determine what is appropriate. Moreover, in the FAQ code, epsilon
represents a relative error, not an absolute one.

The worst thing about that FAQ is the name of the function. It doesn't
compute equality, so isEqual() is not an appropriate name for it - as
the author points out, it isn't even symmetric.
 
M

Marcus Kwok

Richard Herring said:
Fair enough, but note that their "epsilon" is not the function from
numeric_limits, but "some small number such as 1e-5" and the onus is on
you to determine what is appropriate. Moreover, in the FAQ code, epsilon
represents a relative error, not an absolute one.

The worst thing about that FAQ is the name of the function. It doesn't
compute equality, so isEqual() is not an appropriate name for it - as
the author points out, it isn't even symmetric.
Click to expand...

This whole discussion reminds me of an article I read a while ago called
"Comparing Floats - How To Determine if Floating Quantities Are Close
Enough Once a Tolerance Has Been Reached" in the March 2000 C++ Report.
It used to be online at http://www.adtmag.com/joop/crarticle.asp?ID=396
but not anymore. In the process of trying to find the above-cited
article, I found this, which may be helpful for the OP:

http://www.boost.org/libs/test/doc/components/test_tools/floating_point_comparison.html
 
D

dan2online

Richard said:
In message <[email protected]>,
dan2online <[email protected]> writes
Click to expand...
Regardless, it's still the wrong thing to use, because it only
approximates "least detectable difference" for values near to 1.0.
Click to expand...

First, you can test the example of the orignal post and print out
epsilon like:


template <class T> inline bool isEqual( const T& a, const T& b,
const T epsilon = std::numeric_limits<T>::epsilon() )
{
const T diff = a - b;

+ cout << epsilon <<endl;
return ( diff <= epsilon ) && ( diff >= -epsilon );


}

so you will see what value the "epsilon" is.
epsilon = 2.22045e-016 for IA32 architecture.

Second, "only approximates "least detectable difference" for values
near to 1.0" doesn't mean that epsilon is near to 1.0. You can test
the following code for the epsilon algorithm.

#include <iostream>
int main()
{
double base = 1.0f;
double test = base;
double epsilon = test;

while (1)
{
double tmp = base + test;
if (tmp == base)
break;
else
epsilon = test;
test /= 2.0f;
}

std::cout << "epsilon = "<<epsilon << "\n";
return 0;
}


Third, the original post uses gcc 2.96 that cannot handle the epsilon
properly.
That is a defect with gcc 2.96.
 
D

dan2online

What do you mean "it doesn't work"? Equality is perfectly well defined
for floating types. It just isn't always what you want to test.
Click to expand...

Only if two floating point numbers have the same bit pattern in memory,
you can say they are *perfectly* equal.
 
P

Phlip

dan2online said:
Only if two floating point numbers have the same bit pattern in memory,
you can say they are *perfectly* equal.
Click to expand...

RH said "perfectly well defined", not "perfectly equal".

I don't know if The Standard says two floats with the same bit pattern must
compare equal. I will not rely on that, and I can envision a CPU with an
arithmetic logic unit that optimizes such a comparison in some inscrutable
way that breaks equality. So I would prefer my compiler to produce fast
opcodes, not opcodes that hold my hand.
 
R

Richard Herring

First, you can test the example of the orignal post and print out
epsilon like:


template <class T> inline bool isEqual( const T& a, const T& b,
const T epsilon = std::numeric_limits<T>::epsilon() )
{
const T diff = a - b;

+ cout << epsilon <<endl;
return ( diff <= epsilon ) && ( diff >= -epsilon );


}

so you will see what value the "epsilon" is.
epsilon = 2.22045e-016 for IA32 architecture.
Click to expand...

Yes. So what?
Second, "only approximates "least detectable difference" for values
near to 1.0" doesn't mean that epsilon is near to 1.0.
Click to expand...

You appear to be reading things I never posted. Of course epsilon is
nowhere near 1.0, it's *the values being compared* that have to be near
1.0

[...]
Third, the original post uses gcc 2.96 that cannot handle the epsilon
properly.
That is a defect with gcc 2.96.
Click to expand...
Since numeric_limits::epsilon() is the wrong quantity to use here, why
does that matter?
 

Ask a Question

Want to reply to this thread or ask your own question?

You'll need to choose a username for the site, which only take a couple of moments. After that, you can post your question and our members will help you out.

Ask a Question

Similar Threads

Floating Point comparison problem 26
Help optimize nbody bench program (c++ sse2 intrinsics) 3
IEEE floating point 2
Portable marshalling of floating point data 3
Is there such thing as invalid floating point ? 4
String comparison? 1
Somone's SO question: "Is there an existing library for dynamically-determineddimensional array in c 1
Indexed list 0

Members online

No members online now.
Total: 33 (members: 1, guests: 32)
Robots: 402

Forum statistics

Threads
473,768
Messages
2,569,574
Members
45,051
Latest member
CarleyMcCr

Latest Threads

Top