Comparing doubles properly using aliasing with integer representations and ULPs

Question

I have tried to avoid epsilon comparisons for comparing floating point types. The best solution I could come up with used the difference in ULPs (Unit in the last place), though this article had a much better solution using integer representations (/// indicate my own comments):

/* See
https://randomascii.wordpress.com/2012/01/11/tricks-with-the-floating-point-format/
for the potential portability problems with the union and bit-fields below.
*/

#include <stdint.h> // For int32_t, etc.

union Float_t
{
    Float_t(float num = 0.0f) : f(num) {}
    // Portable extraction of components.
    bool Negative() const { return i < 0; }
    int32_t RawMantissa() const { return i & ((1 << 23) - 1); }
    int32_t RawExponent() const { return (i >> 23) & 0xFF; }

    int32_t i; /// Perhaps overflow when using doubles?
    float f;

    #ifdef _DEBUG
    struct
    {   // Bitfields for exploration. Do not use in production code.
        uint32_t mantissa : 23; /// 52 for double?
        uint32_t exponent : 8; /// 11 for double?
        uint32_t sign : 1;
    } parts;
    #endif
};

bool AlmostEqualUlps(float A, float B, int maxUlpsDiff)
{
    Float_t uA(A);
    Float_t uB(B);

    // Different signs means they do not match.
    if (uA.Negative() != uB.Negative())
    {
        // Check for equality to make sure +0==-0
        if (A == B)
            return true;
        return false;
    }

    // Find the difference in ULPs.
    int ulpsDiff = abs(uA.i - uB.i);
    if (ulpsDiff <= maxUlpsDiff)
        return true;

    return false;
}

However, I can't seem to reformat the code in such a way that it supports doubles. I even read up on the explanation, found here.

Does anyone know what would be the best way to tackle this?

---

Before anyone decides to mark this as a duplicate: don't, because the only question that was similar was meant for javascript, and the C++ answer was:

bool IsAlmostEqual(double A, double B)
{
    //http://www.cygnus-software.com/papers/comparingfloats/comparingfloats.htm
    long long aInt = reinterpret_cast<long long&>(A);
    if (aInt < 0) aInt = -9223372036854775808LL - aInt;
    long long bInt = reinterpret_cast<long long&>(B);
    if (bInt < 0) bInt = -9223372036854775808LL - bInt;
    return (std::abs(aInt - bInt) <= 10000);
}

Which doesn't use ULPs, but some kind of bound, and I'm not sure what -9223372036854775808LL - aInt is at all (perhaps where int64 overflows).

Answer 1

I don't think your code work at all. Below is just one example where it is going to get it wrong. (This is not an answer, but the explanation is too long to fit into a comment)

int main()
{
    Float_t x;
    Float_t y;
    x.i = 0x7F800000;
    y.i = 0x7F7FFFFF;
    AlmostEqualUlps(x.f, y.f, 1); // return true
}

However, x.f is actually infinity and y.f is FLT_MAX. The difference by any definition is infinity. Unless this is your intended behavior, that is consider a finite number and infinity being almost equal. Your implementation of ULP is completely off the mark. In fact for the two numbers above, ULP is not even well defined.

Another example will be 0x7F800000 (or any number close to this, depending on your maxULP) and 0x7F800001 (NaN). Unlike the example above, there's not even an argument to made that they should be consider "almost equal".

On the other hand, you reject any pairs with different signs as not close enough, while in truth there are a lot subnormals between -FLT_MIN and FLT_MIN, which one might considered to be "almost equal". For example, 0x80000001 and 0x1 differ by 2ULP, but if you set maxULP to 2 in your function, it will return false.

If you can rule out denormals, infintities, NaNs, then to deal with double you just need to replace uint32_t with uint64_t as mentioned in others' comment.