Revisit use of FLT_EPSILON

[phillipjohnston]

Upon further research, FLT_EPSILON is fraught with problems and not useful comparisons on numbers < 1 or > 2, and I should probably take a different approach to float comparisons.

[phillipjohnston]

https://bitbashing.io/comparing-floats.html

[phillipjohnston]

ULP might be suitable, example from above article. Probably also want this in the embvm-core as a utility.

int32_t ulpsDistance(const float a, const float b)
{
    // We can skip all the following work if they're equal.
    if (a == b) return 0;

    const auto max = numeric_limits<int32_t>::max();

    // We first check if the values are NaN.
    // If this is the case, they're inherently unequal;
    // return the maximum distance between the two.
    if (isnan(a) || isnan(b)) return max;

    // If one's infinite, and they're not equal,
    // return the max distance between the two.
    if (isinf(a) || isinf(b)) return max;

    // At this point we know that the floating-point values aren't equal and
    // aren't special values (infinity/NaN).
    // Because of how IEEE754 floats are laid out
    // (sign bit, then exponent, then mantissa), we can examine the bits
    // as if they were integers to get the distance between them in units
    // of least precision (ULPs).
    static_assert(sizeof(float) == sizeof(int32_t), "What size is float?");

    // memcpy to get around the strict aliasing rule.
    // The compiler knows what we're doing and will just transfer the float
    // values into integer registers.
    int32_t ia, ib;
    memcpy(&ia, &a, sizeof(float));
    memcpy(&ib, &b, sizeof(float));

    // If the signs of the two values aren't the same,
    // return the maximum distance between the two.
    // This is done to avoid integer overflow, and because the bit layout of
    // floats is closer to sign-magnitude than it is to two's complement.
    // This *also* means that if you're checking if a value is close to zero,
    // you should probably just use a fixed epsilon instead of this function.
    if ((ia < 0) != (ib < 0)) return max;

    // If we've satisfied all our caveats above, just subtract the values.
    // The result is the distance between the values in ULPs.
    int32_t distance = ia - ib;
    if (distance < 0) distance = -distance;
    return distance;
}

[phillipjohnston]

This got bumped up - became a problem with M1, since tests are failing due to epsilon comparisons.

[phillipjohnston]

M1 problem was related to the use of long and long long instead of specific fixed-width types, but support for ULP-based comparisons is included. ~It was, however, not necessary, but I am glad I went through the exercise.~ ULP-comparisons revealed what appears to be bug in the underlying gdtoa library that isn't caught with an epsilon-based comparison.

[phillipjohnston]

Although, we may want to play with it in the gdtoa files - that's where most of the comparisons happen.