libc/src/__support/FPUtil/ManipulationFunctions.h

   1 //===-- Floating-point manipulation functions -------------------*- 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_MANIPULATIONFUNCTIONS_H
  10 #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H
  11
  12 #include "FPBits.h"
  13 #include "NearestIntegerOperations.h"
  14 #include "NormalFloat.h"
  15 #include "cast.h"
  16 #include "dyadic_float.h"
  17 #include "rounding_mode.h"
  18
  19 #include "hdr/math_macros.h"
  20 #include "src/__support/CPP/bit.h"
  21 #include "src/__support/CPP/limits.h" // INT_MAX, INT_MIN
  22 #include "src/__support/CPP/type_traits.h"
  23 #include "src/__support/FPUtil/FEnvImpl.h"
  24 #include "src/__support/macros/attributes.h"
  25 #include "src/__support/macros/config.h"
  26 #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
  27
  28 namespace LIBC_NAMESPACE_DECL {
  29 namespace fputil {
  30
  31 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
  32 LIBC_INLINE T frexp(T x, int &exp) {
  33   FPBits<T> bits(x);
  34   if (bits.is_inf_or_nan()) {
  35 #ifdef LIBC_FREXP_INF_NAN_EXPONENT
  36     // The value written back to the second parameter when calling
  37     // frexp/frexpf/frexpl` with `+/-Inf`/`NaN` is unspecified in the standard.
  38     // Set the exp value for Inf/NaN inputs explicitly to
  39     // LIBC_FREXP_INF_NAN_EXPONENT if it is defined.
  40     exp = LIBC_FREXP_INF_NAN_EXPONENT;
  41 #endif // LIBC_FREXP_INF_NAN_EXPONENT
  42     return x;
  43   }
  44   if (bits.is_zero()) {
  45     exp = 0;
  46     return x;
  47   }
  48
  49   NormalFloat<T> normal(bits);
  50   exp = normal.exponent + 1;
  51   normal.exponent = -1;
  52   return normal;
  53 }
  54
  55 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
  56 LIBC_INLINE T modf(T x, T &iptr) {
  57   FPBits<T> bits(x);
  58   if (bits.is_zero() || bits.is_nan()) {
  59     iptr = x;
  60     return x;
  61   } else if (bits.is_inf()) {
  62     iptr = x;
  63     return FPBits<T>::zero(bits.sign()).get_val();
  64   } else {
  65     iptr = trunc(x);
  66     if (x == iptr) {
  67       // If x is already an integer value, then return zero with the right
  68       // sign.
  69       return FPBits<T>::zero(bits.sign()).get_val();
  70     } else {
  71       return x - iptr;
  72     }
  73   }
  74 }
  75
  76 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
  77 LIBC_INLINE T copysign(T x, T y) {
  78   FPBits<T> xbits(x);
  79   xbits.set_sign(FPBits<T>(y).sign());
  80   return xbits.get_val();
  81 }
  82
  83 template <typename T> struct IntLogbConstants;
  84
  85 template <> struct IntLogbConstants<int> {
  86   LIBC_INLINE_VAR static constexpr int FP_LOGB0 = FP_ILOGB0;
  87   LIBC_INLINE_VAR static constexpr int FP_LOGBNAN = FP_ILOGBNAN;
  88   LIBC_INLINE_VAR static constexpr int T_MAX = INT_MAX;
  89   LIBC_INLINE_VAR static constexpr int T_MIN = INT_MIN;
  90 };
  91
  92 template <> struct IntLogbConstants<long> {
  93   LIBC_INLINE_VAR static constexpr long FP_LOGB0 = FP_ILOGB0;
  94   LIBC_INLINE_VAR static constexpr long FP_LOGBNAN = FP_ILOGBNAN;
  95   LIBC_INLINE_VAR static constexpr long T_MAX = LONG_MAX;
  96   LIBC_INLINE_VAR static constexpr long T_MIN = LONG_MIN;
  97 };
  98
  99 template <typename T, typename U>
 100 LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<U>, T>
 101 intlogb(U x) {
 102   FPBits<U> bits(x);
 103   if (LIBC_UNLIKELY(bits.is_zero() || bits.is_inf_or_nan())) {
 104     set_errno_if_required(EDOM);
 105     raise_except_if_required(FE_INVALID);
 106
 107     if (bits.is_zero())
 108       return IntLogbConstants<T>::FP_LOGB0;
 109     if (bits.is_nan())
 110       return IntLogbConstants<T>::FP_LOGBNAN;
 111     // bits is inf.
 112     return IntLogbConstants<T>::T_MAX;
 113   }
 114
 115   DyadicFloat<FPBits<U>::STORAGE_LEN> normal(bits.get_val());
 116   int exponent = normal.get_unbiased_exponent();
 117   // The C standard does not specify the return value when an exponent is
 118   // out of int range. However, XSI conformance required that INT_MAX or
 119   // INT_MIN are returned.
 120   // NOTE: It is highly unlikely that exponent will be out of int range as
 121   // the exponent is only 15 bits wide even for the 128-bit floating point
 122   // format.
 123   if (LIBC_UNLIKELY(exponent > IntLogbConstants<T>::T_MAX ||
 124                     exponent < IntLogbConstants<T>::T_MIN)) {
 125     set_errno_if_required(ERANGE);
 126     raise_except_if_required(FE_INVALID);
 127     return exponent > 0 ? IntLogbConstants<T>::T_MAX
 128                         : IntLogbConstants<T>::T_MIN;
 129   }
 130
 131   return static_cast<T>(exponent);
 132 }
 133
 134 template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
 135 LIBC_INLINE constexpr T logb(T x) {
 136   FPBits<T> bits(x);
 137   if (LIBC_UNLIKELY(bits.is_zero() || bits.is_inf_or_nan())) {
 138     if (bits.is_nan())
 139       return x;
 140
 141     raise_except_if_required(FE_DIVBYZERO);
 142
 143     if (bits.is_zero()) {
 144       set_errno_if_required(ERANGE);
 145       return FPBits<T>::inf(Sign::NEG).get_val();
 146     }
 147     // bits is inf.
 148     return FPBits<T>::inf().get_val();
 149   }
 150
 151   DyadicFloat<FPBits<T>::STORAGE_LEN> normal(bits.get_val());
 152   return static_cast<T>(normal.get_unbiased_exponent());
 153 }
 154
 155 template <typename T, typename U>
 156 LIBC_INLINE constexpr cpp::enable_if_t<
 157     cpp::is_floating_point_v<T> && cpp::is_integral_v<U>, T>
 158 ldexp(T x, U exp) {
 159   FPBits<T> bits(x);
 160   if (LIBC_UNLIKELY((exp == 0) || bits.is_zero() || bits.is_inf_or_nan()))
 161     return x;
 162
 163   // NormalFloat uses int32_t to store the true exponent value. We should ensure
 164   // that adding |exp| to it does not lead to integer rollover. But, if |exp|
 165   // value is larger the exponent range for type T, then we can return infinity
 166   // early. Because the result of the ldexp operation can be a subnormal number,
 167   // we need to accommodate the (mantissaWidth + 1) worth of shift in
 168   // calculating the limit.
 169   constexpr int EXP_LIMIT =
 170       FPBits<T>::MAX_BIASED_EXPONENT + FPBits<T>::FRACTION_LEN + 1;
 171   // Make sure that we can safely cast exp to int when not returning early.
 172   static_assert(EXP_LIMIT <= INT_MAX && -EXP_LIMIT >= INT_MIN);
 173   if (LIBC_UNLIKELY(exp > EXP_LIMIT)) {
 174     int rounding_mode = quick_get_round();
 175     Sign sign = bits.sign();
 176
 177     if ((sign == Sign::POS && rounding_mode == FE_DOWNWARD) ||
 178         (sign == Sign::NEG && rounding_mode == FE_UPWARD) ||
 179         (rounding_mode == FE_TOWARDZERO))
 180       return FPBits<T>::max_normal(sign).get_val();
 181
 182     set_errno_if_required(ERANGE);
 183     raise_except_if_required(FE_OVERFLOW);
 184     return FPBits<T>::inf(sign).get_val();
 185   }
 186
 187   // Similarly on the negative side we return zero early if |exp| is too small.
 188   if (LIBC_UNLIKELY(exp < -EXP_LIMIT)) {
 189     int rounding_mode = quick_get_round();
 190     Sign sign = bits.sign();
 191
 192     if ((sign == Sign::POS && rounding_mode == FE_UPWARD) ||
 193         (sign == Sign::NEG && rounding_mode == FE_DOWNWARD))
 194       return FPBits<T>::min_subnormal(sign).get_val();
 195
 196     set_errno_if_required(ERANGE);
 197     raise_except_if_required(FE_UNDERFLOW);
 198     return FPBits<T>::zero(sign).get_val();
 199   }
 200
 201   // For all other values, NormalFloat to T conversion handles it the right way.
 202   DyadicFloat<FPBits<T>::STORAGE_LEN> normal(bits.get_val());
 203   normal.exponent += static_cast<int>(exp);
 204   // TODO: Add tests for exceptions.
 205   return normal.template as<T, /*ShouldRaiseExceptions=*/true>();
 206 }
 207
 208 template <typename T, typename U,
 209           cpp::enable_if_t<cpp::is_floating_point_v<T> &&
 210                                cpp::is_floating_point_v<U> &&
 211                                (sizeof(T) <= sizeof(U)),
 212                            int> = 0>
 213 LIBC_INLINE T nextafter(T from, U to) {
 214   FPBits<T> from_bits(from);
 215   if (from_bits.is_nan())
 216     return from;
 217
 218   FPBits<U> to_bits(to);
 219   if (to_bits.is_nan())
 220     return cast<T>(to);
 221
 222   // NOTE: This would work only if `U` has a greater or equal precision than
 223   // `T`. Otherwise `from` could loose its precision and the following statement
 224   // could incorrectly evaluate to `true`.
 225   if (cast<U>(from) == to)
 226     return cast<T>(to);
 227
 228   using StorageType = typename FPBits<T>::StorageType;
 229   if (from != T(0)) {
 230     if ((cast<U>(from) < to) == (from > T(0))) {
 231       from_bits = FPBits<T>(StorageType(from_bits.uintval() + 1));
 232     } else {
 233       from_bits = FPBits<T>(StorageType(from_bits.uintval() - 1));
 234     }
 235   } else {
 236     from_bits = FPBits<T>::min_subnormal(to_bits.sign());
 237   }
 238
 239   if (from_bits.is_subnormal())
 240     raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
 241   else if (from_bits.is_inf())
 242     raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
 243
 244   return from_bits.get_val();
 245 }
 246
 247 template <bool IsDown, typename T,
 248           cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>
 249 LIBC_INLINE constexpr T nextupdown(T x) {
 250   constexpr Sign sign = IsDown ? Sign::NEG : Sign::POS;
 251
 252   FPBits<T> xbits(x);
 253   if (xbits.is_nan() || xbits == FPBits<T>::max_normal(sign) ||
 254       xbits == FPBits<T>::inf(sign))
 255     return x;
 256
 257   using StorageType = typename FPBits<T>::StorageType;
 258   if (x != T(0)) {
 259     if (xbits.sign() == sign) {
 260       xbits = FPBits<T>(StorageType(xbits.uintval() + 1));
 261     } else {
 262       xbits = FPBits<T>(StorageType(xbits.uintval() - 1));
 263     }
 264   } else {
 265     xbits = FPBits<T>::min_subnormal(sign);
 266   }
 267
 268   return xbits.get_val();
 269 }
 270
 271 } // namespace fputil
 272 } // namespace LIBC_NAMESPACE_DECL
 273
 274 #ifdef LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80
 275 #include "x86_64/NextAfterLongDouble.h"
 276 #include "x86_64/NextUpDownLongDouble.h"
 277 #endif // LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80
 278
 279 #endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H