libc/src/__support/FPUtil/NearestIntegerOperations.h

   1 //===-- Nearest integer floating-point operations ---------------*- C++ -*-===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8
   9 #ifndef LLVM_LIBC_SRC_SUPPORT_FPUTIL_NEAREST_INTEGER_OPERATIONS_H
  10 #define LLVM_LIBC_SRC_SUPPORT_FPUTIL_NEAREST_INTEGER_OPERATIONS_H
  11
  12 #include "FEnvImpl.h"
  13 #include "FPBits.h"
  14 #include "rounding_mode.h"
  15
  16 #include "src/__support/CPP/type_traits.h"
  17 #include "src/__support/common.h"
  18
  19 #include <math.h>
  20
  21 namespace __llvm_libc {
  22 namespace fputil {
  23
  24 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
  25 LIBC_INLINE T trunc(T x) {
  26   FPBits<T> bits(x);
  27
  28   // If x is infinity or NaN, return it.
  29   // If it is zero also we should return it as is, but the logic
  30   // later in this function takes care of it. But not doing a zero
  31   // check, we improve the run time of non-zero values.
  32   if (bits.is_inf_or_nan())
  33     return x;
  34
  35   int exponent = bits.get_exponent();
  36
  37   // If the exponent is greater than the most negative mantissa
  38   // exponent, then x is already an integer.
  39   if (exponent >= static_cast<int>(MantissaWidth<T>::VALUE))
  40     return x;
  41
  42   // If the exponent is such that abs(x) is less than 1, then return 0.
  43   if (exponent <= -1) {
  44     if (bits.get_sign())
  45       return T(-0.0);
  46     else
  47       return T(0.0);
  48   }
  49
  50   int trim_size = MantissaWidth<T>::VALUE - exponent;
  51   bits.set_mantissa((bits.get_mantissa() >> trim_size) << trim_size);
  52   return T(bits);
  53 }
  54
  55 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
  56 LIBC_INLINE T ceil(T x) {
  57   FPBits<T> bits(x);
  58
  59   // If x is infinity NaN or zero, return it.
  60   if (bits.is_inf_or_nan() || bits.is_zero())
  61     return x;
  62
  63   bool is_neg = bits.get_sign();
  64   int exponent = bits.get_exponent();
  65
  66   // If the exponent is greater than the most negative mantissa
  67   // exponent, then x is already an integer.
  68   if (exponent >= static_cast<int>(MantissaWidth<T>::VALUE))
  69     return x;
  70
  71   if (exponent <= -1) {
  72     if (is_neg)
  73       return T(-0.0);
  74     else
  75       return T(1.0);
  76   }
  77
  78   uint32_t trim_size = MantissaWidth<T>::VALUE - exponent;
  79   bits.set_mantissa((bits.get_mantissa() >> trim_size) << trim_size);
  80   T trunc_value = T(bits);
  81
  82   // If x is already an integer, return it.
  83   if (trunc_value == x)
  84     return x;
  85
  86   // If x is negative, the ceil operation is equivalent to the trunc operation.
  87   if (is_neg)
  88     return trunc_value;
  89
  90   return trunc_value + T(1.0);
  91 }
  92
  93 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
  94 LIBC_INLINE T floor(T x) {
  95   FPBits<T> bits(x);
  96   if (bits.get_sign()) {
  97     return -ceil(-x);
  98   } else {
  99     return trunc(x);
 100   }
 101 }
 102
 103 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
 104 LIBC_INLINE T round(T x) {
 105   using UIntType = typename FPBits<T>::UIntType;
 106   FPBits<T> bits(x);
 107
 108   // If x is infinity NaN or zero, return it.
 109   if (bits.is_inf_or_nan() || bits.is_zero())
 110     return x;
 111
 112   bool is_neg = bits.get_sign();
 113   int exponent = bits.get_exponent();
 114
 115   // If the exponent is greater than the most negative mantissa
 116   // exponent, then x is already an integer.
 117   if (exponent >= static_cast<int>(MantissaWidth<T>::VALUE))
 118     return x;
 119
 120   if (exponent == -1) {
 121     // Absolute value of x is greater than equal to 0.5 but less than 1.
 122     if (is_neg)
 123       return T(-1.0);
 124     else
 125       return T(1.0);
 126   }
 127
 128   if (exponent <= -2) {
 129     // Absolute value of x is less than 0.5.
 130     if (is_neg)
 131       return T(-0.0);
 132     else
 133       return T(0.0);
 134   }
 135
 136   uint32_t trim_size = MantissaWidth<T>::VALUE - exponent;
 137   bool half_bit_set =
 138       bool(bits.get_mantissa() & (UIntType(1) << (trim_size - 1)));
 139   bits.set_mantissa((bits.get_mantissa() >> trim_size) << trim_size);
 140   T trunc_value = T(bits);
 141
 142   // If x is already an integer, return it.
 143   if (trunc_value == x)
 144     return x;
 145
 146   if (!half_bit_set) {
 147     // Franctional part is less than 0.5 so round value is the
 148     // same as the trunc value.
 149     return trunc_value;
 150   } else {
 151     return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
 152   }
 153 }
 154
 155 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
 156 LIBC_INLINE T round_using_current_rounding_mode(T x) {
 157   using UIntType = typename FPBits<T>::UIntType;
 158   FPBits<T> bits(x);
 159
 160   // If x is infinity NaN or zero, return it.
 161   if (bits.is_inf_or_nan() || bits.is_zero())
 162     return x;
 163
 164   bool is_neg = bits.get_sign();
 165   int exponent = bits.get_exponent();
 166   int rounding_mode = quick_get_round();
 167
 168   // If the exponent is greater than the most negative mantissa
 169   // exponent, then x is already an integer.
 170   if (exponent >= static_cast<int>(MantissaWidth<T>::VALUE))
 171     return x;
 172
 173   if (exponent <= -1) {
 174     switch (rounding_mode) {
 175     case FE_DOWNWARD:
 176       return is_neg ? T(-1.0) : T(0.0);
 177     case FE_UPWARD:
 178       return is_neg ? T(-0.0) : T(1.0);
 179     case FE_TOWARDZERO:
 180       return is_neg ? T(-0.0) : T(0.0);
 181     case FE_TONEAREST:
 182       if (exponent <= -2 || bits.get_mantissa() == 0)
 183         return is_neg ? T(-0.0) : T(0.0); // abs(x) <= 0.5
 184       else
 185         return is_neg ? T(-1.0) : T(1.0); // abs(x) > 0.5
 186     default:
 187       __builtin_unreachable();
 188     }
 189   }
 190
 191   uint32_t trim_size = MantissaWidth<T>::VALUE - exponent;
 192   FPBits<T> new_bits = bits;
 193   new_bits.set_mantissa((bits.get_mantissa() >> trim_size) << trim_size);
 194   T trunc_value = T(new_bits);
 195
 196   // If x is already an integer, return it.
 197   if (trunc_value == x)
 198     return x;
 199
 200   UIntType trim_value = bits.get_mantissa() & ((UIntType(1) << trim_size) - 1);
 201   UIntType half_value = (UIntType(1) << (trim_size - 1));
 202   // If exponent is 0, trimSize will be equal to the mantissa width, and
 203   // truncIsOdd` will not be correct. So, we handle it as a special case
 204   // below.
 205   UIntType trunc_is_odd = new_bits.get_mantissa() & (UIntType(1) << trim_size);
 206
 207   switch (rounding_mode) {
 208   case FE_DOWNWARD:
 209     return is_neg ? trunc_value - T(1.0) : trunc_value;
 210   case FE_UPWARD:
 211     return is_neg ? trunc_value : trunc_value + T(1.0);
 212   case FE_TOWARDZERO:
 213     return trunc_value;
 214   case FE_TONEAREST:
 215     if (trim_value > half_value) {
 216       return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
 217     } else if (trim_value == half_value) {
 218       if (exponent == 0)
 219         return is_neg ? T(-2.0) : T(2.0);
 220       if (trunc_is_odd)
 221         return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
 222       else
 223         return trunc_value;
 224     } else {
 225       return trunc_value;
 226     }
 227   default:
 228     __builtin_unreachable();
 229   }
 230 }
 231
 232 namespace internal {
 233
 234 template <typename F, typename I,
 235           cpp::enable_if_t<cpp::is_floating_point_v<F> && cpp::is_integral_v<I>,
 236                            int> = 0>
 237 LIBC_INLINE I rounded_float_to_signed_integer(F x) {
 238   constexpr I INTEGER_MIN = (I(1) << (sizeof(I) * 8 - 1));
 239   constexpr I INTEGER_MAX = -(INTEGER_MIN + 1);
 240   FPBits<F> bits(x);
 241   auto set_domain_error_and_raise_invalid = []() {
 242     set_errno_if_required(EDOM);
 243     raise_except_if_required(FE_INVALID);
 244   };
 245
 246   if (bits.is_inf_or_nan()) {
 247     set_domain_error_and_raise_invalid();
 248     return bits.get_sign() ? INTEGER_MIN : INTEGER_MAX;
 249   }
 250
 251   int exponent = bits.get_exponent();
 252   constexpr int EXPONENT_LIMIT = sizeof(I) * 8 - 1;
 253   if (exponent > EXPONENT_LIMIT) {
 254     set_domain_error_and_raise_invalid();
 255     return bits.get_sign() ? INTEGER_MIN : INTEGER_MAX;
 256   } else if (exponent == EXPONENT_LIMIT) {
 257     if (bits.get_sign() == 0 || bits.get_mantissa() != 0) {
 258       set_domain_error_and_raise_invalid();
 259       return bits.get_sign() ? INTEGER_MIN : INTEGER_MAX;
 260     }
 261     // If the control reaches here, then it means that the rounded
 262     // value is the most negative number for the signed integer type I.
 263   }
 264
 265   // For all other cases, if `x` can fit in the integer type `I`,
 266   // we just return `x`. static_cast will convert the floating
 267   // point value to the exact integer value.
 268   return static_cast<I>(x);
 269 }
 270
 271 } // namespace internal
 272
 273 template <typename F, typename I,
 274           cpp::enable_if_t<cpp::is_floating_point_v<F> && cpp::is_integral_v<I>,
 275                            int> = 0>
 276 LIBC_INLINE I round_to_signed_integer(F x) {
 277   return internal::rounded_float_to_signed_integer<F, I>(round(x));
 278 }
 279
 280 template <typename F, typename I,
 281           cpp::enable_if_t<cpp::is_floating_point_v<F> && cpp::is_integral_v<I>,
 282                            int> = 0>
 283 LIBC_INLINE I round_to_signed_integer_using_current_rounding_mode(F x) {
 284   return internal::rounded_float_to_signed_integer<F, I>(
 285       round_using_current_rounding_mode(x));
 286 }
 287
 288 } // namespace fputil
 289 } // namespace __llvm_libc
 290
 291 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_NEAREST_INTEGER_OPERATIONS_H