comparing floating points

[Mortal]

hello.

first of all, apologies for the question wasn't directly related to SFML, but i was so frustrated about what's the proper way to compare floating points.  if someone had similar problem how can we do it.

i read alot of online papers about this subject. the latest one was about ULP which is used the integer comparison instead of comparing the float/double data. i've implemented this already but it seems not working either.
i made minimal example about ULPs method. which is, by the way taken from google units test. and i have changed it a little bit to suite my engine that i'm working on.
the example shows the failure of comparing two matrices for assign operator operator== when i tested the sf::Transform against A A-1 = A-1 A = I.

(click to show/hide)

#include <type_traits>
#include <cctype>
#include <cfloat>
#include <limits>
#include <bitset>
#include <iostream>
#include <iomanip>
#include <SFML/Graphics.hpp>

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

typedef Types<4>::int_type int32;
typedef Types<4>::uint_type uint32;
typedef Types<8>::int_type int64;
typedef Types<8>::uint_type uint64;

template <typename T>
class Float
{
        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);

        union
        {
                value_type value;
                bit_type bit;
                struct {
                        bit_type fraction : fraction_count;
                        bit_type exponent : exponent_count;
                        bit_type sign     : 1;
                };
        };

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

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

        enum { max_ulps = 4 };

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

        bit_type exponent_bits() const { return (exponent | 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 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);
        }

};

typedef Float<float> Float32;
typedef Float<double> Float64;

template <typename T>
struct is_float : std::false_type {};
template <typename T>
struct is_float<Float<T>> : std::true_type {};

// -- debug prints --
template <typename T> static inline
std::enable_if_t<is_float<T>::value, std::ostream&>
operator<<(std::ostream& os, const 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;
}

static inline std::ostream& operator<<(std::ostream& os, const sf::Transform& m)
{
        const float* data = m.getMatrix();
        os << std::fixed << std::setprecision(5) << std::right;
        for (unsigned int c = 0; c < 4; ++c) {
                os << "[";
                for (unsigned int r = 0; r < 4; ++r) {
                        os << std::setw(10) << data[c + 4 * r];
                }
                os << "]\n";
        }
        os << std::resetiosflags(std::ios_base::fixed | std::ios_base::floatfield);
        return os;
}

template <typename T> static inline bool
operator==(const T& m1, const T& m2)
{
        const float* data1 = m1.getMatrix();
        const float* data2 = m2.getMatrix();
        bool result = true;
        for (unsigned int i = 0; i < 16; i++) {
                result &= almost_equals(data1[i], data2[i]);
        }
        return result;
}

#define INFINITY_TEST 0
#define NAN_TEST 0
#define SIGN_TEST 0
#define ULPS_TEST 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);

        std::cout << f1 << '\n' << f2 << '\n';
#if INFINITY_TEST
        if (f1.is_infinity() || f2.is_infinity())
                return f1.data_float() == f2.data_float();
#endif
#if NAN_TEST
        if (f1.is_nan() || f2.is_nan()) return false;
#endif
#if SIGN_TEST
        if (f1.sign_bits() || f2.sign_bits())
                return f1.data_float() == f2.data_float();
#endif
#if ULPS_TEST
        return Float<T>::distance(f1.data_bits(), f2.data_bits()) <= Float<T>::max_ulps;
#endif
}

int main()
{
        sf::Transform m1(2.0f, 1.0f, 1.0f,
                                         2.0f, 1.0f, 2.0f,
                                         3.0f, 3.0f, 1.0f);
        sf::Transform inv = m1.getInverse();

        // A *A^-1 = A^-1 * A = I
        sf::Transform m2 = m1 * inv;
        bool result = (m2 == sf::Transform::Identity);

        std::cout << " matrix:\n" <<m1 << '\n'
                          << "inverse:\n" << inv << '\n'
                          << "m2 = m1 * inv:\n" << m2 << '\n'
                          << "identity:\n" << sf::Transform::Identity << '\n'
                          << "m2 == indentity : "<< std::boolalpha << result << "!\n";

        getchar();
}  

[eXpl0it3r]

Use an epsilon to compare. See also this.

[Mortal]

thanks, this will fix the matrix example. but i was looking for proper technique to compare any arbitrary float numbers that could be zero or non-zero and whether they close to each other or not. and of course it will handle over/under-flow since in floating points there are many entities like (+0-0) and (+inf-inf) and nan... etc

[Mortal]

here my humble conclusion about this subject. it seems that there are no other methods to compare floating points except the epsilon and ULPs comparisons which is used the integer representation of floating point. it works similarly like relative epsilon and both methods will fail if the difference of two values or both of them are closed to zero.

this doesn't give me choice rather than run multiple tests to single pair once with fixed epsilon and others with relative epsilon and ULPs for double checking.

here the update code after i fixed a bug with matrix example.

(click to show/hide)

#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());

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

        bool under_flow = abs_f1 < std::numeric_limits<T>::min() || abs_f2 < std::numeric_limits<T>::min();
        bool over_flow = 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 = distance < 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";
}