//===-- Floating-point manipulation functions -------------------*- C++ -*-===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //===----------------------------------------------------------------------===// #ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H #define LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H #include "FPBits.h" #include "NearestIntegerOperations.h" #include "NormalFloat.h" #include "cast.h" #include "dyadic_float.h" #include "rounding_mode.h" #include "hdr/math_macros.h" #include "src/__support/CPP/bit.h" #include "src/__support/CPP/limits.h" // INT_MAX, INT_MIN #include "src/__support/CPP/type_traits.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/macros/attributes.h" #include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" // LIBC_UNLIKELY namespace LIBC_NAMESPACE_DECL { namespace fputil { template , int> = 0> LIBC_INLINE T frexp(T x, int &exp) { FPBits bits(x); if (bits.is_inf_or_nan()) { #ifdef LIBC_FREXP_INF_NAN_EXPONENT // The value written back to the second parameter when calling // frexp/frexpf/frexpl` with `+/-Inf`/`NaN` is unspecified in the standard. // Set the exp value for Inf/NaN inputs explicitly to // LIBC_FREXP_INF_NAN_EXPONENT if it is defined. exp = LIBC_FREXP_INF_NAN_EXPONENT; #endif // LIBC_FREXP_INF_NAN_EXPONENT return x; } if (bits.is_zero()) { exp = 0; return x; } NormalFloat normal(bits); exp = normal.exponent + 1; normal.exponent = -1; return normal; } template , int> = 0> LIBC_INLINE T modf(T x, T &iptr) { FPBits bits(x); if (bits.is_zero() || bits.is_nan()) { iptr = x; return x; } else if (bits.is_inf()) { iptr = x; return FPBits::zero(bits.sign()).get_val(); } else { iptr = trunc(x); if (x == iptr) { // If x is already an integer value, then return zero with the right // sign. return FPBits::zero(bits.sign()).get_val(); } else { return x - iptr; } } } template , int> = 0> LIBC_INLINE T copysign(T x, T y) { FPBits xbits(x); xbits.set_sign(FPBits(y).sign()); return xbits.get_val(); } template struct IntLogbConstants; template <> struct IntLogbConstants { LIBC_INLINE_VAR static constexpr int FP_LOGB0 = FP_ILOGB0; LIBC_INLINE_VAR static constexpr int FP_LOGBNAN = FP_ILOGBNAN; LIBC_INLINE_VAR static constexpr int T_MAX = INT_MAX; LIBC_INLINE_VAR static constexpr int T_MIN = INT_MIN; }; template <> struct IntLogbConstants { LIBC_INLINE_VAR static constexpr long FP_LOGB0 = FP_ILOGB0; LIBC_INLINE_VAR static constexpr long FP_LOGBNAN = FP_ILOGBNAN; LIBC_INLINE_VAR static constexpr long T_MAX = LONG_MAX; LIBC_INLINE_VAR static constexpr long T_MIN = LONG_MIN; }; template LIBC_INLINE constexpr cpp::enable_if_t, T> intlogb(U x) { FPBits bits(x); if (LIBC_UNLIKELY(bits.is_zero() || bits.is_inf_or_nan())) { set_errno_if_required(EDOM); raise_except_if_required(FE_INVALID); if (bits.is_zero()) return IntLogbConstants::FP_LOGB0; if (bits.is_nan()) return IntLogbConstants::FP_LOGBNAN; // bits is inf. return IntLogbConstants::T_MAX; } DyadicFloat::STORAGE_LEN> normal(bits.get_val()); int exponent = normal.get_unbiased_exponent(); // The C standard does not specify the return value when an exponent is // out of int range. However, XSI conformance required that INT_MAX or // INT_MIN are returned. // NOTE: It is highly unlikely that exponent will be out of int range as // the exponent is only 15 bits wide even for the 128-bit floating point // format. if (LIBC_UNLIKELY(exponent > IntLogbConstants::T_MAX || exponent < IntLogbConstants::T_MIN)) { set_errno_if_required(ERANGE); raise_except_if_required(FE_INVALID); return exponent > 0 ? IntLogbConstants::T_MAX : IntLogbConstants::T_MIN; } return static_cast(exponent); } template , int> = 0> LIBC_INLINE constexpr T logb(T x) { FPBits bits(x); if (LIBC_UNLIKELY(bits.is_zero() || bits.is_inf_or_nan())) { if (bits.is_nan()) return x; raise_except_if_required(FE_DIVBYZERO); if (bits.is_zero()) { set_errno_if_required(ERANGE); return FPBits::inf(Sign::NEG).get_val(); } // bits is inf. return FPBits::inf().get_val(); } DyadicFloat::STORAGE_LEN> normal(bits.get_val()); return static_cast(normal.get_unbiased_exponent()); } template LIBC_INLINE constexpr cpp::enable_if_t< cpp::is_floating_point_v && cpp::is_integral_v, T> ldexp(T x, U exp) { FPBits bits(x); if (LIBC_UNLIKELY((exp == 0) || bits.is_zero() || bits.is_inf_or_nan())) return x; // NormalFloat uses int32_t to store the true exponent value. We should ensure // that adding |exp| to it does not lead to integer rollover. But, if |exp| // value is larger the exponent range for type T, then we can return infinity // early. Because the result of the ldexp operation can be a subnormal number, // we need to accommodate the (mantissaWidth + 1) worth of shift in // calculating the limit. constexpr int EXP_LIMIT = FPBits::MAX_BIASED_EXPONENT + FPBits::FRACTION_LEN + 1; // Make sure that we can safely cast exp to int when not returning early. static_assert(EXP_LIMIT <= INT_MAX && -EXP_LIMIT >= INT_MIN); if (LIBC_UNLIKELY(exp > EXP_LIMIT)) { int rounding_mode = quick_get_round(); Sign sign = bits.sign(); if ((sign == Sign::POS && rounding_mode == FE_DOWNWARD) || (sign == Sign::NEG && rounding_mode == FE_UPWARD) || (rounding_mode == FE_TOWARDZERO)) return FPBits::max_normal(sign).get_val(); set_errno_if_required(ERANGE); raise_except_if_required(FE_OVERFLOW); return FPBits::inf(sign).get_val(); } // Similarly on the negative side we return zero early if |exp| is too small. if (LIBC_UNLIKELY(exp < -EXP_LIMIT)) { int rounding_mode = quick_get_round(); Sign sign = bits.sign(); if ((sign == Sign::POS && rounding_mode == FE_UPWARD) || (sign == Sign::NEG && rounding_mode == FE_DOWNWARD)) return FPBits::min_subnormal(sign).get_val(); set_errno_if_required(ERANGE); raise_except_if_required(FE_UNDERFLOW); return FPBits::zero(sign).get_val(); } // For all other values, NormalFloat to T conversion handles it the right way. DyadicFloat::STORAGE_LEN> normal(bits.get_val()); normal.exponent += static_cast(exp); // TODO: Add tests for exceptions. return normal.template as(); } template && cpp::is_floating_point_v && (sizeof(T) <= sizeof(U)), int> = 0> LIBC_INLINE T nextafter(T from, U to) { FPBits from_bits(from); if (from_bits.is_nan()) return from; FPBits to_bits(to); if (to_bits.is_nan()) return cast(to); // NOTE: This would work only if `U` has a greater or equal precision than // `T`. Otherwise `from` could loose its precision and the following statement // could incorrectly evaluate to `true`. if (cast(from) == to) return cast(to); using StorageType = typename FPBits::StorageType; if (from != T(0)) { if ((cast(from) < to) == (from > T(0))) { from_bits = FPBits(StorageType(from_bits.uintval() + 1)); } else { from_bits = FPBits(StorageType(from_bits.uintval() - 1)); } } else { from_bits = FPBits::min_subnormal(to_bits.sign()); } if (from_bits.is_subnormal()) raise_except_if_required(FE_UNDERFLOW | FE_INEXACT); else if (from_bits.is_inf()) raise_except_if_required(FE_OVERFLOW | FE_INEXACT); return from_bits.get_val(); } template , int> = 0> LIBC_INLINE constexpr T nextupdown(T x) { constexpr Sign sign = IsDown ? Sign::NEG : Sign::POS; FPBits xbits(x); if (xbits.is_nan() || xbits == FPBits::max_normal(sign) || xbits == FPBits::inf(sign)) return x; using StorageType = typename FPBits::StorageType; if (x != T(0)) { if (xbits.sign() == sign) { xbits = FPBits(StorageType(xbits.uintval() + 1)); } else { xbits = FPBits(StorageType(xbits.uintval() - 1)); } } else { xbits = FPBits::min_subnormal(sign); } return xbits.get_val(); } } // namespace fputil } // namespace LIBC_NAMESPACE_DECL #ifdef LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80 #include "x86_64/NextAfterLongDouble.h" #include "x86_64/NextUpDownLongDouble.h" #endif // LIBC_TYPES_LONG_DOUBLE_IS_X86_FLOAT80 #endif // LLVM_LIBC_SRC___SUPPORT_FPUTIL_MANIPULATIONFUNCTIONS_H