numerics.hpp

Go to the documentation of this file.

#pragma once
 
#include <bit>
#include <cmath>
#include <cstdint>
#include <type_traits>
#include "attributes/common.hpp"
 
#ifdef __AVX512F__
#define IFAVX512(x,y) x
#else
#define IFAVX512(x,y) y
#endif
 
#include "types.hpp"
#include "numerics/bf16.hpp"
#include "numerics/fp16.hpp"
 
namespace ein {
 
template <auto N, auto ... candidates>
concept one_of = ((N==candidates) || ... || false);
 
template <auto N, auto ... candidates>
concept not_one_of = (!one_of<N,candidates...>);
 
template <size_t N> struct integer_traits {};
template <> struct integer_traits<8>  { using signed_t = int8_t;  using unsigned_t = uint8_t; };
template <> struct integer_traits<16> { using signed_t = int16_t; using unsigned_t = uint16_t; };
template <> struct integer_traits<32> { using signed_t = int32_t; using unsigned_t = uint32_t; };
template <> struct integer_traits<64> { using signed_t = int64_t; using unsigned_t = uint64_t; };
 
template <typename T>
requires one_of<sizeof(T),1,2,4,8>
using int_t = typename integer_traits<sizeof(T)*8>::signed_t;
 
template <typename T>
requires one_of<sizeof(T),1,2,4,8>
using uint_t = typename integer_traits<sizeof(T)*8>::unsigned_t;
 
template <size_t N>
struct imm_t {
  static constexpr size_t value = N;
  ein_nodiscard ein_inline ein_const ein_artificial
  consteval operator size_t () const noexcept { return N; }
};
 
template <size_t N>
constinit imm_t<N> imm {};
 
template <typename T>
ein_nodiscard ein_inline ein_pure
constexpr bool cmp_unord(T a, T b) noexcept {
  return isnan(a) || isnan(b);
}
 
extern template bool cmp_unord(float,float) noexcept;
extern template bool cmp_unord(double,double) noexcept;
 
template <typename T>
ein_nodiscard ein_inline ein_pure
constexpr bool cmp_ord(T a, T b) noexcept {
  return !isnan(a) && !isnan(b);
}
 
extern template bool cmp_ord(float,float) noexcept;
extern template bool cmp_ord(double,double) noexcept;
 
template <one_of_t<float,double> T>
ein_nodiscard ein_inline ein_pure
constexpr T scalef(T x, T y) noexcept {
  if consteval {
    // Constexpr path using bit manipulation
 
    using uint_type = uint_t<T>; //  std::conditional_t<std::is_same_v<T, float>, uint32_t, uint64_t>;
    constexpr int exponent_bias = std::is_same_v<T, float> ? 127 : 1023;
    constexpr uint_type exponent_mask = std::is_same_v<T, float> ? uint_type(0x7F800000) : uint_type(0x7FF0000000000000ull);
    constexpr uint_type mantissa_mask = std::is_same_v<T, float> ? uint_type(0x007FFFFF) : uint_type(0x000FFFFFFFFFFFFFull);
 
    if (x == 0.0 || isnan(x) || std::isinf(x))
      return x; // Handle special cases
 
    // Get the raw bits of the floating-point number
    uint_type x_bits = std::bit_cast<uint_type>(x);
 
    // Extract exponent and mantissa
    int exponent = static_cast<int>(((x_bits & exponent_mask) >> (std::is_same_v<T, float> ? 23 : 52)) - exponent_bias);
    uint_type mantissa = x_bits & mantissa_mask;
 
    // Normalize the mantissa (implicit leading 1 for normalized numbers)
    if (exponent != -exponent_bias)
      mantissa |= std::is_same_v<T, float> ? (1u << 23) : (1ull << 52); // Set implicit leading bit
 
    // Scale the exponent by adding y
    exponent += static_cast<int>(y);
 
    // Check for overflow and underflow
    if (exponent > std::numeric_limits<T>::max_exponent)
      return std::numeric_limits<T>::infinity();
 
    if (exponent < std::numeric_limits<T>::min_exponent)
      return static_cast<T>(0.0); // Underflow to zero
 
    // Rebuild the floating-point number from the new exponent and mantissa
    x_bits =
      (x_bits & (std::is_same_v<T, float> ? 0x80000000u : 0x8000000000000000ull)) | // Preserve the sign bit
      (static_cast<uint_type>(exponent + exponent_bias) << (std::is_same_v<T, float> ? 23 : 52)) | // Apply new exponent
      (mantissa & mantissa_mask); // Mask mantissa bits
 
    return std::bit_cast<T>(x_bits);
  } else {
    // Runtime path using fast standard library call
    return std::scalbn(x, static_cast<int>(y));
  }
}
 
enum class ein_nodiscard CMPINT : size_t {
  EQ    = 0x0uz 
, LT    = 0x1uz 
, LE    = 0x2uz 
, FALSE = 0x3uz 
, NE    = 0x4uz 
, NLT   = 0x5uz 
, NLE   = 0x6uz 
, TRUE  = 0x7uz 
};
 
template <CMPINT imm8, typename T>
requires (one_of_t<T,uint8_t,int8_t,uint16_t,int16_t,uint32_t,int32_t,uint64_t,int64_t> && (size_t(imm8) < 8uz))
ein_nodiscard ein_inline ein_const
constexpr bool cmpint(T a, T b) noexcept {
  using enum CMPINT;
  if      constexpr (imm8 == TRUE)  return -1;
  else if constexpr (imm8 == FALSE) return 0;
  else if constexpr (imm8 == LT)    return a < b;
  else if constexpr (imm8 == NLT)   return a >= b;
  else if constexpr (imm8 == LE)    return a <= b;
  else if constexpr (imm8 == NLE)   return a > b;
  else if constexpr (imm8 == EQ)    return a == b;
  else if constexpr (imm8 == NE)    return a != b;
  else static_assert(false);
}
 
constexpr size_t max_fp_comparison_predicate = IFAVX512(32,8);
 
enum class ein_nodiscard CMP : size_t {
  EQ_OQ     = 0x00uz  
, LT_OS     = 0x01uz  
, LE_OS     = 0x02uz  
, UNORD_Q   = 0x03uz  
, NEQ_UQ    = 0x04uz  
, NLT_US    = 0x05uz  
, NLE_US    = 0x06uz  
, ORD_Q     = 0x07uz  
, EQ_UQ     = 0x08uz  
, NGE_US    = 0x09uz  
, NGT_US    = 0x0Auz  
, FALSE_OQ  = 0x0Buz  
, NEQ_OQ    = 0x0Cuz  
, GE_OS     = 0x0Duz  
, GT_OS     = 0x0Euz  
, TRUE_UQ   = 0x0Fuz  
, EQ_OS     = 0x10uz  
, LT_OQ     = 0x11uz  
, LE_OQ     = 0x12uz  
, UNORD_S   = 0x13uz  
, NEQ_US    = 0x14uz  
, NLT_UQ    = 0x15uz  
, NLE_UQ    = 0x16uz  
, ORD_S     = 0x17uz  
, EQ_US     = 0x18uz  
, NGE_UQ    = 0x19uz  
, NGT_UQ    = 0x1Auz  
, FALSE_OS  = 0x1Buz  
, NEQ_OS    = 0x1Cuz  
, GE_OQ     = 0x1Duz  
, GT_OQ     = 0x1Euz  
, TRUE_US   = 0x1Fuz  
};
 
template <CMP imm8, typename T>
requires (one_of_t<T,float,double> && (size_t(imm8) < 32uz))
ein_nodiscard ein_inline ein_pure
constexpr bool cmp(T a, T b) noexcept {
  using enum CMP;
  if      constexpr (imm8 == EQ_OQ)    return cmp_ord(a, b) && (a == b);
  else if constexpr (imm8 == LT_OS)    return cmp_ord(a, b) && (a < b);
  else if constexpr (imm8 == LE_OS)    return cmp_ord(a, b) && (a <= b);
  else if constexpr (imm8 == UNORD_Q)  return cmp_unord(a, b);
  else if constexpr (imm8 == NEQ_UQ)   return cmp_unord(a, b) || (a != b);
  else if constexpr (imm8 == NLT_US)   return cmp_unord(a, b) || !(a < b);
  else if constexpr (imm8 == NLE_US)   return cmp_unord(a, b) || !(a <= b);
  else if constexpr (imm8 == ORD_Q)    return cmp_ord(a, b);
  else if constexpr (imm8 == EQ_UQ)    return cmp_unord(a, b) || (a == b);
  else if constexpr (imm8 == NGE_US)   return cmp_unord(a, b) || !(a >= b);
  else if constexpr (imm8 == NGT_US)   return cmp_unord(a, b) || !(a > b);
  else if constexpr (imm8 == FALSE_OQ) return 0;
  else if constexpr (imm8 == NEQ_OQ)   return cmp_ord(a, b) && (a != b);
  else if constexpr (imm8 == GE_OS)    return cmp_ord(a, b) && (a >= b);
  else if constexpr (imm8 == GT_OS)    return cmp_ord(a, b) && (a > b);
  else if constexpr (imm8 == TRUE_UQ)  return -1;
  else if constexpr (imm8 == EQ_OS)    return cmp_ord(a, b) && (a == b);
  else if constexpr (imm8 == LT_OQ)    return cmp_ord(a, b) && (a < b);
  else if constexpr (imm8 == LE_OQ)    return cmp_ord(a, b) && (a <= b);
  else if constexpr (imm8 == UNORD_S)  return cmp_unord(a, b);
  else if constexpr (imm8 == NEQ_US)   return cmp_unord(a, b) || (a != b);
  else if constexpr (imm8 == NLT_UQ)   return cmp_unord(a, b) || !(a < b);
  else if constexpr (imm8 == NLE_UQ)   return cmp_unord(a, b) || !(a <= b);
  else if constexpr (imm8 == ORD_S)    return cmp_ord(a, b);
  else if constexpr (imm8 == EQ_US)    return cmp_unord(a, b) || (a == b);
  else if constexpr (imm8 == NGE_UQ)   return cmp_unord(a, b) || !(a >= b);
  else if constexpr (imm8 == NGT_UQ)   return cmp_unord(a, b) || !(a > b);
  else if constexpr (imm8 == FALSE_OS) return 0;
  else if constexpr (imm8 == NEQ_OS)   return cmp_ord(a, b) && (a != b);
  else if constexpr (imm8 == GE_OQ)    return cmp_ord(a, b) && (a >= b);
  else if constexpr (imm8 == GT_OQ)    return cmp_ord(a, b) && (a > b);
  else if constexpr (imm8 == TRUE_US)  return -1;
  else static_assert(false);
}
 
#define X extern
#include "numerics.x"
#undef X
 
}
 
#if defined(EIN_TESTING) || defined(EIN_TESTING_NUMERICS)
TEST_CASE("numerics","[numerics]") {
  using namespace ein;
  using Catch::Approx;
 
  SECTION("Concepts: one_of and not_one_of") {
    CHECK(one_of<1, 1, 2, 3>);
    CHECK_FALSE(one_of<4, 1, 2, 3>);
    CHECK(not_one_of<4, 1, 2, 3>);
    CHECK_FALSE(not_one_of<1, 1, 2, 3>);
  }
 
  SECTION("integer_traits") {
    CHECK(std::is_same_v<integer_traits<8>::signed_t, int8_t>);
    CHECK(std::is_same_v<integer_traits<16>::unsigned_t, uint16_t>);
    CHECK(std::is_same_v<integer_traits<32>::signed_t, int32_t>);
    CHECK(std::is_same_v<integer_traits<64>::unsigned_t, uint64_t>);
  }
 
  SECTION("int_t and uint_t") {
    CHECK(std::is_same_v<int_t<int8_t>, int8_t>);
    CHECK(std::is_same_v<uint_t<int16_t>, uint16_t>);
    CHECK(std::is_same_v<int_t<int32_t>, int32_t>);
    CHECK(std::is_same_v<uint_t<int64_t>, uint64_t>);
  }
 
  SECTION("imm_t compile-time constant") {
    constexpr auto imm4 = imm<4>;
    CHECK(imm4.value == 4);
    CHECK(static_cast<size_t>(imm4) == 4);
  }
 
  SECTION("cmp_unord and cmp_ord") {
    CHECK(cmp_unord(NAN, 1.0f));
    CHECK(cmp_unord(1.0f, NAN));
    CHECK_FALSE(cmp_unord(1.0f, 2.0f));
 
    CHECK(cmp_ord(1.0f, 2.0f));
    CHECK_FALSE(cmp_ord(NAN, 1.0f));
    CHECK_FALSE(cmp_ord(1.0f, NAN));
  }
 
  SECTION("scalef") {
    constexpr float x = 2.0f;
    constexpr float y = 3.0f;
    CHECK(scalef(x, y) == Approx(16.0f));
 
    constexpr double a = 4.0;
    constexpr double b = -2.0;
    CHECK(scalef(a, b) == Approx(1.0));
  }
 
  SECTION("CMPINT comparison") {
    CHECK(cmpint<CMPINT::EQ>(5, 5));
    CHECK(cmpint<CMPINT::NE>(5, 4));
    CHECK_FALSE(cmpint<CMPINT::LT>(5, 4));
    CHECK(cmpint<CMPINT::LE>(5, 5));
    CHECK(cmpint<CMPINT::NLT>(5, 5));
  }
 
#ifdef __AVX512F__
  SECTION("CMP floating-point comparison with AVX512") {
    CHECK(cmp<CMP::EQ_OQ>(1.0f, 1.0f));
    CHECK(cmp<CMP::LT_OS>(1.0f, 2.0f));
    CHECK_FALSE(cmp<CMP::GT_OS>(1.0f, 2.0f));
    //CHECK(cmp<CMP::NEQ_UQ>(1.0f, NAN));
  }
#endif // __AVX512F__
}
#endif