floating_eq.hpp

Go to the documentation of this file.

#pragma once
 
#include <algorithm>
#include <bit>
#include <cmath>
#include <complex>
#include <cstdint>
#include <cstdlib>
#include <limits>
#include <type_traits>
#include <uni20/core/numeric_limits.hpp>
#include <uni20/core/scalar_concepts.hpp>
#include <uni20/core/scalar_traits.hpp>
#if UNI20_HAS_MPREAL
#include "detail/mpreal_ulp.hpp"
#endif
 
namespace uni20::check
{
 
namespace detail
{
 
template <typename T>
inline constexpr bool has_ieee_binary_interchange_layout_v =
    (sizeof(T) == 4 && uni20::numeric_limits<T>::digits() == 24 && uni20::numeric_limits<T>::max_exponent == 128) ||
    (sizeof(T) == 8 && uni20::numeric_limits<T>::digits() == 53 && uni20::numeric_limits<T>::max_exponent == 1024) ||
    (sizeof(T) == 16 && uni20::numeric_limits<T>::digits() == 113 && uni20::numeric_limits<T>::max_exponent == 16384);
 
} // namespace detail
 
template <typename T>
concept IeeeBinaryReal =
    uni20::Real<std::remove_cvref_t<T>> && uni20::numeric_limits<std::remove_cvref_t<T>>::is_iec559 &&
    uni20::numeric_limits<std::remove_cvref_t<T>>::radix == 2 && std::is_trivially_copyable_v<std::remove_cvref_t<T>> &&
    detail::has_ieee_binary_interchange_layout_v<std::remove_cvref_t<T>>;
 
template <typename T>
concept UlpOrderedReal =
    IeeeBinaryReal<T> ||
    (std::same_as<std::remove_cvref_t<T>, long double> && uni20::numeric_limits<long double>::is_iec559 &&
     uni20::numeric_limits<long double>::radix == 2 && uni20::numeric_limits<long double>::digits() == 64 &&
     uni20::numeric_limits<long double>::min_exponent == -16381 &&
     uni20::numeric_limits<long double>::max_exponent == 16384 &&
     uni20::numeric_limits<long double>::denorm_min() > 0 &&
     uni20::numeric_limits<long double>::denorm_min() < uni20::numeric_limits<long double>::min());
 
namespace detail
{
 
template <UlpOrderedReal T>
using ulp_uint_t =
    std::conditional_t<sizeof(T) == 4, std::uint32_t, std::conditional_t<sizeof(T) == 8, std::uint64_t, __uint128_t>>;
 
template <UlpOrderedReal T> constexpr ulp_uint_t<T> ordered_ulp_key(T value)
{
  using UInt = ulp_uint_t<T>;
  if constexpr (IeeeBinaryReal<T>)
  {
    constexpr UInt SignBit = UInt{1} << (sizeof(UInt) * 8 - 1);
    UInt const bits = std::bit_cast<UInt>(value);
 
    // Collapse the two signed-zero encodings and order negative values before
    // positive values while preserving adjacency between finite values.
    return (bits & SignBit) != 0 ? ~bits + UInt{1} : bits | SignBit;
  }
  else
  {
    // Finite fp80 values only: callers handle infinities and NaNs first.
    // frexp/ldexp recover the integer significand exactly, without touching
    // padding or depending on byte order or the explicit integer-bit encoding.
    // Each normal binade contains 2^63 values; the subnormal binade uses the
    // same spacing as the first normal one.
    constexpr UInt SignBit = UInt{1} << 78;
    if (value == T{0}) return SignBit;
    int exponent;
    T const fraction = std::frexp(std::abs(value), &exponent);
    UInt significand = static_cast<UInt>(std::ldexp(fraction, 64));
    constexpr int MinExponent = uni20::numeric_limits<T>::min_exponent;
    UInt magnitude;
    if (exponent < MinExponent)
      magnitude = significand >> (MinExponent - exponent);
    else
      magnitude = (UInt(exponent - MinExponent) << 63) + significand;
    return std::signbit(value) ? SignBit - magnitude : SignBit + magnitude;
  }
}
 
template <UlpOrderedReal T> constexpr ulp_uint_t<T> ulp_distance_magnitude(T a, T b)
{
  auto const ai = ordered_ulp_key(a);
  auto const bi = ordered_ulp_key(b);
  return ai < bi ? bi - ai : ai - bi;
}
 
} // namespace detail
 
template <UlpOrderedReal T> inline long long float_distance(T a, T b)
{
  if (std::isnan(a) || std::isnan(b))
  {
    return std::numeric_limits<long long>::max();
  }
  if (std::isinf(a) || std::isinf(b))
  {
    return (a == b) ? 0 : std::numeric_limits<long long>::max();
  }
  if (a == b)
  {
    return 0; // handles +0 == -0
  }
 
  using UInt = detail::ulp_uint_t<T>;
  auto const ai = detail::ordered_ulp_key(a);
  auto const bi = detail::ordered_ulp_key(b);
  UInt const magnitude = ai < bi ? bi - ai : ai - bi;
  constexpr auto MaxDistance = std::numeric_limits<long long>::max();
  if (magnitude >= static_cast<UInt>(MaxDistance))
  {
    return bi > ai ? MaxDistance : -MaxDistance;
  }
 
  auto const signed_magnitude = static_cast<long long>(magnitude);
  return bi > ai ? signed_magnitude : -signed_magnitude;
}
 
template <typename T> struct FloatingULP;
 
template <typename T>
  requires UlpOrderedReal<T>
struct FloatingULP<T>
{
    static bool eq(T a, T b, std::int64_t max_ulps = 4)
    {
      if (max_ulps < 0 || std::isnan(a) || std::isnan(b))
      {
        return false;
      }
      if (a == b)
      {
        return true;
      }
      if (std::isinf(a) || std::isinf(b))
      {
        return false;
      }
 
      using UInt = detail::ulp_uint_t<T>;
      return detail::ulp_distance_magnitude(a, b) <= static_cast<UInt>(max_ulps);
    }
};
 
#if UNI20_HAS_MPREAL
template <> struct FloatingULP<mpreal>
{
    static bool eq(mpreal const& a, mpreal const& b, std::int64_t max_ulps = 4)
    {
      return detail::compare_mpreal(a, b, max_ulps).equal;
    }
};
 
inline long long float_distance(mpreal const& a, mpreal const& b) { return detail::compare_mpreal(a, b, 0).distance; }
 
inline long long float_abs_distance(mpreal const& a, mpreal const& b) { return std::llabs(float_distance(a, b)); }
#endif
 
template <typename T>
  requires uni20::Complex<T> && UlpOrderedReal<uni20::make_real_t<T>>
struct FloatingULP<T>
{
    static bool eq(T const& a, T const& b, std::int64_t max_ulps = 4)
    {
      using S = uni20::make_real_t<T>;
      return FloatingULP<S>::eq(a.real(), b.real(), max_ulps) && FloatingULP<S>::eq(a.imag(), b.imag(), max_ulps);
    }
};
 
template <typename T>
concept UlpComparable = requires(T a, T b) {
  { FloatingULP<T>::eq(a, b) } -> std::same_as<bool>;
};
 
template <UlpOrderedReal T> inline long long float_abs_distance(T a, T b)
{
  auto dist = float_distance(a, b);
  return (dist == std::numeric_limits<long long>::max()) ? std::numeric_limits<long long>::max() : std::llabs(dist);
}
 
template <uni20::Complex T>
  requires UlpOrderedReal<uni20::make_real_t<T>>
inline long long float_abs_distance(T a, T b)
{
  auto dr = float_abs_distance(a.real(), b.real());
  auto di = float_abs_distance(a.imag(), b.imag());
  return std::max(dr, di);
}
 
} // namespace uni20::check