Floating point comparison program

Question

I made a small program to compare floating-point values by running multiple tests to a single pair of floats. I have added extra checks for under/overflow and allowing tolerance.

How can I improve it further?

#include <type_traits>
#include <cctype>
#include <cfloat>
#include <limits>
#include <bitset>
#include <iostream>
#include <iomanip>
#include <algorithm>

template <size_t size>
struct Types {
    typedef void int_type;
};

template <>
struct Types<4> {
    typedef int int_type;
    typedef unsigned int uint_type;
};

template <>
struct Types<8> {
    typedef __int64 int_type;
    typedef unsigned __int64 uint_type;
};

template <typename T>
class Float
{
public:
    typedef typename Types<sizeof(T)>::uint_type bit_type;
    typedef typename T value_type;

    static const bit_type bit_count = 8 * sizeof(value_type);
    static const bit_type fraction_count = std::numeric_limits<value_type>::digits - 1;
    static const bit_type exponent_count = bit_count - 1 - fraction_count;

    static const bit_type sign_mask = static_cast<bit_type>(1) << (bit_count - 1);
    static const bit_type fraction_mask = ~static_cast<bit_type>(0) >> (exponent_count + 1);
    static const bit_type exponent_mask = ~(sign_mask | fraction_mask);
    static const bit_type max_ulps = static_cast<bit_type>(4);

    explicit Float(const T& x) { value = x; }

    const value_type &data_float() const { return value; }
    const bit_type &data_bits() const { return bit; }

    bit_type exponent_bits() const { return (exponent_mask & bit); }
    bit_type sign_bits() const { return sign; }
    bit_type fraction_bits() const { return fraction; }

    bool is_infinity()
    {
        return ((bit & ~sign_mask) == exponent_mask);
    }

    bool is_nan() const {
        bool nan = true;
        nan &= (exponent_mask & bit) == exponent_mask;
        nan &= (fraction_mask & bit) != static_cast<bit_type>(0);
        return nan;
    }

    static bit_type to_biased(bit_type bits) {
        return (sign_mask & bits) ? (~bits + 1) : (sign_mask | bits);
    }

    static bit_type distance(bit_type bits1, bit_type bits2) {
        const bit_type biased1 = to_biased(bits1);
        const bit_type biased2 = to_biased(bits2);
        return (biased1 >= biased2) ? (biased1 - biased2) : (biased2 - biased1);
    }

private:
    union
    {
        value_type value;
        bit_type bit;

        struct {
            bit_type fraction : fraction_count;
            bit_type exponent : exponent_count;
            bit_type sign : 1;
        };
    };
};

#define PRINT_DEBUG_INFO 1

template <typename T> static inline
std::enable_if_t<std::is_floating_point<T>::value, bool>
almost_equals(const T& lhs, const T& rhs)
{
    Float<T> f1(lhs), f2(rhs);

#if PRINT_DEBUG_INFO
    std::cout << std::setfill(' ') << f1 << '\n' << f2;
    std::cout << std::setfill('-') << std::setw(71) << ' ' << std::setfill(' ') << std::endl;
#endif

    const Float<T>::bit_type distance = Float<T>::distance(f1.data_bits(), f2.data_bits());
    const T abs_f1 = std::abs(f1.data_float());
    const T abs_f2 = std::abs(f2.data_float());
    const T diff = std::max(std::abs(abs_f1 - abs_f2), std::numeric_limits<T>::min());
    const T sum = std::min(std::abs(abs_f1 + abs_f2), std::numeric_limits<T>::max());

    const T tolerance = static_cast<T>(0.000001);

    bool under_flow = diff < std::numeric_limits<T>::min() || abs_f1 < std::numeric_limits<T>::min() || abs_f2 < std::numeric_limits<T>::min();
    bool over_flow = diff > std::numeric_limits<T>::max() || abs_f1 > std::numeric_limits<T>::max() || abs_f2 > std::numeric_limits<T>::max();

    bool sign = (f1.sign_bits() ^ f2.sign_bits()) == 1 && !(under_flow ^ over_flow);

    bool inff = (f1.is_infinity() ^ f2.is_infinity()) == 1;
    bool nan = (f1.is_nan() ^ f2.is_nan()) == 1;

    bool assign = f1.data_float() == f2.data_float();
    bool ulp = Float<T>::distance(f1.data_bits(), f2.data_bits()) < Float<T>::max_ulps;
    bool fixed_epsilon = diff < tolerance;
    bool relative_epsilon = diff < std::numeric_limits<T>::epsilon() * sum;

#if PRINT_DEBUG_INFO
    std::cout << "\n"
        << "distance         = " << distance << '\n'
        << "diff             = " << diff << '\n'
        << "sum              = " << sum << '\n'
        << "min              = " << std::numeric_limits<T>::min() << '\n'
        << "max              = " << std::numeric_limits<T>::max() << '\n'
        << "---------------- \n"
        << std::boolalpha
        << std::setw(15) << "under_flow       = " << std::setw(7) << under_flow << '\n'
        << std::setw(15) << "over_flow        = " << std::setw(7) << over_flow << '\n'
        << std::setw(15) << "diff sign        = " << std::setw(7) << sign << '\n'
        << std::setw(15) << "inf              = " << std::setw(7) << inff << '\n'
        << std::setw(15) << "nan              = " << std::setw(7) << nan << '\n'
        << "---------------- \n"
        << std::setw(15) << "assign           = " << std::setw(7) << assign << '\n'
        << std::setw(15) << "ulp              = " << std::setw(7) << ulp << '\n'
        << std::setw(15) << "fixed_epsilon    = " << std::setw(7) << fixed_epsilon << '\n'
        << std::setw(15) << "relative_epsilon = " << std::setw(7) << relative_epsilon << "\n\n";
    std::cout << std::setfill('-') << std::setw(71) << ' ' << std::setfill(' ') << "\n\n";
#endif

    if (sign || nan || inff) return false;

    return assign || ulp || fixed_epsilon || relative_epsilon;
}

// -- debug prints --
template <typename T> static inline std::ostream&
operator<<(std::ostream& os, const Float<T>& f)
{
    os << std::fixed << std::setprecision(25) << std::left;

    os << "float    = " << std::setw(10) << std::dec << f.data_float() << "\n";

    os << "bits     = " << std::setw(10) << std::dec << f.data_bits() << " : 0x"
        << std::setw(10) << std::hex << f.data_bits() << " : "
        << std::setw(32) << std::bitset<32>(f.data_bits()) << "\n";

    os << "sign     = " << std::setw(10) << std::dec << f.sign_bits() << " : 0x"
        << std::setw(10) << std::hex << f.sign_bits() << " : "
        << std::setw(32) << std::bitset<32>(f.sign_bits()) << "\n";

    os << "exponent = " << std::setw(10) << std::dec << f.exponent_bits() << " : 0x"
        << std::setw(10) << std::hex << f.exponent_bits() << " : "
        << std::setw(32) << std::bitset<32>(f.exponent_bits()) << "\n";

    os << "fraction = " << std::setw(10) << std::dec << f.fraction_bits() << " : 0x"
        << std::setw(10) << std::hex << f.fraction_bits() << " : "
        << std::setw(32) << std::bitset<32>(f.fraction_bits()) << "\n\n";

    os << std::resetiosflags(std::ios_base::fixed | std::ios_base::floatfield) << std::dec;
    return os;
}

int main()
{
    float a = 1.0f;
    float b = 3.0f;
    bool result = almost_equals(a, b);
    std::cout << "a == b : " << std::boolalpha << result << "!\n\n";

    float c = a / b;
    float d = b / a;
    result = almost_equals(c,1/d);
    std::cout << "a/b == a/b : " << std::boolalpha << result << "!\n\n";
}

Answer 1

A good effort with lots of useful debug to see the inner working of the task.

---

Biggest design concern: lack of control on how "almost" almost_equals() operates. I'd expect a 3rd, perhaps optional, parameter to gauge "closeness".

From a design stand-point, code is mixing absolute and relative in determining the the return of almost_equals() and is incapable of returning a sensible result based on just one of those. Perhaps 2 functions:

almost_equals_rel(value_type a, value_type b, uint_type ulp_abs_difference);
almost_equals_abs(value_type a, value_type b, value_type max_abs_difference);

---

This and that:

1. Code uses a fixed tolerance, yet uses templated code. A value of 0.000001 is quite arbitrary and not justified in code. Might make sense for float, yet not double. IMO, the tolerance assessment detracts from this code.

   // Why the magic number 0.000001?
   const T tolerance = static_cast<T>(0.000001);
   

2. -0.0 has the same value of 0.0, yet incorrectly fails almost_equals() due to differ signs.

3. Values of 7, 10 and 15 throughout code should 1) be driven by a common named constants and 2) not a magic number. IMO, for typical float and double I would expect 9 and 17 significant digits printed to allow round-tripping the text back to the FP variable. (Of course this is mostly debug code.)

4. bool nan = (f1.is_nan() ^ f2.is_nan()) == 1; goes against the spirit of IEEE math. Even if both f1,f2 are both NaNs with the same bit pattern, they are not equal.

5. Unsigned fixed width types makes sense for double. I'd use int64_t int_type rather than __int64 int_type Potable code should use fixed width for struct Types<4>

   struct Types<4> {
     //typedef int int_type;
     //typedef unsigned int uint_type;
     typedef int32_t int_type;
     typedef uint32_t uint_type;
   };
   

6. For code to function correctly, it assumes the endian of integers and FP is the same. Certainly common, but not specified and exceptions exists.

   union {
     value_type value;
     bit_type bit;
   

7. Further, code assumes the bit packing/bit fields are in the desired order.

   struct {
       bit_type fraction : fraction_count;
       bit_type exponent : exponent_count;
       bit_type sign : 1;
   };
   // or should it be  
   struct {
       bit_type sign : 1;
       bit_type exponent : exponent_count;
       bit_type fraction : fraction_count;
   };
   

8. A comment that applies to C, unsure about C++: bit fields are well defined for integer types unsigned, signed int and maybe bool. So a bit field of type __int64 may lack portability. Higher portable code would use shifts and masks rather than bit fields to control endian, range and padding issues.

9. int main() only exercises float. Recommend adding a double example. Also compare values a,b that are both about 1e30 or 1e-30. The test set is much to small.