TITLE: Floating point equality

PROBLEM: smyth@Jpl.Nasa.Gov (Mike M. Smyth)

I have a question for language lawyers.  Is it required that the
following code

double y=some_value;
double x=y;
if(x ==y)
  cout<<"TRUE";

always prints TRUE?  I realize of course that I can do the test
fabs(x-y)<EPS, but this requires one additional floating point
operation.


RESPONSE: rmartin@rcmcon.com (Robert Martin), 18 Nov 93

It is best not to assume this.  While I would expect most compilers to
copy the double value faithfully, some might throw a "normalize" in
during the copy and lose a few insignificant bits.  In any case,
floating point expressions just a little more complex than the one
above will confound any program that depends upon == for comparison.

Also, fabs(x-y)<EPS is not a very good approach unless x and y are
guaranteed to be larger than EPS but not too much larger.  Consider if
EPS is 1e-8 but x and y are of the order of 1e40.  Even if x and y
differ only in the least significant bit, that bit will be far larger
than 1e-8 so equality of big numbers becomes impossible.

On the other hand if x and y are close to or smaller than EPS then
their difference will always be smaller than EPS no matter how wildly
their bit patterns differ.  Thus small numbers will always compare
equal.

The following C++ module is my solution to the problem of floating
point equality.   I don't contend that it is the most efficient
solution, but it works for my purposes.

[ RCM code moved to end of article. -adc ]


RESPONSE: nagle@netcom.com (John Nagle)

     No.  There is an obscure exception.  For machines with IEEE standard 
floating point (this includes almost all modern computers), there is a
floating point value called "NaN", for "not a number".  There are also
values "infinity" and "negative infinity".  

     If "some_value" is "not a number",  then "x ==y" will be FALSE.

     If you want to play with these, use the "Calculator" desk accessory
on a Macintosh, and try things like "0 / 0", which will yield "not a
number".  Also try "1 / 0" and "-1 / 0", to get "infinity" and "-infinity".

     That [ testing fabs(x-y)<EPS ] should raise an exception in an IEEE
 environment, depending on what math interrupts are enabled.  It is an
 error to apply "<" to "not a number" in some implementations.  "not a number"
is considered unordered, and ordered comparisons are prohibited.

     If you have floating point exceptions enabled, so that you get an
interrupt rather than a "silent NAN" for something like "0.0 / 0.0",
then you should not get NAN values in variables as a result of computation.
But if you don't handle floating point exceptions, NAN and INFINITY values
can appear and propagate through the computation, and you have to handle
them.

     A good explanation of this subject appears in the "Apple Numerics
Manual, 2nd. edition", on page 48.


RESPONSE: matt@physics16.berkeley.edu (Matt Austern)

Equality testing of floating-point numbers is dangerous.  I can't think
of anything in the standard that guarantees that == will ever return 
true, if its arguments are floating-point.

-----------------------------------------------------------------------------
	RCM code
-----------------------------------------------------------------------------

------------------------------feq.h---------------------------------
#ifndef COMPONENT_FEQ_H
#define COMPONENT_FEQ_H

int Feq(double, double, double epsilon=.000001);
int Fle(double, double, double epsilon=.000001);
int Fge(double, double, double epsilon=.000001);

#endif
--------------------------------------------------------------------

----------------------------feq.cpp---------------------------------

#include "componen/feq.h"
template <class T> T abs(const T& a) {return (a<0) ? -a : a;}

//---------------------------------------------
// Name
//  Feq
//
// Description
//  Used to test two floating point numbers to see if
//  they are close enough to be considered equal.  The
//  usual espsilon is 1e-6
//
int Feq(double a, double b, double e)
{
	int retval;
	if (a == 0 && abs(b) < e)
		retval = 1;
	else if (b == 0 && abs(a) < e)
		retval = 1;
	else
	{
		retval = (abs((a-b)) < (e * abs(a)));
        }
	return retval;
}

int Fle(double a, double b, double e)
{
    return ( (a<b) || Feq(a,b,e));
}

int Fge(double a, double b, double e)
{
	return ( (a>b) || Feq(a,b,e));
}
-------------------------------------------------------------------