Exploit IEEE 754 spec's property on comparing non-zero numbers.

float_cmp.cc

#include <algorithm>
#include <cstdint>
#include <cstdio>
#include <limits>

using namespace std;

template <size_t Bits>
struct IntegerTrait {};

template <>
struct IntegerTrait<32> { using UnsignedType = uint32_t; };

template <>
struct IntegerTrait<64> { using UnsignedType = uint64_t; };

template <typename FloatType>
using CorrespondingUnsignedInteger =
    typename IntegerTrait<sizeof(FloatType) * 8>::UnsignedType;

template <typename FloatType>
CorrespondingUnsignedInteger<FloatType> ULPCast(FloatType x) {
  union {
    FloatType f;
    CorrespondingUnsignedInteger<FloatType> i;
  } x_ = { .f = x };
  return x_.i;
}

// cf. https://randomascii.wordpress.com/2012/02/25/comparing-floating-point-numbers-2012-edition/
template <typename FloatType>
uint64_t CountULPs(FloatType a, FloatType b) {
  if (a < b) { swap(a, b); }
  return ULPCast<FloatType>(a) - ULPCast<FloatType>(b);
}

template <typename FloatType>
bool IsNearEnough(FloatType a, FloatType b, uint64_t max_diff) {
  if (a == 0.0 && b == 0.0) { return true; }
  if (a * b < 0.0) { return false; }
  return CountULPs(a, b) <= max_diff;
}

int main() {
  double a = 1.0;
  double b = 1.0 + numeric_limits<double>::epsilon();
  printf("ULPs: %llu\n", CountULPs(a, b));
  printf("(Almost) equal?: %s\n", IsNearEnough(a, b, 10) ? "yes" : "no");
  return 0;
}