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_NEARESTINTEGEROPERATIONS_H
  10 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_NEARESTINTEGEROPERATIONS_H
  11
  12 #include "FEnvImpl.h"
  13 #include "FPBits.h"
  14 #include "rounding_mode.h"
  15
  16 #include "hdr/math_macros.h"
  17 #include "src/__support/CPP/type_traits.h"
  18 #include "src/__support/common.h"
  19 #include "src/__support/macros/config.h"
  20
  21 namespace LIBC_NAMESPACE_DECL {
  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   using StorageType = typename FPBits<T>::StorageType;
  27   FPBits<T> bits(x);
  28
  29   // If x is infinity or NaN, return it.
  30   // If it is zero also we should return it as is, but the logic
  31   // later in this function takes care of it. But not doing a zero
  32   // check, we improve the run time of non-zero values.
  33   if (bits.is_inf_or_nan())
  34     return x;
  35
  36   int exponent = bits.get_exponent();
  37
  38   // If the exponent is greater than the most negative mantissa
  39   // exponent, then x is already an integer.
  40   if (exponent >= static_cast<int>(FPBits<T>::FRACTION_LEN))
  41     return x;
  42
  43   // If the exponent is such that abs(x) is less than 1, then return 0.
  44   if (exponent <= -1)
  45     return FPBits<T>::zero(bits.sign()).get_val();
  46
  47   int trim_size = FPBits<T>::FRACTION_LEN - exponent;
  48   StorageType trunc_mantissa =
  49       static_cast<StorageType>((bits.get_mantissa() >> trim_size) << trim_size);
  50   bits.set_mantissa(trunc_mantissa);
  51   return bits.get_val();
  52 }
  53
  54 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
  55 LIBC_INLINE T ceil(T x) {
  56   using StorageType = typename FPBits<T>::StorageType;
  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.is_neg();
  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>(FPBits<T>::FRACTION_LEN))
  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 = FPBits<T>::FRACTION_LEN - exponent;
  79   StorageType x_u = bits.uintval();
  80   StorageType trunc_u =
  81       static_cast<StorageType>((x_u >> trim_size) << trim_size);
  82
  83   // If x is already an integer, return it.
  84   if (trunc_u == x_u)
  85     return x;
  86
  87   bits.set_uintval(trunc_u);
  88   T trunc_value = bits.get_val();
  89
  90   // If x is negative, the ceil operation is equivalent to the trunc operation.
  91   if (is_neg)
  92     return trunc_value;
  93
  94   return trunc_value + T(1.0);
  95 }
  96
  97 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
  98 LIBC_INLINE T floor(T x) {
  99   FPBits<T> bits(x);
 100   if (bits.is_neg()) {
 101     return -ceil(-x);
 102   } else {
 103     return trunc(x);
 104   }
 105 }
 106
 107 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
 108 LIBC_INLINE T round(T x) {
 109   using StorageType = typename FPBits<T>::StorageType;
 110   FPBits<T> bits(x);
 111
 112   // If x is infinity NaN or zero, return it.
 113   if (bits.is_inf_or_nan() || bits.is_zero())
 114     return x;
 115
 116   int exponent = bits.get_exponent();
 117
 118   // If the exponent is greater than the most negative mantissa
 119   // exponent, then x is already an integer.
 120   if (exponent >= static_cast<int>(FPBits<T>::FRACTION_LEN))
 121     return x;
 122
 123   if (exponent == -1) {
 124     // Absolute value of x is greater than equal to 0.5 but less than 1.
 125     return FPBits<T>::one(bits.sign()).get_val();
 126   }
 127
 128   if (exponent <= -2) {
 129     // Absolute value of x is less than 0.5.
 130     return FPBits<T>::zero(bits.sign()).get_val();
 131   }
 132
 133   uint32_t trim_size = FPBits<T>::FRACTION_LEN - exponent;
 134   bool half_bit_set =
 135       bool(bits.get_mantissa() & (StorageType(1) << (trim_size - 1)));
 136   StorageType x_u = bits.uintval();
 137   StorageType trunc_u =
 138       static_cast<StorageType>((x_u >> trim_size) << trim_size);
 139
 140   // If x is already an integer, return it.
 141   if (trunc_u == x_u)
 142     return x;
 143
 144   bits.set_uintval(trunc_u);
 145   T trunc_value = bits.get_val();
 146
 147   if (!half_bit_set) {
 148     // Franctional part is less than 0.5 so round value is the
 149     // same as the trunc value.
 150     return trunc_value;
 151   } else {
 152     return bits.is_neg() ? trunc_value - T(1.0) : trunc_value + T(1.0);
 153   }
 154 }
 155
 156 template <typename T>
 157 LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<T>, T>
 158 round_using_specific_rounding_mode(T x, int rnd) {
 159   using StorageType = typename FPBits<T>::StorageType;
 160   FPBits<T> bits(x);
 161
 162   // If x is infinity NaN or zero, return it.
 163   if (bits.is_inf_or_nan() || bits.is_zero())
 164     return x;
 165
 166   bool is_neg = bits.is_neg();
 167   int exponent = bits.get_exponent();
 168
 169   // If the exponent is greater than the most negative mantissa
 170   // exponent, then x is already an integer.
 171   if (exponent >= static_cast<int>(FPBits<T>::FRACTION_LEN))
 172     return x;
 173
 174   if (exponent <= -1) {
 175     switch (rnd) {
 176     case FP_INT_DOWNWARD:
 177       return is_neg ? T(-1.0) : T(0.0);
 178     case FP_INT_UPWARD:
 179       return is_neg ? T(-0.0) : T(1.0);
 180     case FP_INT_TOWARDZERO:
 181       return is_neg ? T(-0.0) : T(0.0);
 182     case FP_INT_TONEARESTFROMZERO:
 183       if (exponent < -1)
 184         return is_neg ? T(-0.0) : T(0.0); // abs(x) < 0.5
 185       return is_neg ? T(-1.0) : T(1.0);   // abs(x) >= 0.5
 186     case FP_INT_TONEAREST:
 187     default:
 188       if (exponent <= -2 || bits.get_mantissa() == 0)
 189         return is_neg ? T(-0.0) : T(0.0); // abs(x) <= 0.5
 190       else
 191         return is_neg ? T(-1.0) : T(1.0); // abs(x) > 0.5
 192     }
 193   }
 194
 195   uint32_t trim_size = FPBits<T>::FRACTION_LEN - exponent;
 196   StorageType x_u = bits.uintval();
 197   StorageType trunc_u =
 198       static_cast<StorageType>((x_u >> trim_size) << trim_size);
 199
 200   // If x is already an integer, return it.
 201   if (trunc_u == x_u)
 202     return x;
 203
 204   FPBits<T> new_bits(trunc_u);
 205   T trunc_value = new_bits.get_val();
 206
 207   StorageType trim_value =
 208       bits.get_mantissa() &
 209       static_cast<StorageType>(((StorageType(1) << trim_size) - 1));
 210   StorageType half_value =
 211       static_cast<StorageType>((StorageType(1) << (trim_size - 1)));
 212   // If exponent is 0, trimSize will be equal to the mantissa width, and
 213   // truncIsOdd` will not be correct. So, we handle it as a special case
 214   // below.
 215   StorageType trunc_is_odd =
 216       new_bits.get_mantissa() & (StorageType(1) << trim_size);
 217
 218   switch (rnd) {
 219   case FP_INT_DOWNWARD:
 220     return is_neg ? trunc_value - T(1.0) : trunc_value;
 221   case FP_INT_UPWARD:
 222     return is_neg ? trunc_value : trunc_value + T(1.0);
 223   case FP_INT_TOWARDZERO:
 224     return trunc_value;
 225   case FP_INT_TONEARESTFROMZERO:
 226     if (trim_value >= half_value)
 227       return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
 228     return trunc_value;
 229   case FP_INT_TONEAREST:
 230   default:
 231     if (trim_value > half_value) {
 232       return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
 233     } else if (trim_value == half_value) {
 234       if (exponent == 0)
 235         return is_neg ? T(-2.0) : T(2.0);
 236       if (trunc_is_odd)
 237         return is_neg ? trunc_value - T(1.0) : trunc_value + T(1.0);
 238       else
 239         return trunc_value;
 240     } else {
 241       return trunc_value;
 242     }
 243   }
 244 }
 245
 246 template <typename T>
 247 LIBC_INLINE cpp::enable_if_t<cpp::is_floating_point_v<T>, T>
 248 round_using_current_rounding_mode(T x) {
 249   int rounding_mode = quick_get_round();
 250
 251   switch (rounding_mode) {
 252   case FE_DOWNWARD:
 253     return round_using_specific_rounding_mode(x, FP_INT_DOWNWARD);
 254   case FE_UPWARD:
 255     return round_using_specific_rounding_mode(x, FP_INT_UPWARD);
 256   case FE_TOWARDZERO:
 257     return round_using_specific_rounding_mode(x, FP_INT_TOWARDZERO);
 258   case FE_TONEAREST:
 259     return round_using_specific_rounding_mode(x, FP_INT_TONEAREST);
 260   default:
 261     __builtin_unreachable();
 262   }
 263 }
 264
 265 template <bool IsSigned, typename T>
 266 LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<T>, T>
 267 fromfp(T x, int rnd, unsigned int width) {
 268   using StorageType = typename FPBits<T>::StorageType;
 269
 270   constexpr StorageType EXPLICIT_BIT =
 271       FPBits<T>::SIG_MASK - FPBits<T>::FRACTION_MASK;
 272
 273   if (width == 0U) {
 274     raise_except_if_required(FE_INVALID);
 275     return FPBits<T>::quiet_nan().get_val();
 276   }
 277
 278   FPBits<T> bits(x);
 279
 280   if (bits.is_inf_or_nan()) {
 281     raise_except_if_required(FE_INVALID);
 282     return FPBits<T>::quiet_nan().get_val();
 283   }
 284
 285   T rounded_value = round_using_specific_rounding_mode(x, rnd);
 286
 287   if constexpr (IsSigned) {
 288     // T can't hold a finite number >= 2.0 * 2^EXP_BIAS.
 289     if (width - 1 > FPBits<T>::EXP_BIAS)
 290       return rounded_value;
 291
 292     StorageType range_exp =
 293         static_cast<StorageType>(width - 1 + FPBits<T>::EXP_BIAS);
 294     // rounded_value < -2^(width - 1)
 295     T range_min =
 296         FPBits<T>::create_value(Sign::NEG, range_exp, EXPLICIT_BIT).get_val();
 297     if (rounded_value < range_min) {
 298       raise_except_if_required(FE_INVALID);
 299       return FPBits<T>::quiet_nan().get_val();
 300     }
 301     // rounded_value > 2^(width - 1) - 1
 302     T range_max =
 303         FPBits<T>::create_value(Sign::POS, range_exp, EXPLICIT_BIT).get_val() -
 304         T(1.0);
 305     if (rounded_value > range_max) {
 306       raise_except_if_required(FE_INVALID);
 307       return FPBits<T>::quiet_nan().get_val();
 308     }
 309
 310     return rounded_value;
 311   }
 312
 313   if (rounded_value < T(0.0)) {
 314     raise_except_if_required(FE_INVALID);
 315     return FPBits<T>::quiet_nan().get_val();
 316   }
 317
 318   // T can't hold a finite number >= 2.0 * 2^EXP_BIAS.
 319   if (width > FPBits<T>::EXP_BIAS)
 320     return rounded_value;
 321
 322   StorageType range_exp = static_cast<StorageType>(width + FPBits<T>::EXP_BIAS);
 323   // rounded_value > 2^width - 1
 324   T range_max =
 325       FPBits<T>::create_value(Sign::POS, range_exp, EXPLICIT_BIT).get_val() -
 326       T(1.0);
 327   if (rounded_value > range_max) {
 328     raise_except_if_required(FE_INVALID);
 329     return FPBits<T>::quiet_nan().get_val();
 330   }
 331
 332   return rounded_value;
 333 }
 334
 335 template <bool IsSigned, typename T>
 336 LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<T>, T>
 337 fromfpx(T x, int rnd, unsigned int width) {
 338   T rounded_value = fromfp<IsSigned>(x, rnd, width);
 339   FPBits<T> bits(rounded_value);
 340
 341   if (!bits.is_nan() && rounded_value != x)
 342     raise_except_if_required(FE_INEXACT);
 343
 344   return rounded_value;
 345 }
 346
 347 namespace internal {
 348
 349 template <typename FloatType, typename IntType,
 350           cpp::enable_if_t<cpp::is_floating_point_v<FloatType> &&
 351                                cpp::is_integral_v<IntType>,
 352                            int> = 0>
 353 LIBC_INLINE IntType rounded_float_to_signed_integer(FloatType x) {
 354   constexpr IntType INTEGER_MIN = (IntType(1) << (sizeof(IntType) * 8 - 1));
 355   constexpr IntType INTEGER_MAX = -(INTEGER_MIN + 1);
 356   FPBits<FloatType> bits(x);
 357   auto set_domain_error_and_raise_invalid = []() {
 358     set_errno_if_required(EDOM);
 359     raise_except_if_required(FE_INVALID);
 360   };
 361
 362   if (bits.is_inf_or_nan()) {
 363     set_domain_error_and_raise_invalid();
 364     return bits.is_neg() ? INTEGER_MIN : INTEGER_MAX;
 365   }
 366
 367   int exponent = bits.get_exponent();
 368   constexpr int EXPONENT_LIMIT = sizeof(IntType) * 8 - 1;
 369   if (exponent > EXPONENT_LIMIT) {
 370     set_domain_error_and_raise_invalid();
 371     return bits.is_neg() ? INTEGER_MIN : INTEGER_MAX;
 372   } else if (exponent == EXPONENT_LIMIT) {
 373     if (bits.is_pos() || bits.get_mantissa() != 0) {
 374       set_domain_error_and_raise_invalid();
 375       return bits.is_neg() ? INTEGER_MIN : INTEGER_MAX;
 376     }
 377     // If the control reaches here, then it means that the rounded
 378     // value is the most negative number for the signed integer type IntType.
 379   }
 380
 381   // For all other cases, if `x` can fit in the integer type `IntType`,
 382   // we just return `x`. static_cast will convert the floating
 383   // point value to the exact integer value.
 384   return static_cast<IntType>(x);
 385 }
 386
 387 } // namespace internal
 388
 389 template <typename FloatType, typename IntType,
 390           cpp::enable_if_t<cpp::is_floating_point_v<FloatType> &&
 391                                cpp::is_integral_v<IntType>,
 392                            int> = 0>
 393 LIBC_INLINE IntType round_to_signed_integer(FloatType x) {
 394   return internal::rounded_float_to_signed_integer<FloatType, IntType>(
 395       round(x));
 396 }
 397
 398 template <typename FloatType, typename IntType,
 399           cpp::enable_if_t<cpp::is_floating_point_v<FloatType> &&
 400                                cpp::is_integral_v<IntType>,
 401                            int> = 0>
 402 LIBC_INLINE IntType
 403 round_to_signed_integer_using_current_rounding_mode(FloatType x) {
 404   return internal::rounded_float_to_signed_integer<FloatType, IntType>(
 405       round_using_current_rounding_mode(x));
 406 }
 407
 408 } // namespace fputil
 409 } // namespace LIBC_NAMESPACE_DECL
 410
 411 #endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_NEARESTINTEGEROPERATIONS_H