Closed phillipjohnston closed 9 months ago
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;
}
This got bumped up - became a problem with M1, since tests are failing due to epsilon comparisons.
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.
Although, we may want to play with it in the gdtoa files - that's where most of the comparisons happen.
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.