[libc] Add platform independent floating point rounding mode checks.

[lntue]

Summary

Many math functions need to check for floating point rounding modes to
return correct values. Currently most of them use the internal implementation
of fegetround, which is platform-dependent and blocking math functions to be
enabled on platforms with unimplemented fegetround. In this change, we add
platform independent rounding mode checks and switching math functions to use
them instead. https://github.com/llvm/llvm-project/issues/63016

[michaelrj]

libc/src/__support/FPUtil/FEnvImpl.h

130

I like the idea, but I'd like to have a flag to use the old method as well for targets where we're using -ffastmath or other optimizations that might precompute this value unexpectedly.

[lntue]

libc/src/__support/FPUtil/FEnvImpl.h

130

For math functions, I think it is ok be not correctly rounded in all rounding modes with -ffastmath, but ideally, string functions should not be affected by it. On the other hand, -ffast-math flags break the floating point standards in many ways, such as flush denormal numbers to 0, so maybe losing the last bit of accuracy in some rounding modes with -ffast-math is acceptable?

[sivachandra]

libc/src/__support/FPUtil/FEnvImpl.h

130

What do you think of this:

1. Eliminate get_round from FEnvImpl.h.
2. Move these new functions to FPUtil/rounding_mode.h.
3. Make rounding_mode.h::get_round the only implementation of get_round in the libc.

[lntue]

libc/src/__support/FPUtil/FEnvImpl.h

130

I that would be great, but unfortunately I think fegetround still have to return the correct one from fesetround, even with -ffast-math and fno-rounding-math: https://godbolt.org/z/8ca6djYvh

So we still need the current get_round somewhere. Only the internal functions that are not expected to be correctly rounded with -ffast-math or -fno-rounding-mode should use this function.

[sivachandra]

libc/src/__support/FPUtil/FEnvImpl.h

130

I do not know enough about -ffast-math and friends that I can agree or disagree with this. But, I think moving this to a separate header file is still good for distinction of what this FEnvImpl.h is for (it an implementation of a uniform interface to observe and manipulate the hardware FP env).

[lntue]

Move the new functions to a separate header.

[lntue]

libc/src/__support/FPUtil/FEnvImpl.h

130

I've moved these functions to rounding_mode.h

libc/src/__support/FPUtil/FEnvImpl.h

LIBC_INLINE int raise_except_if_required(int excepts) {
   return 0;
 }
 
 LIBC_INLINE void set_errno_if_required(int err) {
   if (math_errhandling & MATH_ERRNO)
     libc_errno = err;
 }
 
+// Quick free-standing test whether fegetround() == FE_UPWARD.
+// Using the following observation:
+//   1.0f + 2^-25 = 1.0f        for FE_TONEAREST, FE_DOWNWARD, FE_TOWARDZERO
+//                = 0x1.000002f for FE_UPWARD.
+LIBC_INLINE bool fenv_is_round_up() {
+  volatile float x = 0x1.0p-25f;
+  return (1.0f + x != 1.0f);
+}
+
+// Quick free-standing test whether fegetround() == FE_DOWNWARD.
+// Using the following observation:
+//   -1.0f - 2^-25 = -1.0f        for FE_TONEAREST, FE_UPWARD, FE_TOWARDZERO
+//                 = -0x1.000002f for FE_DOWNWARD.
+LIBC_INLINE bool fenv_is_round_down() {
+  volatile float x = 0x1.0p-25f;
+  return (-1.0f - x != -1.0f);
+}
+
+// Quick free-standing test whether fegetround() == FE_TONEAREST.
+// Using the following observation:
+//   1.5f + 2^-24 = 1.5f           for FE_TONEAREST, FE_DOWNWARD, FE_TOWARDZERO
+//                = 0x1.100002p0f  for FE_UPWARD,
+//   1.5f - 2^-24 = 1.5f           for FE_TONEAREST, FE_UPWARD
+//                = 0x1.0ffffep-1f for FE_DOWNWARD, FE_TOWARDZERO
+LIBC_INLINE bool fenv_is_round_to_nearest() {
+  static volatile float x = 0x1.0p-24f;
+  float y = x;
+  return (1.5f + y == 1.5f - y);
+}
+
+// Quick free-standing test whether fegetround() == FE_TOWARDZERO.
+// Using the following observation:
+//   1.0f + 2^-23 + 2^-24 = 0x1.000002p0f for FE_DOWNWARD, FE_TOWARDZERO
+//                        = 0x1.000004p0f for FE_TONEAREST, FE_UPWARD,
+//  -1.0f - 2^-24 = -1.0f          for FE_TONEAREST, FE_UPWARD, FE_TOWARDZERO
+//                = -0x1.000002p0f for FE_DOWNWARD
+// So:
+// (0x1.000002p0f + 2^-24) + (-1.0f - 2^-24) = 2^-23 for FE_TOWARDZERO
+//                                           = 2^-22 for FE_TONEAREST, FE_UPWARD
+//                                           = 0 for FE_DOWNWARD
+LIBC_INLINE bool fenv_is_round_to_zero() {
+  static volatile float x = 0x1.0p-24f;
+  float y = x;
+  return ((0x1.000002p0f + y) + (-1.0f - y) == 0x1.0p-23f);
+}
+
+// Quick free standing get rounding mode based on the above observations.
+LIBC_INLINE int quick_get_round() {
+  static volatile float x = 0x1.0p-24f;
+  float y = x;
+  float z = (0x1.000002p0f + y) + (-1.0f - y);
+
+  if (z == 0.0f)
+    return FE_DOWNWARD;
+  if (z == 0x1.0p-23f)
+    return FE_TOWARDZERO;
+  return (2.0f + y == 2.0f) ? FE_TONEAREST : FE_UPWARD;
+}
+
 } // namespace __llvm_libc::fputil
 
 #endif // LLVM_LIBC_SRC_SUPPORT_FPUTIL_FENVIMPL_H

libc/src/__support/FPUtil/Hypot.h

if (sum >= (DUIntType(1) << (2 * y_mant_width + 2))) {
     // a^2 + b^2 >= 4* leading_one^2, so we will need an extra bit to the left.
     if (leading_one == ONE) {
       // For normal result, we discard the last 2 bits of the sum and increase
       // the exponent.
       sticky_bits = sticky_bits || ((sum & 0x3U) != 0);
       sum >>= 2;
       ++out_exp;
       if (out_exp >= FPBits_t::MAX_EXPONENT) {
-        if (int round_mode = get_round();
+        if (int round_mode = quick_get_round();
             round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
           return T(FPBits_t::inf());
         return T(FPBits_t(FPBits_t::MAX_NORMAL));
       }
     } else {
       // For denormal result, we simply move the leading bit of the result to
       // the left by 1.
       leading_one <<= 1;
if (y_new >= ONE) {
     if (out_exp == 0) {
       out_exp = 1;
     }
   }
 
   y_new >>= 1;
 
   // Round to the nearest, tie to even.
-  int round_mode = get_round();
+  int round_mode = quick_get_round();
   switch (round_mode) {
   case FE_TONEAREST:
     // Round to nearest, ties to even
     if (round_bit && (lsb || sticky_bits || (r != 0)))
       ++y_new;
     break;
   case FE_UPWARD:
     if (round_bit || sticky_bits || (r != 0))

libc/src/__support/FPUtil/NearestIntegerOperations.h

LIBC_INLINE T round_using_current_rounding_mode(T x) {
   FPBits<T> bits(x);
 
   // If x is infinity NaN or zero, return it.
   if (bits.is_inf_or_nan() || bits.is_zero())
     return x;
 
   bool is_neg = bits.get_sign();
   int exponent = bits.get_exponent();
-  int rounding_mode = get_round();
+  int rounding_mode = quick_get_round();
 
   // If the exponent is greater than the most negative mantissa
   // exponent, then x is already an integer.
   if (exponent >= static_cast<int>(MantissaWidth<T>::VALUE))
     return x;
 
   if (exponent <= -1) {
     switch (rounding_mode) {

libc/src/__support/FPUtil/except_value_utils.h

template <typename T, size_t N> struct ExceptValues {
   };
 
   Mapping values[N];
 
   LIBC_INLINE constexpr cpp::optional<T> lookup(UIntType x_bits) const {
     for (size_t i = 0; i < N; ++i) {
       if (LIBC_UNLIKELY(x_bits == values[i].input)) {
         UIntType out_bits = values[i].rnd_towardzero_result;
-        switch (fputil::get_round()) {
+        switch (fputil::quick_get_round()) {
         case FE_UPWARD:
           out_bits += values[i].rnd_upward_offset;
           break;
         case FE_DOWNWARD:
           out_bits += values[i].rnd_downward_offset;
           break;
         case FE_TONEAREST:
           out_bits += values[i].rnd_tonearest_offset;
           break;
         }
         return FPBits<T>(out_bits).get_val();
       }
     }
     return cpp::nullopt;
   }
 
   LIBC_INLINE constexpr cpp::optional<T> lookup_odd(UIntType x_abs,
                                                     bool sign) const {
     for (size_t i = 0; i < N; ++i) {
       if (LIBC_UNLIKELY(x_abs == values[i].input)) {
         UIntType out_bits = values[i].rnd_towardzero_result;
-        switch (fputil::get_round()) {
+        switch (fputil::quick_get_round()) {
         case FE_UPWARD:
           out_bits += sign ? values[i].rnd_downward_offset
                            : values[i].rnd_upward_offset;
           break;
         case FE_DOWNWARD:
           out_bits += sign ? values[i].rnd_upward_offset
                            : values[i].rnd_downward_offset;
           break;

libc/src/__support/FPUtil/generic/FMA.h

if (r_exp > 0) {
       r_exp = 0;
     }
   } else {
     // Return +0.0 when there is exact cancellation, i.e., x*y == -z exactly.
     prod_sign = false;
   }
 
   // Finalize the result.
-  int round_mode = fputil::get_round();
+  int round_mode = fputil::quick_get_round();
   if (LIBC_UNLIKELY(r_exp >= FPBits::MAX_EXPONENT)) {
     if ((round_mode == FE_TOWARDZERO) ||
         (round_mode == FE_UPWARD && prod_sign) ||
         (round_mode == FE_DOWNWARD && !prod_sign)) {
       result = FPBits::MAX_NORMAL;
       return prod_sign ? -cpp::bit_cast<double>(result)
                        : cpp::bit_cast<double>(result);
     }

libc/src/__support/FPUtil/generic/sqrt.h

if (bits.is_inf_or_nan()) {
         rb = true;
       }
 
       // Remove hidden bit and append the exponent field.
       x_exp = ((x_exp >> 1) + FPBits<T>::EXPONENT_BIAS);
 
       y = (y - ONE) | (static_cast<UIntType>(x_exp) << MantissaWidth<T>::VALUE);
 
-      switch (get_round()) {
+      switch (quick_get_round()) {
       case FE_TONEAREST:
         // Round to nearest, ties to even
         if (rb && (lsb || (r != 0)))
           ++y;
         break;
       case FE_UPWARD:
         if (rb || (r != 0))
           ++y;

libc/src/__support/FPUtil/generic/sqrt_80_bit_long_double.h

if (r >= tmp) {
       rb = true;
     }
 
     // Append the exponent field.
     x_exp = ((x_exp >> 1) + FPBits<long double>::EXPONENT_BIAS);
     y |= (static_cast<UIntType>(x_exp)
           << (MantissaWidth<long double>::VALUE + 1));
 
-    switch (get_round()) {
+    switch (quick_get_round()) {
     case FE_TONEAREST:
       // Round to nearest, ties to even
       if (rb && (lsb || (r != 0)))
         ++y;
       break;
     case FE_UPWARD:
       if (rb || (r != 0))
         ++y;

libc/src/__support/str_to_float.h

if (isdigit(src[index]) || src[index] == DECIMAL_POINT) { // regular number
     if (is_float_hex_start(src + index, DECIMAL_POINT)) {
       base = 16;
       index += 2;
       seen_digit = true;
     }
 
     RoundDirection round_direction = RoundDirection::Nearest;
 
-    switch (fputil::get_round()) {
+    switch (fputil::quick_get_round()) {
     case FE_TONEAREST:
       round_direction = RoundDirection::Nearest;
       break;
     case FE_UPWARD:
       if (sign == '+') {
         round_direction = RoundDirection::Up;
       } else {
         round_direction = RoundDirection::Down;

libc/src/math/generic/atanf.cpp

LLVM_LIBC_FUNCTION(float, atanf, (float x)) {
   if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
     if (xbits.is_inf())
       return static_cast<float>(opt_barrier(sign ? -M_MATH_PI_2 : M_MATH_PI_2));
     else
       return x;
   }
   // |x| == 0.06905200332403183
   if (LIBC_UNLIKELY(xbits.uintval() == 0x3d8d6b23U)) {
-    if (fputil::get_round() == FE_TONEAREST) {
+    if (fputil::fenv_is_round_to_nearest()) {
       // 0.06894256919622421
       FPBits br(0x3d8d31c3U);
       br.set_sign(sign);
       return br.get_val();
     }
   }
 
   // |x| == 1.8670953512191772
   if (LIBC_UNLIKELY(xbits.uintval() == 0x3feefcfbU)) {
-    int rounding_mode = fputil::get_round();
+    int rounding_mode = fputil::quick_get_round();
     if (sign) {
       if (rounding_mode == FE_DOWNWARD) {
         // -1.0790828466415405
         return FPBits(0xbf8a1f63U).get_val();
       }
     } else {
       if (rounding_mode == FE_UPWARD) {
         // 1.0790828466415405

libc/src/math/generic/coshf.cpp

if (LIBC_UNLIKELY(x_u <= 0x3280'0000U)) {
     return 1.0f + x;
   }
 
   // When |x| >= 90, or x is inf or nan
   if (LIBC_UNLIKELY(x_u >= 0x42b4'0000U)) {
     if (xbits.is_inf_or_nan())
       return x + FPBits::inf().get_val();
 
-    int rounding = fputil::get_round();
+    int rounding = fputil::quick_get_round();
     if (LIBC_UNLIKELY(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO))
       return FPBits(FPBits::MAX_NORMAL).get_val();
 
     fputil::set_errno_if_required(ERANGE);
     fputil::raise_except_if_required(FE_OVERFLOW);
 
     return x + FPBits::inf().get_val();
   }

libc/src/math/generic/exp10f.cpp

if (LIBC_UNLIKELY(x_abs >= 0x421a'209bU)) {
     // When x < log10(2^-150) or nan
     if (x_u > 0xc234'9e35U) {
       // exp(-Inf) = 0
       if (xbits.is_inf())
         return 0.0f;
       // exp(nan) = nan
       if (xbits.is_nan())
         return x;
-      if (fputil::get_round() == FE_UPWARD)
+      if (fputil::fenv_is_round_up())
         return static_cast<float>(FPBits(FPBits::MIN_SUBNORMAL));
       fputil::set_errno_if_required(ERANGE);
       fputil::raise_except_if_required(FE_UNDERFLOW);
       return 0.0f;
     }
     // x >= log10(2^128) or nan
     if (!xbits.get_sign() && (x_u >= 0x421a'209bU)) {
       // x is finite
       if (x_u < 0x7f80'0000U) {
-        int rounding = fputil::get_round();
+        int rounding = fputil::quick_get_round();
         if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
           return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
 
         fputil::set_errno_if_required(ERANGE);
         fputil::raise_except_if_required(FE_OVERFLOW);
       }
       // x is +inf or nan
       return x + static_cast<float>(FPBits::inf());
     }
   }
 
   // When |x| <= log10(2)*2^-6
   if (LIBC_UNLIKELY(x_abs <= 0x3b9a'209bU)) {
     if (LIBC_UNLIKELY(x_u == 0xb25e'5bd9U)) { // x = -0x1.bcb7b2p-27f
-      if (fputil::get_round() == FE_TONEAREST)
+      if (fputil::fenv_is_round_to_nearest())
         return 0x1.fffffep-1f;
     }
     // |x| < 2^-25
     // 10^x ~ 1 + log(10) * x
     if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
       return fputil::multiply_add(x, 0x1.26bb1cp+1f, 1.0f);
     }
 
     return static_cast<float>(Exp10Base::powb_lo(x));
   }
 
   // Exceptional value.
   if (LIBC_UNLIKELY(x_u == 0x3d14'd956U)) { // x = 0x1.29b2acp-5f
-    if (fputil::get_round() == FE_UPWARD)
+    if (fputil::fenv_is_round_up())
       return 0x1.1657c4p+0f;
   }
 
   // Exact outputs when x = 1, 2, ..., 10.
   // Quick check mask: 0x800f'ffffU = ~(bits of 1.0f | ... | bits of 10.0f)
   if (LIBC_UNLIKELY((x_u & 0x800f'ffffU) == 0)) {
     switch (x_u) {
     case 0x3f800000U: // x = 1.0f

libc/src/math/generic/exp2f.cpp

LLVM_LIBC_FUNCTION(float, exp2f, (float x)) {
 
   // // When |x| >= 128, or x is nan
   if (LIBC_UNLIKELY(x_abs >= 0x4300'0000U)) {
 
     // x >= 128
     if (!xbits.get_sign()) {
       // x is finite
       if (x_u < 0x7f80'0000U) {
-        int rounding = fputil::get_round();
+        int rounding = fputil::quick_get_round();
         if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
           return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
 
         fputil::set_errno_if_required(ERANGE);
         fputil::raise_except_if_required(FE_OVERFLOW);
       }
       // x is +inf or nan
       return x + FPBits::inf().get_val();
     }
     // x <= -150
     if (x_u >= 0xc316'0000U) {
       // exp(-Inf) = 0
       if (xbits.is_inf())
         return 0.0f;
       // exp(nan) = nan
       if (xbits.is_nan())
         return x;
-      if (fputil::get_round() == FE_UPWARD)
+      if (fputil::fenv_is_round_up())
         return FPBits(FPBits::MIN_SUBNORMAL).get_val();
       if (x != 0.0f) {
         fputil::set_errno_if_required(ERANGE);
         fputil::raise_except_if_required(FE_UNDERFLOW);
       }
       return 0.0f;
     }
   }
 
   // Check exceptional values.
   if (LIBC_UNLIKELY((x_u & EXVAL_MASK) == EXVAL_MASK)) {
     if (LIBC_UNLIKELY(x_u == EXVAL1)) { // x = 0x1.853a6ep-9f
-      if (fputil::get_round() == FE_TONEAREST)
+      if (fputil::fenv_is_round_to_nearest())
         return 0x1.00870ap+0f;
     } else if (LIBC_UNLIKELY(x_u == EXVAL2)) { // x = -0x1.e7526ep-6f
-      if (fputil::get_round() == FE_TONEAREST)
+      if (fputil::fenv_is_round_to_nearest())
         return 0x1.f58d62p-1f;
     }
   }
 
   // For -150 < x < 128, to compute 2^x, we perform the following range
   // reduction: find hi, mid, lo such that:
   //   x = hi + mid + lo, in which
   //     hi is an integer,

libc/src/math/generic/expf.cpp

if (LIBC_UNLIKELY(x_abs >= 0x42b2'0000U || x_abs <= 0x3280'0000U)) {
     // When x < log(2^-150) or nan
     if (xbits.uintval() >= 0xc2cf'f1b5U) {
       // exp(-Inf) = 0
       if (xbits.is_inf())
         return 0.0f;
       // exp(nan) = nan
       if (xbits.is_nan())
         return x;
-      if (fputil::get_round() == FE_UPWARD)
+      if (fputil::fenv_is_round_up())
         return static_cast<float>(FPBits(FPBits::MIN_SUBNORMAL));
       fputil::set_errno_if_required(ERANGE);
       fputil::raise_except_if_required(FE_UNDERFLOW);
       return 0.0f;
     }
     // x >= 89 or nan
     if (!xbits.get_sign() && (xbits.uintval() >= 0x42b2'0000)) {
       // x is finite
       if (xbits.uintval() < 0x7f80'0000U) {
-        int rounding = fputil::get_round();
+        int rounding = fputil::quick_get_round();
         if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
           return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
 
         fputil::set_errno_if_required(ERANGE);
         fputil::raise_except_if_required(FE_OVERFLOW);
       }
       // x is +inf or nan
       return x + static_cast<float>(FPBits::inf());

libc/src/math/generic/expm1f.cpp

LLVM_LIBC_FUNCTION(float, expm1f, (float x)) {
   using FPBits = typename fputil::FPBits<float>;
   FPBits xbits(x);
 
   uint32_t x_u = xbits.uintval();
   uint32_t x_abs = x_u & 0x7fff'ffffU;
 
   // Exceptional value
   if (LIBC_UNLIKELY(x_u == 0x3e35'bec5U)) { // x = 0x1.6b7d8ap-3f
-    int round_mode = fputil::get_round();
+    int round_mode = fputil::quick_get_round();
     if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
       return 0x1.8dbe64p-3f;
     return 0x1.8dbe62p-3f;
   }
 
 #if !defined(LIBC_TARGET_CPU_HAS_FMA)
   if (LIBC_UNLIKELY(x_u == 0xbdc1'c6cbU)) { // x = -0x1.838d96p-4f
-    int round_mode = fputil::get_round();
+    int round_mode = fputil::quick_get_round();
     if (round_mode == FE_TONEAREST || round_mode == FE_DOWNWARD)
       return -0x1.71c884p-4f;
     return -0x1.71c882p-4f;
   }
 #endif // LIBC_TARGET_CPU_HAS_FMA
 
   // When |x| > 25*log(2), or nan
   if (LIBC_UNLIKELY(x_abs >= 0x418a'a123U)) {
     // x < log(2^-25)
     if (xbits.get_sign()) {
       // exp(-Inf) = 0
       if (xbits.is_inf())
         return -1.0f;
       // exp(nan) = nan
       if (xbits.is_nan())
         return x;
-      int round_mode = fputil::get_round();
+      int round_mode = fputil::quick_get_round();
       if (round_mode == FE_UPWARD || round_mode == FE_TOWARDZERO)
         return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f
       return -1.0f;
     } else {
       // x >= 89 or nan
       if (xbits.uintval() >= 0x42b2'0000) {
         if (xbits.uintval() < 0x7f80'0000U) {
-          int rounding = fputil::get_round();
+          int rounding = fputil::quick_get_round();
           if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
             return static_cast<float>(FPBits(FPBits::MAX_NORMAL));
 
           fputil::set_errno_if_required(ERANGE);
           fputil::raise_except_if_required(FE_OVERFLOW);
         }
         return x + static_cast<float>(FPBits::inf());
       }

libc/src/math/generic/sincosf.cpp

#endif // LIBC_TARGET_CPU_HAS_FMA
   }
 
   // Check exceptional values.
   for (int i = 0; i < N_EXCEPTS; ++i) {
     if (LIBC_UNLIKELY(x_abs == EXCEPT_INPUTS[i])) {
       uint32_t s = EXCEPT_OUTPUTS_SIN[i][0]; // FE_TOWARDZERO
       uint32_t c = EXCEPT_OUTPUTS_COS[i][0]; // FE_TOWARDZERO
       bool x_sign = x < 0;
-      switch (fputil::get_round()) {
+      switch (fputil::quick_get_round()) {
       case FE_UPWARD:
         s += x_sign ? EXCEPT_OUTPUTS_SIN[i][2] : EXCEPT_OUTPUTS_SIN[i][1];
         c += EXCEPT_OUTPUTS_COS[i][1];
         break;
       case FE_DOWNWARD:
         s += x_sign ? EXCEPT_OUTPUTS_SIN[i][1] : EXCEPT_OUTPUTS_SIN[i][2];
         c += EXCEPT_OUTPUTS_COS[i][2];
         break;

libc/src/math/generic/sinf.cpp

#endif // LIBC_TARGET_CPU_HAS_FMA
     double result =
         fputil::polyeval(xsq, 1.0, -0x1.55555555554c6p-3, 0x1.1111111085e65p-7,
                          -0x1.a019f70fb4d4fp-13, 0x1.718d179815e74p-19);
     return static_cast<float>(xd * result);
   }
 
   if (LIBC_UNLIKELY(x_abs == 0x4619'9998U)) { // x = 0x1.33333p13
     float r = -0x1.63f4bap-2f;
-    int rounding = fputil::get_round();
+    int rounding = fputil::quick_get_round();
     bool sign = xbits.get_sign();
     if ((rounding == FE_DOWNWARD && !sign) || (rounding == FE_UPWARD && sign))
       r = -0x1.63f4bcp-2f;
     return xbits.get_sign() ? -r : r;
   }
 
   if (LIBC_UNLIKELY(x_abs >= 0x7f80'0000U)) {
     if (x_abs == 0x7f80'0000U) {

libc/src/math/generic/sinhf.cpp

LLVM_LIBC_FUNCTION(float, sinhf, (float x)) {
   // When |x| >= 90, or x is inf or nan
   if (LIBC_UNLIKELY(x_abs >= 0x42b4'0000U)) {
     if (xbits.is_nan())
       return x + 1.0f; // sNaN to qNaN + signal
 
     if (xbits.is_inf())
       return x;
 
-    int rounding = fputil::get_round();
+    int rounding = fputil::quick_get_round();
     if (sign) {
       if (LIBC_UNLIKELY(rounding == FE_UPWARD || rounding == FE_TOWARDZERO))
         return FPBits(FPBits::MAX_NORMAL | FPBits::FloatProp::SIGN_MASK)
             .get_val();
     } else {
       if (LIBC_UNLIKELY(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO))
         return FPBits(FPBits::MAX_NORMAL).get_val();
     }
 
     fputil::set_errno_if_required(ERANGE);
     fputil::raise_except_if_required(FE_OVERFLOW);
 
     return x + FPBits::inf(sign).get_val();
   }
 
   // |x| <= 0.078125
   if (LIBC_UNLIKELY(x_abs <= 0x3da0'0000U)) {
     // |x| = 0.0005589424981735646724700927734375
     if (LIBC_UNLIKELY(x_abs == 0x3a12'85ffU)) {
-      if (fputil::get_round() == FE_TONEAREST)
+      if (fputil::fenv_is_round_to_nearest())
         return x;
     }
 
     double xdbl = x;
     double x2 = xdbl * xdbl;
     // Sollya: fpminimax(sinh(x),[|3,5,7|],[|D...|],[-1/16-1/64;1/16+1/64],x);
     // Sollya output: x * (0x1p0 + x^0x1p1 * (0x1.5555555556583p-3 + x^0x1p1
     //                  * (0x1.111110d239f1fp-7

libc/src/math/generic/tanhf.cpp

if (LIBC_UNLIKELY(x_abs <= 0x3da0'0000U)) {
     // Pure Taylor series.
     double pe = fputil::polyeval(x2, 0.0, -0x1.5555555555555p-2,
                                  0x1.1111111111111p-3, -0x1.ba1ba1ba1ba1cp-5,
                                  0x1.664f4882c10fap-6, -0x1.226e355e6c23dp-7);
     return static_cast<float>(fputil::multiply_add(xdbl, pe, xdbl));
   }
 
   if (LIBC_UNLIKELY(xbits.bits == 0x4058'e0a3U)) {
-    if (fputil::get_round() == FE_DOWNWARD)
+    if (fputil::fenv_is_round_down())
       return FPBits(0x3f7f'6ad9U).get_val();
   }
 
   // Range reduction: e^(2x) = 2^(mid + hi) * e^lo
   auto ep = exp_b_range_reduc<ExpBase>(2.0f * x); // exp(2 * x)
   double r = ExpBase::powb_lo(ep.lo);
   // tanh(x) = (exp(2x) - 1) / (exp(2x) + 1)
 #if defined(LIBC_TARGET_CPU_HAS_FMA)

libc/src/stdio/printf_core/float_dec_converter.h

for (; i < blocks; ++i) {
         // Is m * 10^(additionalDigits + 1) / 2^(-exponent) integer?
         const int32_t requiredTwos =
             -(exponent - MANT_WIDTH) - static_cast<int32_t>(precision) - 1;
         const bool trailingZeros =
             requiredTwos <= 0 ||
             (requiredTwos < 60 &&
              multiple_of_power_of_2(float_bits.get_explicit_mantissa(),
                                     static_cast<uint32_t>(requiredTwos)));
-        switch (fputil::get_round()) {
+        switch (fputil::quick_get_round()) {
         case FE_TONEAREST:
           // Round to nearest, if it's exactly halfway then round to even.
           if (last_digit != 5) {
             round = last_digit > 5 ? RoundDirection::Up : RoundDirection::Down;
           } else {
             round = trailingZeros ? RoundDirection::Even : RoundDirection::Up;
           }
           break;
LIBC_INLINE int convert_float_dec_exp_typed(Writer *writer,
   // Is m * 10^(additionalDigits + 1) / 2^(-exponent) integer?
   const int32_t requiredTwos =
       -(exponent - MANT_WIDTH) - static_cast<int32_t>(precision) - 1;
   const bool trailingZeros =
       requiredTwos <= 0 ||
       (requiredTwos < 60 &&
        multiple_of_power_of_2(float_bits.get_explicit_mantissa(),
                               static_cast<uint32_t>(requiredTwos)));
-  switch (fputil::get_round()) {
+  switch (fputil::quick_get_round()) {
   case FE_TONEAREST:
     // Round to nearest, if it's exactly halfway then round to even.
     if (last_digit != 5) {
       round = last_digit > 5 ? RoundDirection::Up : RoundDirection::Down;
     } else {
       round = trailingZeros ? RoundDirection::Even : RoundDirection::Up;
     }
     break;
LIBC_INLINE int convert_float_dec_auto_typed(Writer *writer,
   const int32_t requiredTwos =
       -(exponent - MANT_WIDTH) - static_cast<int32_t>(exp_precision) - 1;
   // TODO: rename this variable to remove confusion with trailing_zeroes
   const bool trailingZeros =
       requiredTwos <= 0 ||
       (requiredTwos < 60 &&
        multiple_of_power_of_2(float_bits.get_explicit_mantissa(),
                               static_cast<uint32_t>(requiredTwos)));
-  switch (fputil::get_round()) {
+  switch (fputil::quick_get_round()) {
   case FE_TONEAREST:
     // Round to nearest, if it's exactly halfway then round to even.
     if (last_digit != 5) {
       round = last_digit > 5 ? RoundDirection::Up : RoundDirection::Down;
     } else {
       round = trailingZeros ? RoundDirection::Even : RoundDirection::Up;
     }
     break;

libc/src/stdio/printf_core/float_hex_converter.h

const size_t shift_amount =
         (mant_len - intended_digits) * BITS_IN_HEX_DIGIT;
 
     const MantissaInt truncated_bits =
         mantissa & ((MantissaInt(1) << shift_amount) - 1);
     const MantissaInt halfway_const = MantissaInt(1) << (shift_amount - 1);
 
     mantissa >>= shift_amount;
 
-    switch (fputil::get_round()) {
+    switch (fputil::quick_get_round()) {
     case FE_TONEAREST:
       // Round to nearest, if it's exactly halfway then round to even.
       if (truncated_bits > halfway_const)
         ++mantissa;
       else if (truncated_bits == halfway_const)
         mantissa = mantissa + (mantissa & 1);
       break;
     case FE_DOWNWARD:

libc/test/src/__support/FPUtil/CMakeLists.txt

add_libc_test(
   fpbits_test
   SUITE
     libc-fputil-tests
   SRCS
     fpbits_test.cpp
   DEPENDS
     libc.src.__support.FPUtil.fp_bits
 )
+
+add_fp_unittest(
+  fenvimpl_test
+  SUITE
+    libc-fputil-tests
+  SRCS
+    fenvimpl_test.cpp
+  DEPENDS
+    libc.src.__support.FPUtil.fenv_impl
+)

libc/test/src/__support/FPUtil/fenvimpl_test.cpp

+//===-- Unittests for the FEnvImpl ----------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <math.h>
+
+using __llvm_libc::testing::mpfr::ForceRoundingMode;
+using __llvm_libc::testing::mpfr::RoundingMode;
+
+TEST(LlvmLibcFEnvImplTest, QuickRoundingUpTest) {
+  using __llvm_libc::fputil::fenv_is_round_up;
+  {
+    ForceRoundingMode __r(RoundingMode::Upward);
+    ASSERT_TRUE(fenv_is_round_up());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Downward);
+    ASSERT_FALSE(fenv_is_round_up());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Nearest);
+    ASSERT_FALSE(fenv_is_round_up());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::TowardZero);
+    ASSERT_FALSE(fenv_is_round_up());
+  }
+}
+
+TEST(LlvmLibcFEnvImplTest, QuickRoundingDownTest) {
+  using __llvm_libc::fputil::fenv_is_round_down;
+  {
+    ForceRoundingMode __r(RoundingMode::Upward);
+    ASSERT_FALSE(fenv_is_round_down());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Downward);
+    ASSERT_TRUE(fenv_is_round_down());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Nearest);
+    ASSERT_FALSE(fenv_is_round_down());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::TowardZero);
+    ASSERT_FALSE(fenv_is_round_down());
+  }
+}
+
+TEST(LlvmLibcFEnvImplTest, QuickRoundingNearestTest) {
+  using __llvm_libc::fputil::fenv_is_round_to_nearest;
+  {
+    ForceRoundingMode __r(RoundingMode::Upward);
+    ASSERT_FALSE(fenv_is_round_to_nearest());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Downward);
+    ASSERT_FALSE(fenv_is_round_to_nearest());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Nearest);
+    ASSERT_TRUE(fenv_is_round_to_nearest());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::TowardZero);
+    ASSERT_FALSE(fenv_is_round_to_nearest());
+  }
+}
+
+TEST(LlvmLibcFEnvImplTest, QuickRoundingTowardZeroTest) {
+  using __llvm_libc::fputil::fenv_is_round_to_zero;
+  {
+    ForceRoundingMode __r(RoundingMode::Upward);
+    ASSERT_FALSE(fenv_is_round_to_zero());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Downward);
+    ASSERT_FALSE(fenv_is_round_to_zero());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Nearest);
+    ASSERT_FALSE(fenv_is_round_to_zero());
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::TowardZero);
+    ASSERT_TRUE(fenv_is_round_to_zero());
+  }
+}
+
+TEST(LlvmLibcFEnvImplTest, QuickGetRoundTest) {
+  using __llvm_libc::fputil::quick_get_round;
+  {
+    ForceRoundingMode __r(RoundingMode::Upward);
+    ASSERT_EQ(quick_get_round(), FE_UPWARD);
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Downward);
+    ASSERT_EQ(quick_get_round(), FE_DOWNWARD);
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::Nearest);
+    ASSERT_EQ(quick_get_round(), FE_TONEAREST);
+  }
+  {
+    ForceRoundingMode __r(RoundingMode::TowardZero);
+    ASSERT_EQ(quick_get_round(), FE_TOWARDZERO);
+  }
+}