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 280