include/marnav/math/floatingpoint_ulps.hpp

   1 #ifndef MARNAV_MATH_FLOATINGPONIT_ULPS_HPP
   2 #define MARNAV_MATH_FLOATINGPONIT_ULPS_HPP
   3
   4 // implementation based on https://bitbashing.io/comparing-floats.html
   5 // by matt <at> bitbashing.io
   6 //
   7 // this file is distributed under following license: CC BY-SA 4.0
   8 // https://creativecommons.org/licenses/by-sa/4.0/
   9
  10 #include <limits>
  11 #include <cstdint>
  12 #include <cmath>
  13 #include <cstring>
  14
  15 namespace marnav
  16 {
  17 namespace math
  18 {
  19 namespace detail
  20 {
  21 // clang-format off
  22 template <class TF> struct corresponding_integer         {};
  23 template <>         struct corresponding_integer<float>  { using type = int32_t; };
  24 template <>         struct corresponding_integer<double> { using type = int64_t; };
  25
  26 template <class TF>
  27 using corresponding_integer_t = typename corresponding_integer<TF>::type;
  28
  29 template <class TF> struct default_ulps         {};
  30 template <>         struct default_ulps<float>  { static constexpr corresponding_integer_t<float>  value = 1e4; };
  31 template <>         struct default_ulps<double> { static constexpr corresponding_integer_t<double> value = 1e6; };
  32 // clang-format on
  33
  34 template <class TF, class TI = corresponding_integer_t<TF>>
  35 TI ulps_distance(const TF a, const TF b)
  36 {
  37         // save work if the floats are equal, also handles +0 == -0
  38         if (a == b)
  39                 return 0;
  40
  41         // max distance for NaN
  42         if (std::isnan(a) || std::isnan(b))
  43                 return std::numeric_limits<TI>::max();
  44
  45         // max distance for INF
  46         if (std::isinf(a) || std::isinf(b))
  47                 return std::numeric_limits<TI>::max();
  48
  49         TI ia;
  50         TI ib;
  51         ::memcpy(&ia, &a, sizeof(TF));
  52         ::memcpy(&ib, &b, sizeof(TF));
  53
  54         // don't compare differently-signed floats.
  55         if ((ia < 0) != (ib < 0))
  56                 return std::numeric_limits<TI>::max();
  57
  58         // return the absolute value of the distance in ULPs.
  59         const TI distance = ia - ib;
  60         return (distance < 0) ? -distance : distance;
  61 }
  62 }
  63
  64 template <class TF, class TI = detail::corresponding_integer_t<TF>>
  65 bool nearly_equal(TF a, TF b, TF fixed_epsilon = std::numeric_limits<TF>::epsilon(),
  66         TI ulps_epsilon = detail::default_ulps<TF>::value)
  67 {
  68         const auto difference = std::abs(a - b);
  69         if (difference <= fixed_epsilon)
  70                 return true;
  71         return detail::ulps_distance(a, b) <= ulps_epsilon;
  72 }
  73 }
  74 }
  75
  76 #endif