Generic and accurate floating point "equality"

Question

Like those who have come before me, I foolishly have sought to implement a generic, efficient, idiomatic, and most importantly correct method of comparing floating point numbers for equality. Like those who have come before me, I have most likely failed in at least one of those regards. I'd appreciate comments on any aspect of the following code

#include <cmath>
#include <cstddef>
#include <cstdint>
#include <limits>
#include <type_traits>
template<bool B, typename T, typename DEFAULT = void>
struct cond
{
 static constexpr bool value = B;
 using default_type = DEFAULT;
 using type = T;
};
template<typename Head, typename... Tail>
struct select
{
 using type = 
 typename std::conditional<Head::value, 
 typename Head::type, 
 typename select<Tail...>::type>::type;
};
template <typename T>
struct select<T>
{
 using type = T;
};
template <bool B, typename T, typename DEFAULT>
struct select<cond<B, T, DEFAULT>>
{
 using type = typename std::conditional<B, T, DEFAULT>::type;
};
template <typename Head, typename... Tail>
using select_t = typename select<Head, Tail...>::type;
template <typename FloatingPoint>
struct SelectIntegralOfSameSize {
 using type = select_t<
 cond<sizeof(FloatingPoint) == sizeof(int64_t), int64_t>,
 cond<sizeof(FloatingPoint) == sizeof(int32_t), int32_t>,
 cond<sizeof(FloatingPoint) == sizeof(int16_t), int16_t>,
 cond<sizeof(FloatingPoint) == sizeof(int8_t), int8_t>
 >;
 static_assert(!std::is_void<type>::value, "Could not find sane integral type");
};
template <typename FP, typename I = typename SelectIntegralOfSameSize<FP>::type>
union Fp_Int_Union {
 FP fp_;
 I i_;
};
template <typename FP>
bool nearlyEqual(FP left, FP right, FP maxDiff, std::size_t maxUlpsDiff) {
 // handle NaNs and infinities
 if (!std::isfinite(left) || !std::isfinite(right)) {
 return false;
 }
 // ULP comparison breaks for different signs
 if ((left < 0) != (right < 0)) {
 // Handle -0 == +0
 return left == right;
 }
 // ULP comparison fails near zero
 if (std::abs(left - right) <= maxDiff) {
 return true;
 }
 Fp_Int_Union<FP> l, r;
 l.fp_ = left;
 r.fp_ = right;
 return std::abs(l.i_ - r.i_) <= maxUlpsDiff;
}

• Quuxplusone: Can you explain a bit more in your question what your code is trying to do? Specifically, what would you say to someone who asked "why don't you just use the < operator?"?
• Dan Oberlam: @Quuxplusone do you mean the == operator? And I'd say that given that most real numbers cannot be exactly represented as a float or double we usually don't care about exact equality, but approximate equality.
• Deduplicator: And how do you come to your threshold?
• Dan Oberlam: @Deduplicator That'll depend on your domain and how tolerant your application is. It's hard to come up with a number that works for all situations, so I require users to pick.
• bipll: Please note that, in C++, a read from a member of a union after a write to another member is UB.

Answer 1

Consider using std::nextafter()

The math library includes the function std::nextafter() that will give you the next floating point number with the smallest possible, strictly positive difference with the given number. This is much safer than union-casting a float into an int, because there is no guarantee that two close floating point numbers will also be close as integers. In particular, that is only true if the exponent of both floating point numbers is the same. For example:

bool closeEnough = false;
if (left <= right) {
 do {
 if (left >= right) {
 closeEnough = true;
 break;
 }
 left = std::nextafter(left);
 } while (maxUlpsDiff--);
} else {
 // handle the left > right case
}
return closeEnough;

You'd obviously only want to use this method for small values of maxUlpsDiff.

Is it really what you need?

You are trying to check whether two floating point numbers are close, for two possible ways of closeness: floating point difference, and ULP distance. Is that really what you need? Are you sure that the values you choose for maxDiff and maxUlpsDiff are exactly right? If not, you might get false positives or negatives.

Can you reformulate your floating point calculations into integer calculations? Integers may have a limitted range, but at least they are exact. Maybe you can use them to perform fixed-point arithmetic, or perhaps you can use two integers to form a rational number.