2
\$\begingroup\$

I have designed a class to perform float compare. I used Knuth's strong compare to do it. Float compare is really really tricky, I'm sure there's something wrong :) Any comments?

#include <limits>
#include <cmath>
#include <assert.h>

class FloatUtils {
private:
    FloatUtils();
    template<typename T>
    static T safeDivision( const T& f1, const T& f2 );
public:
    template<typename T>
    static bool almostEqual( const T& a, const T& b, const T& relEpsilon = std::numeric_limits<T>::epsilon(),const T& absEpsilon = std::numeric_limits<T>::epsilon()) noexcept;
    template<typename T>
    static bool almostZero( const T& a, const T& absEpsilon = std::numeric_limits<T>::epsilon() ) noexcept;
};

template<typename T>
inline T FloatUtils::safeDivision( const T& f1, const T& f2 ) {
    // Avoid overflow.
    if ( (f2 < static_cast<T>(1)) && (f1 > f2 * std::numeric_limits<T>::max()) )
        return std::numeric_limits<T>::max();

    // Avoid underflow.
    if ( (f1 == static_cast<T>(0)) || ((f2 > static_cast<T>(1)) && (f1 < f2 * std::numeric_limits<T>::min())) )
        return static_cast<T>(0);

    return f1 / f2;
}

template<typename T>
inline bool FloatUtils::almostEqual( const T& a, const T& b, const T& relEpsilon, const T& absEpsilon ) noexcept {
    static_assert(!std::is_integral<T>::value, "T can be only no integral type");
    static_assert(std::numeric_limits<T>::is_iec559, "IEEE754 supported only");
  
    if ( !std::isfinite(a) || !std::isfinite(b) ) {
        return false;
    }

    /**
     * avoid division by zero checking exactly 0 value
     */
    if ( a == 0 ) {
        return almostZero(b, absEpsilon);
    }

    /**
     * avoid division by zero checking exactly 0 value
     */
    if ( b == 0 ) {
        return almostZero(a, absEpsilon);
    }

    /**
     * We use Knuth's formula modified, i.e. the right side of inequations may cause underflow,
     * to prevent this, the formula scales the difference rather than epsilon.
     * We don't nanage 0 or Inf here so we need to add the special cases above.
     */
    T diff = std::abs(a - b);
    T d1 = safeDivision<T>(diff, std::abs(a));
    T d2 = safeDivision<T>(diff, std::abs(b));

    return d1 <= relEpsilon && d2 <= relEpsilon;
}

template<typename T>
inline bool FloatUtils::almostZero( const T& a, const T& absEpsilon ) noexcept {
    if ( std::abs(a) <= absEpsilon )
        return true;
    return false;
}

Test cases:

int main() {
assert(FloatUtils::almostEqual(3.145646, 3.145646) == true);
assert(FloatUtils::almostEqual(145697.78965, 145553.09186, 0.001) == true);
assert(FloatUtils::almostEqual(145697.78965, 145552.09186, 0.001) == false);
assert(FloatUtils::almostEqual(0.00001, 0.0, 0.1, 0.001) == true);
assert(FloatUtils::almostEqual(1.0, 1.0, 0.0) == true);
assert(FloatUtils::almostEqual(0.0, 1E-20, 1E-5) == true);
assert(FloatUtils::almostEqual(0.0, 1E-30, 1E-5) == true);
assert(FloatUtils::almostEqual(0.0, -1E-10, 0.1, 1E-9) == true);
assert(FloatUtils::almostEqual(0.0, -1E-10, 0.1, 1E-11) == false);
assert(FloatUtils::almostEqual(0.123456, 0.123457, 1E-6) == false);
assert(FloatUtils::almostEqual(0.123456, 0.123457, 1E-3) == true);
assert(FloatUtils::almostEqual(0.123456, -0.123457, 1E-3) == false);
assert(FloatUtils::almostEqual(1.23456E28, 1.23457E28, 1E-3) == true);
assert(FloatUtils::almostEqual(1.23456E-10, 1.23457E-10, 1E-3) == true);
assert(FloatUtils::almostEqual(1.111E-10, 1.112E-10, 0.000899) == false);
assert(FloatUtils::almostEqual(1.111E-10, 1.112E-10, 0.1) == true);
assert(FloatUtils::almostEqual(1.0, 1.0001, 1.1E-2) == true);
assert(FloatUtils::almostEqual(1.0002, 1.0001, 1.1E-2) == true);
assert(FloatUtils::almostEqual(1.0, 1.0002, 1.1E-4) == false);
return 0;
}
\$\endgroup\$
0

1 Answer 1

5
\$\begingroup\$

First, note that anyone trying to "fuzzy compare" floating-point numbers for "approximate equality" is already in a state of sin. In real code, either you care about actual (bit-for-bit) equality, or else you have a specific error-bar in mind which you've carried through the whole calculation. There's pretty much no situation where the programmer ever wants "x is approximately equal to y, based on some fudge-factor epsilons I don't even care enough to specify." That's just not a thing people want to do.


class FloatUtils {
private:
    FloatUtils();

You're using a struct instead of a namespace because you don't like namespaces, right? But then you're also giving it a private, non-defined constructor so that the user can't accidentally create an instance of the struct type itself? First of all, this is massive overkill and you should just write either

namespace FloatUtils {

or

struct FloatUtils {

But, if you really must prevent people from constructing instances of your type, the correct way to do that is not with a private constructor but instead with a deleted constructor:

struct FloatUtils {
    FloatUtils() = delete;

I don't see why you need safeDivision to be private. Couldn't it profitably be made public?


template<typename T>
static bool almostEqual( const T& a, const T& b, const T& relEpsilon = std::numeric_limits<T>::epsilon(),const T& absEpsilon = std::numeric_limits<T>::epsilon()) noexcept;

This line is comically long, and also default function arguments are the devil. It would be strictly better to write one declaration per signature to which you intend the user to have access:

template<class T>
static bool almostEqual(const T& a, const T& b) {
    auto eps = std::numeric_limits<T>::epsilon();
    return almostEqual(a, b, eps, eps);
}

template<class T>
static bool almostEqual(const T& a, const T& b, const T& relEps) {
    return almostEqual(a, b, relEps, std::numeric_limits<T>::epsilon());
}

template<class T>
static bool almostEqual(const T& a, const T& b, const T& relEps, const T& absEps);

Marking these functions noexcept doesn't serve any purpose; nobody is going to be making compile-time decisions based on whether noexcept(almostEqual(a, b)).


It is mildly strange that you consider almostEqual(INFINITY, INFINITY) == false. They're bitwise the same value; shouldn't they be considered equal?


static_assert(!std::is_integral<T>::value, "T can be only no integral type");

I suspect you mean

static_assert(std::is_floating_point_v<T>, "T must be a floating-point type");

Why do you need safeDivision to protect against overflow? It might well be that abs(a - b) is INFINITY, and so diff / abs(a) is also INFINITY... but that's definitely going to be greater than absEpsilon, so who cares?

Vice versa, I suggest that if a > 0 && b < 0, the numbers are definitely not approximately equal... well, not for many applications, anyway. Depends if you're using them as operands to + or as operands to *, I guess. Again, anyone comparing floats for fuzzy equality is already in a state of sin.


T diff = std::abs(a - b);
T d1 = safeDivision<T>(diff, std::abs(a));
T d2 = safeDivision<T>(diff, std::abs(b));

return d1 <= relEpsilon && d2 <= relEpsilon;

It seems like you could save a division by doing

T diff = std::abs(a - b);
T denominator = std::min(std::abs(a), std::abs(b));
return safeDivision(diff, denominator) <= relEpsilon;

Notice that I've dropped the T from safeDivision. In general, you should call function templates as if they were ordinary functions, without any scary angle brackets. Let type deduction do its thing.

If you really want to force the programmer to specify the T argument to safeDivision, then you should make it non-deducible. In C++20 you can do this with type_identity_t:

template<class T>
static T safeDivision(const std::type_identity_t<T>& f1, const std::type_identity_t<T>& f2);

safeDivision(1.0, 2.0);  // error, can't deduce T
safeDivision<double>(1.0, 2.0);  // OK
safeDivision<float>(1.0, 2.0);  // also OK

By the way, since T must be a floating-point type, you're losing performance via unnecessary pass-by-reference. Prefer pass-by-value for trivial primitive types like float and long double:

template<class T>
static T safeDivision(T a, T b);
\$\endgroup\$

You must log in to answer this question.

Start asking to get answers

Find the answer to your question by asking.

Ask question

Explore related questions

See similar questions with these tags.