280 lines · c
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-exception6//7//===----------------------------------------------------------------------===//8 9#ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H10#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H11 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_MIN22#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_UNLIKELY27 28namespace LIBC_NAMESPACE_DECL {29namespace fputil {30 31template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>32LIBC_INLINE constexpr T frexp(T x, int &exp) {33 FPBits<T> bits(x);34 if (bits.is_inf_or_nan()) {35#ifdef LIBC_FREXP_INF_NAN_EXPONENT36 // The value written back to the second parameter when calling37 // frexp/frexpf/frexpl` with `+/-Inf`/`NaN` is unspecified in the standard.38 // Set the exp value for Inf/NaN inputs explicitly to39 // LIBC_FREXP_INF_NAN_EXPONENT if it is defined.40 exp = LIBC_FREXP_INF_NAN_EXPONENT;41#endif // LIBC_FREXP_INF_NAN_EXPONENT42 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 55template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>56LIBC_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 right68 // sign.69 return FPBits<T>::zero(bits.sign()).get_val();70 } else {71 return x - iptr;72 }73 }74}75 76template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>77LIBC_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 83template <typename T> struct IntLogbConstants;84 85template <> 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 92template <> 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 99template <typename T, typename U>100LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<U>, T>101intlogb(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 is118 // out of int range. However, XSI conformance required that INT_MAX or119 // INT_MIN are returned.120 // NOTE: It is highly unlikely that exponent will be out of int range as121 // the exponent is only 15 bits wide even for the 128-bit floating point122 // 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_MAX128 : IntLogbConstants<T>::T_MIN;129 }130 131 return static_cast<T>(exponent);132}133 134template <typename T, cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>135LIBC_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 155template <typename T, typename U>156LIBC_INLINE constexpr cpp::enable_if_t<157 cpp::is_floating_point_v<T> && cpp::is_integral_v<U>, T>158ldexp(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 ensure164 // 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 infinity166 // early. Because the result of the ldexp operation can be a subnormal number,167 // we need to accommodate the (mantissaWidth + 1) worth of shift in168 // 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 208template <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>213LIBC_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 than223 // `T`. Otherwise `from` could loose its precision and the following statement224 // 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 247template <bool IsDown, typename T,248 cpp::enable_if_t<cpp::is_floating_point_v<T>, int> = 0>249LIBC_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 fputil272} // namespace LIBC_NAMESPACE_DECL273 274#ifdef LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80275#include "x86_64/NextAfterLongDouble.h"276#include "x86_64/NextUpDownLongDouble.h"277#endif // LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80278 279#endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H280