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

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.