Comparing IEEE floats and doubles for equality

What is the best method for comparing IEEE floats and doubles for equality? I have heard of several methods, but I wanted to see what the community thought.

The best approach I think is to compare ULPs.

bool is_nan(float f)
{
    return (*reinterpret_cast<unsigned __int32*>(&f) & 0x7f800000) == 0x7f800000 && (*reinterpret_cast<unsigned __int32*>(&f) & 0x007fffff) != 0;
}

bool is_finite(float f)
{
    return (*reinterpret_cast<unsigned __int32*>(&f) & 0x7f800000) != 0x7f800000;
}

// if this symbol is defined, NaNs are never equal to anything (as is normal in IEEE floating point)
// if this symbol is not defined, NaNs are hugely different from regular numbers, but might be equal to each other
#define UNEQUAL_NANS 1
// if this symbol is defined, infinites are never equal to finite numbers (as they're unimaginably greater)
// if this symbol is not defined, infinities are 1 ULP away from +/- FLT_MAX
#define INFINITE_INFINITIES 1

// test whether two IEEE floats are within a specified number of representable values of each other
// This depends on the fact that IEEE floats are properly ordered when treated as signed magnitude integers
bool equal_float(float lhs, float rhs, unsigned __int32 max_ulp_difference)
{
#ifdef UNEQUAL_NANS
    if(is_nan(lhs) || is_nan(rhs))
    {
        return false;
    }
#endif
#ifdef INFINITE_INFINITIES
    if((is_finite(lhs) && !is_finite(rhs)) || (!is_finite(lhs) && is_finite(rhs)))
    {
        return false;
    }
#endif
    signed __int32 left(*reinterpret_cast<signed __int32*>(&lhs));
    // transform signed magnitude ints into 2s complement signed ints
    if(left < 0)
    {
        left = 0x80000000 - left;
    }
    signed __int32 right(*reinterpret_cast<signed __int32*>(&rhs));
    // transform signed magnitude ints into 2s complement signed ints
    if(right < 0)
    {
        right = 0x80000000 - right;
    }
    if(static_cast<unsigned __int32>(std::abs(left - right)) <= max_ulp_difference)
    {
        return true;
    }
    return false;
}

A similar technique can be used for doubles. The trick is to convert the floats so that they're ordered (as if integers) and then just see how different they are.

I have no idea why this damn thing is screwing up my underscores. Edit: Oh, perhaps that is just an artefact of the preview. That's OK then.

The current version I am using is this

bool is_equals(float A, float B,
               float maxRelativeError, float maxAbsoluteError)
{

  if (fabs(A - B) < maxAbsoluteError)
    return true;

  float relativeError;
  if (fabs(B) > fabs(A))
    relativeError = fabs((A - B) / B);
  else
    relativeError = fabs((A - B) / A);

  if (relativeError <= maxRelativeError)
    return true;

  return false;
}

This seems to take care of most problems by combining relative and absolute error tolerance. Is the ULP approach better? If so, why?

@DrPizza: I am no performance guru but I would expect fixed point operations to be quicker than floating point operations (in most cases).

It rather depends on what you are doing with them. A fixed-point type with the same range as an IEEE float would be many many times slower (and many times larger).

Things suitable for floats:

3D graphics, physics/engineering, simulation, climate simulation....

In numerical software you often want to test whether two floating point numbers are exactly equal. LAPACK is full of examples for such cases. Sure, the most common case is where you want to test whether a floating point number equals "Zero", "One", "Two", "Half". If anyone is interested I can pick some algorithms and go more into detail.

Also in BLAS you often want to check whether a floating point number is exactly Zero or One. For example, the routine dgemv can compute operations of the form

• y = beta*y + alpha*A*x
• y = beta*y + alpha*A^T*x
• y = beta*y + alpha*A^H*x

So if beta equals One you have an "plus assignment" and for beta equals Zero a "simple assignment". So you certainly can cut the computational cost if you give these (common) cases a special treatment.

Sure, you could design the BLAS routines in such a way that you can avoid exact comparisons (e.g. using some flags). However, the LAPACK is full of examples where it is not possible.

P.S.:

• There are certainly many cases where you don't want check for "is exactly equal". For many people this even might be the only case they ever have to deal with. All I want to point out is that there are other cases too.

• Although LAPACK is written in Fortran the logic is the same if you are using other programming languages for numerical software.