[libc-commits] [libc] [libc][math] Adding LIBC_MATH_ALWAYS_ROUND_NEAREST option (PR #201154)
Hoàng Minh Thiên via libc-commits
libc-commits at lists.llvm.org
Tue Jun 2 21:27:01 PDT 2026
https://github.com/hmthien050209 updated https://github.com/llvm/llvm-project/pull/201154
>From 247e359193afebd93576b6caa3c503008f587fc3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Tue, 2 Jun 2026 21:40:23 +0700
Subject: [PATCH 1/4] feat: declared the macro and updated Hypot.h
---
libc/src/__support/FPUtil/Hypot.h | 16 ++++++++++++++++
libc/src/__support/FPUtil/rounding_mode.h | 6 ++++++
libc/src/__support/macros/optimization.h | 5 +++++
3 files changed, 27 insertions(+)
diff --git a/libc/src/__support/FPUtil/Hypot.h b/libc/src/__support/FPUtil/Hypot.h
index 292140065c754..480dbb3be195e 100644
--- a/libc/src/__support/FPUtil/Hypot.h
+++ b/libc/src/__support/FPUtil/Hypot.h
@@ -199,10 +199,14 @@ LIBC_INLINE T hypot(T x, T y) {
sum >>= 2;
++out_exp;
if (out_exp >= FPBits_t::MAX_BIASED_EXPONENT) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
if (int round_mode = quick_get_round();
round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
return FPBits_t::inf().get_val();
return FPBits_t::max_normal().get_val();
+#else
+ return FPBits_t::inf().get_val();
+#endif // LIBC_HAS_ALWAYS_ROUND_NEAREST
}
} else {
// For denormal result, we simply move the leading bit of the result to
@@ -241,6 +245,17 @@ LIBC_INLINE T hypot(T x, T y) {
y_new >>= 1;
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ // Round to the nearest, tie to even.
+ if (round_bit && (lsb || sticky_bits || (r != 0)))
+ ++y_new;
+ if (y_new >= (ONE >> 1)) {
+ y_new -= ONE >> 1;
+ ++out_exp;
+ if (out_exp >= FPBits_t::MAX_BIASED_EXPONENT)
+ return FPBits_t::inf().get_val();
+ }
+#else
// Round to the nearest, tie to even.
int round_mode = quick_get_round();
switch (round_mode) {
@@ -264,6 +279,7 @@ LIBC_INLINE T hypot(T x, T y) {
return FPBits_t::max_normal().get_val();
}
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
y_new |= static_cast<StorageType>(out_exp) << FPBits_t::FRACTION_LEN;
diff --git a/libc/src/__support/FPUtil/rounding_mode.h b/libc/src/__support/FPUtil/rounding_mode.h
index 6ce693d41da50..3a4a41d8beacc 100644
--- a/libc/src/__support/FPUtil/rounding_mode.h
+++ b/libc/src/__support/FPUtil/rounding_mode.h
@@ -67,6 +67,7 @@ LIBC_INLINE bool fenv_is_round_to_zero() {
// Quick free standing get rounding mode based on the above observations.
LIBC_INLINE int quick_get_round() {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
static volatile float x = 0x1.0p-24f;
float y = x;
float z = (0x1.000002p0f + y) + (-1.0f - y);
@@ -76,6 +77,11 @@ LIBC_INLINE int quick_get_round() {
if (z == 0x1.0p-23f)
return FE_TOWARDZERO;
return (2.0f + y == 2.0f) ? FE_TONEAREST : FE_UPWARD;
+#else
+ // the optimization should've happened in each user of this function,
+ // but just in case
+ return FE_TONEAREST;
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
}
} // namespace generic
diff --git a/libc/src/__support/macros/optimization.h b/libc/src/__support/macros/optimization.h
index 411653439c88c..ecfadc8a86ca0 100644
--- a/libc/src/__support/macros/optimization.h
+++ b/libc/src/__support/macros/optimization.h
@@ -56,6 +56,7 @@ LIBC_INLINE constexpr bool expects_bool_condition(T value, T expected) {
(LIBC_MATH_SKIP_ACCURATE_PASS | LIBC_MATH_SMALL_TABLES | \
LIBC_MATH_NO_ERRNO | LIBC_MATH_NO_EXCEPT)
#define LIBC_MATH_INTERMEDIATE_COMP_IN_FLOAT 0x10
+#define LIBC_MATH_ALWAYS_ROUND_NEAREST 0x20
#ifndef LIBC_MATH
#define LIBC_MATH 0
@@ -81,4 +82,8 @@ LIBC_INLINE constexpr bool expects_bool_condition(T value, T expected) {
#define LIBC_MATH_HAS_NO_EXCEPT
#endif
+#if ((LIBC_MATH) & LIBC_MATH_ALWAYS_ROUND_NEAREST)
+#define LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+#endif
+
#endif // LLVM_LIBC_SRC___SUPPORT_MACROS_OPTIMIZATION_H
>From 5b15a60869adc0a4338ebac23f7d283c10bb6467 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Tue, 2 Jun 2026 21:54:45 +0700
Subject: [PATCH 2/4] implemented on all all dependant of
`fputil::quick_get_round()`
---
.../__support/FPUtil/ManipulationFunctions.h | 3 +++
.../FPUtil/NearestIntegerOperations.h | 4 ++++
libc/src/__support/FPUtil/dyadic_float.h | 21 +++++++++++++++++++
.../src/__support/FPUtil/except_value_utils.h | 8 +++++++
libc/src/__support/FPUtil/generic/FMA.h | 6 +++++-
libc/src/__support/FPUtil/generic/add_sub.h | 8 +++++++
.../FPUtil/generic/sqrt_80_bit_long_double.h | 7 ++++++-
libc/src/__support/FPUtil/rounding_mode.h | 9 +++-----
libc/src/__support/math/asinbf16.h | 2 ++
libc/src/__support/math/asinf16.h | 2 ++
libc/src/__support/math/asinhf16.h | 2 ++
libc/src/__support/math/asinpi.h | 5 +++++
libc/src/__support/math/coshf.h | 2 ++
libc/src/__support/math/coshf16.h | 6 ++++++
libc/src/__support/math/exp.h | 4 ++++
libc/src/__support/math/exp10.h | 4 ++++
libc/src/__support/math/exp10f.h | 6 ++++++
libc/src/__support/math/exp10f16.h | 6 ++++++
libc/src/__support/math/exp10m1f.h | 13 ++++++++++++
libc/src/__support/math/exp10m1f16.h | 14 +++++++++++--
libc/src/__support/math/exp2.h | 4 ++++
libc/src/__support/math/exp2f.h | 2 ++
libc/src/__support/math/exp2f16.h | 6 ++++++
libc/src/__support/math/exp2m1f.h | 4 ++++
libc/src/__support/math/exp2m1f16.h | 10 +++++++++
libc/src/__support/math/expf.h | 2 ++
libc/src/__support/math/expf16.h | 10 +++++++++
libc/src/__support/math/expm1.h | 2 ++
libc/src/__support/math/expm1f.h | 8 +++++++
libc/src/__support/math/expm1f16.h | 6 ++++++
libc/src/__support/math/sin.h | 2 ++
libc/src/__support/math/sincos.h | 2 ++
libc/src/__support/math/sincosf.h | 5 +++++
libc/src/__support/math/sinf.h | 2 ++
libc/src/__support/math/sinf16.h | 7 +++++++
libc/src/__support/math/sinhf.h | 2 ++
libc/src/__support/math/sinhf16.h | 6 ++++++
libc/src/__support/math/tan.h | 3 +++
libc/src/__support/math/tanf16.h | 2 ++
libc/src/__support/math/tanhf16.h | 10 +++++++++
libc/src/__support/str_to_float.h | 4 ++++
.../stdio/printf_core/float_dec_converter.h | 8 +++++++
.../stdio/printf_core/float_hex_converter.h | 8 +++++++
43 files changed, 237 insertions(+), 10 deletions(-)
diff --git a/libc/src/__support/FPUtil/ManipulationFunctions.h b/libc/src/__support/FPUtil/ManipulationFunctions.h
index 16dc094deda47..a8c034cfaa118 100644
--- a/libc/src/__support/FPUtil/ManipulationFunctions.h
+++ b/libc/src/__support/FPUtil/ManipulationFunctions.h
@@ -170,6 +170,8 @@ ldexp(T x, U exp) {
FPBits<T>::MAX_BIASED_EXPONENT + FPBits<T>::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);
+
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
if (LIBC_UNLIKELY(exp > EXP_LIMIT)) {
int rounding_mode = quick_get_round();
Sign sign = bits.sign();
@@ -197,6 +199,7 @@ ldexp(T x, U exp) {
raise_except_if_required(FE_UNDERFLOW);
return FPBits<T>::zero(sign).get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
// For all other values, NormalFloat to T conversion handles it the right way.
DyadicFloat<FPBits<T>::STORAGE_LEN> normal(bits.get_val());
diff --git a/libc/src/__support/FPUtil/NearestIntegerOperations.h b/libc/src/__support/FPUtil/NearestIntegerOperations.h
index 4637b46bc9473..9329d2985ac6e 100644
--- a/libc/src/__support/FPUtil/NearestIntegerOperations.h
+++ b/libc/src/__support/FPUtil/NearestIntegerOperations.h
@@ -246,6 +246,10 @@ round_using_specific_rounding_mode(T x, int rnd) {
template <typename T>
LIBC_INLINE constexpr cpp::enable_if_t<cpp::is_floating_point_v<T>, T>
round_using_current_rounding_mode(T x) {
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ return round_using_specific_rounding_mode(x, FP_INT_TONEAREST);
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+
int rounding_mode = quick_get_round();
switch (rounding_mode) {
diff --git a/libc/src/__support/FPUtil/dyadic_float.h b/libc/src/__support/FPUtil/dyadic_float.h
index 7e4ea2f1a81be..4f436a71d2876 100644
--- a/libc/src/__support/FPUtil/dyadic_float.h
+++ b/libc/src/__support/FPUtil/dyadic_float.h
@@ -46,6 +46,18 @@ rounding_direction(const LIBC_NAMESPACE::UInt<Bits> &value, size_t rshift,
(rshift >= Bits && value == 0))
return 0; // exact
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ if (rshift > 0 && rshift <= Bits && value.get_bit(rshift - 1)) {
+ // We round up, unless the value is an exact halfway case and
+ // the bit that will end up in the units place is 0, in which
+ // case tie-break-to-even says round down.
+ bool round_bit = rshift < Bits ? value.get_bit(rshift) : 0;
+ return round_bit != 0 || (value << (Bits - rshift + 1)) != 0 ? +1 : -1;
+ } else {
+ return -1;
+ }
+#endif
+
switch (quick_get_round()) {
case FE_TONEAREST:
if (rshift > 0 && rshift <= Bits && value.get_bit(rshift - 1)) {
@@ -190,6 +202,10 @@ template <size_t Bits> struct DyadicFloat {
raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
}
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ return FPBits::inf(sign).get_val();
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+
switch (quick_get_round()) {
case FE_TONEAREST:
return FPBits::inf(sign).get_val();
@@ -248,6 +264,10 @@ template <size_t Bits> struct DyadicFloat {
StorageType result =
FPBits::create_value(sign, out_biased_exp, out_mantissa).uintval();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ if (round && (lsb || sticky))
+ ++result;
+#else
switch (quick_get_round()) {
case FE_TONEAREST:
if (round && (lsb || sticky))
@@ -264,6 +284,7 @@ template <size_t Bits> struct DyadicFloat {
default:
break;
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
if (ShouldSignalExceptions && (round || sticky)) {
int excepts = FE_INEXACT;
diff --git a/libc/src/__support/FPUtil/except_value_utils.h b/libc/src/__support/FPUtil/except_value_utils.h
index 37ecc52784c16..a125c3aafa827 100644
--- a/libc/src/__support/FPUtil/except_value_utils.h
+++ b/libc/src/__support/FPUtil/except_value_utils.h
@@ -58,6 +58,9 @@ template <typename T, size_t N> struct ExceptValues {
for (size_t i = 0; i < N; ++i) {
if (LIBC_UNLIKELY(x_bits == values[i].input)) {
StorageType out_bits = values[i].rnd_towardzero_result;
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ out_bits += values[i].rnd_tonearest_offset;
+#else
switch (fputil::quick_get_round()) {
case FE_UPWARD:
out_bits += values[i].rnd_upward_offset;
@@ -69,6 +72,7 @@ template <typename T, size_t N> struct ExceptValues {
out_bits += values[i].rnd_tonearest_offset;
break;
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
return FPBits<T>(out_bits).get_val();
}
}
@@ -80,6 +84,9 @@ template <typename T, size_t N> struct ExceptValues {
for (size_t i = 0; i < N; ++i) {
if (LIBC_UNLIKELY(x_abs == values[i].input)) {
StorageType out_bits = values[i].rnd_towardzero_result;
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ out_bits += values[i].rnd_tonearest_offset;
+#else
switch (fputil::quick_get_round()) {
case FE_UPWARD:
if (sign)
@@ -100,6 +107,7 @@ template <typename T, size_t N> struct ExceptValues {
out_bits += values[i].rnd_tonearest_offset;
break;
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
T result = FPBits<T>(out_bits).get_val();
if (sign)
result = -result;
diff --git a/libc/src/__support/FPUtil/generic/FMA.h b/libc/src/__support/FPUtil/generic/FMA.h
index 0de4f98969b72..978583870629a 100644
--- a/libc/src/__support/FPUtil/generic/FMA.h
+++ b/libc/src/__support/FPUtil/generic/FMA.h
@@ -18,7 +18,7 @@
#include "src/__support/FPUtil/dyadic_float.h"
#include "src/__support/FPUtil/rounding_mode.h"
#include "src/__support/big_int.h"
-#include "src/__support/macros/attributes.h" // LIBC_INLINE
+#include "src/__support/macros/attributes.h" // LIBC_INLINE
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
@@ -275,10 +275,14 @@ fma(InType x, InType y, InType z) {
if (prod_mant == 0) {
// When there is exact cancellation, i.e., x*y == -z exactly, return -0.0 if
// rounding downward and +0.0 for other rounding modes.
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ prod_sign = Sign::POS;
+#else
if (fputil::quick_get_round() == FE_DOWNWARD)
prod_sign = Sign::NEG;
else
prod_sign = Sign::POS;
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
}
DyadicFloat result(prod_sign, prod_lsb_exp - InFPBits::EXP_BIAS, prod_mant);
diff --git a/libc/src/__support/FPUtil/generic/add_sub.h b/libc/src/__support/FPUtil/generic/add_sub.h
index 1dc16e4d4fce4..3142344e65abb 100644
--- a/libc/src/__support/FPUtil/generic/add_sub.h
+++ b/libc/src/__support/FPUtil/generic/add_sub.h
@@ -98,12 +98,16 @@ add_or_sub(InType x, InType y) {
if (y_bits.is_zero()) {
if (is_effectively_add)
return OutFPBits::zero(x_bits.sign()).get_val();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ return OutFPBits::zero(Sign::POS).get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_DOWNWARD:
return OutFPBits::zero(Sign::NEG).get_val();
default:
return OutFPBits::zero(Sign::POS).get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
}
if constexpr (cpp::is_same_v<InType, bfloat16> &&
@@ -136,12 +140,16 @@ add_or_sub(InType x, InType y) {
InType y_abs = y_bits.abs().get_val();
if (x_abs == y_abs && !is_effectively_add) {
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ return OutFPBits::zero(Sign::POS).get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_DOWNWARD:
return OutFPBits::zero(Sign::NEG).get_val();
default:
return OutFPBits::zero(Sign::POS).get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
}
Sign result_sign = Sign::POS;
diff --git a/libc/src/__support/FPUtil/generic/sqrt_80_bit_long_double.h b/libc/src/__support/FPUtil/generic/sqrt_80_bit_long_double.h
index cc692d0fc4c26..ad383e9eaf9d8 100644
--- a/libc/src/__support/FPUtil/generic/sqrt_80_bit_long_double.h
+++ b/libc/src/__support/FPUtil/generic/sqrt_80_bit_long_double.h
@@ -110,7 +110,11 @@ LIBC_INLINE LIBC_CONSTEXPR_DEFAULT long double sqrt(long double x) {
// Append the exponent field.
x_exp = ((x_exp >> 1) + LDBits::EXP_BIAS);
y |= (static_cast<StorageType>(x_exp) << (LDBits::FRACTION_LEN + 1));
-
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ // Round to nearest, ties to even
+ if (rb && (lsb || (r != 0)))
+ ++y;
+#else
switch (quick_get_round()) {
case FE_TONEAREST:
// Round to nearest, ties to even
@@ -122,6 +126,7 @@ LIBC_INLINE LIBC_CONSTEXPR_DEFAULT long double sqrt(long double x) {
++y;
break;
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
// Extract output
FPBits<long double> out(0.0L);
diff --git a/libc/src/__support/FPUtil/rounding_mode.h b/libc/src/__support/FPUtil/rounding_mode.h
index 3a4a41d8beacc..83ad911ed17f2 100644
--- a/libc/src/__support/FPUtil/rounding_mode.h
+++ b/libc/src/__support/FPUtil/rounding_mode.h
@@ -44,6 +44,9 @@ LIBC_INLINE bool fenv_is_round_down() {
// 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() {
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ return true;
+#endif
static volatile float x = 0x1.0p-24f;
float y = 1.5f + x;
return (y == 1.5f - x);
@@ -67,7 +70,6 @@ LIBC_INLINE bool fenv_is_round_to_zero() {
// Quick free standing get rounding mode based on the above observations.
LIBC_INLINE int quick_get_round() {
-#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
static volatile float x = 0x1.0p-24f;
float y = x;
float z = (0x1.000002p0f + y) + (-1.0f - y);
@@ -77,11 +79,6 @@ LIBC_INLINE int quick_get_round() {
if (z == 0x1.0p-23f)
return FE_TOWARDZERO;
return (2.0f + y == 2.0f) ? FE_TONEAREST : FE_UPWARD;
-#else
- // the optimization should've happened in each user of this function,
- // but just in case
- return FE_TONEAREST;
-#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
}
} // namespace generic
diff --git a/libc/src/__support/math/asinbf16.h b/libc/src/__support/math/asinbf16.h
index b291852762d73..b0edce0577d40 100644
--- a/libc/src/__support/math/asinbf16.h
+++ b/libc/src/__support/math/asinbf16.h
@@ -45,11 +45,13 @@ LIBC_INLINE LIBC_CONSTEXPR bfloat16 asinbf16(bfloat16 x) {
return x; // with sign
if (LIBC_UNLIKELY(x_abs <= 0x3D00)) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if ((xbits.is_pos() && rounding == FE_UPWARD) ||
(xbits.is_neg() && rounding == FE_DOWNWARD)) {
return fputil::cast<bfloat16>(fputil::multiply_add(xf, 0x1.0p-9f, xf));
}
+#endif
return x;
}
diff --git a/libc/src/__support/math/asinf16.h b/libc/src/__support/math/asinf16.h
index e9c9c6fca9184..86bedbdaba752 100644
--- a/libc/src/__support/math/asinf16.h
+++ b/libc/src/__support/math/asinf16.h
@@ -73,11 +73,13 @@ LIBC_INLINE constexpr float16 asinf16(float16 x) {
// else, in other rounding modes,
// asin(x) = x
if (LIBC_UNLIKELY(x_abs <= 0x1a1e)) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if ((xbits.is_pos() && rounding == FE_UPWARD) ||
(xbits.is_neg() && rounding == FE_DOWNWARD))
return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
+#endif
return x;
}
diff --git a/libc/src/__support/math/asinhf16.h b/libc/src/__support/math/asinhf16.h
index 39118d9357043..aaeec0dabb0f2 100644
--- a/libc/src/__support/math/asinhf16.h
+++ b/libc/src/__support/math/asinhf16.h
@@ -87,11 +87,13 @@ LIBC_INLINE constexpr float16 asinhf16(float16 x) {
// when |x| < 0x1.718p-5, asinhf16(x) = x. Adjust by 1 ULP for certain
// rounding types.
if (LIBC_UNLIKELY(x_abs < 0x29c6)) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if ((rounding == FE_UPWARD || rounding == FE_TOWARDZERO) && xf < 0)
return fputil::cast<float16>(xf + 0x1p-24f);
if ((rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) && xf > 0)
return fputil::cast<float16>(xf - 0x1p-24f);
+#endif
return fputil::cast<float16>(xf);
}
diff --git a/libc/src/__support/math/asinpi.h b/libc/src/__support/math/asinpi.h
index 2ee855d20455f..7979e7db9f939 100644
--- a/libc/src/__support/math/asinpi.h
+++ b/libc/src/__support/math/asinpi.h
@@ -84,6 +84,10 @@ LIBC_INLINE double asinpi(double x) {
MantT sticky_mask = (MantT(1) << (SHIFT_53 - 1)) - 1;
bool sticky = (r.mantissa & sticky_mask) != 0;
bool lsb = static_cast<bool>(m53 & 1);
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ // Carry if round_bit && (lsb || sticky) (round half to even).
+ raise_underflow = !(round_bit && (lsb || sticky));
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
// Carry if round_bit && (lsb || sticky) (round half to even).
@@ -100,6 +104,7 @@ LIBC_INLINE double asinpi(double x) {
raise_underflow = true; // truncation never carries
break;
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
}
}
if (raise_underflow)
diff --git a/libc/src/__support/math/coshf.h b/libc/src/__support/math/coshf.h
index d80d86af301ab..a993d29aae0bd 100644
--- a/libc/src/__support/math/coshf.h
+++ b/libc/src/__support/math/coshf.h
@@ -40,9 +40,11 @@ LIBC_INLINE float coshf(float x) {
if (xbits.is_inf_or_nan())
return x + FPBits::inf().get_val();
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if (LIBC_UNLIKELY(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO))
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
diff --git a/libc/src/__support/math/coshf16.h b/libc/src/__support/math/coshf16.h
index dc39944d89f85..ca3fda873329a 100644
--- a/libc/src/__support/math/coshf16.h
+++ b/libc/src/__support/math/coshf16.h
@@ -90,6 +90,11 @@ LIBC_INLINE constexpr float16 coshf16(float16 x) {
if (x_bits.is_inf())
return FPBits::inf().get_val();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+ return FPBits::inf().get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_UPWARD:
@@ -99,6 +104,7 @@ LIBC_INLINE constexpr float16 coshf16(float16 x) {
default:
return FPBits::max_normal().get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
}
}
diff --git a/libc/src/__support/math/exp.h b/libc/src/__support/math/exp.h
index 6c9f5a210a6a4..4a0247850a4c6 100644
--- a/libc/src/__support/math/exp.h
+++ b/libc/src/__support/math/exp.h
@@ -211,8 +211,10 @@ LIBC_INLINE double set_exceptional(double x) {
if (xbits.is_nan())
return x;
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
if (fputil::quick_get_round() == FE_UPWARD)
return FPBits::min_subnormal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_UNDERFLOW);
return 0.0;
@@ -221,9 +223,11 @@ LIBC_INLINE double set_exceptional(double x) {
// x >= round(log(MAX_NORMAL), D, RU) = 0x1.62e42fefa39fp+9 or +inf/nan
// x is finite
if (x_u < 0x7ff0'0000'0000'0000ULL) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
diff --git a/libc/src/__support/math/exp10.h b/libc/src/__support/math/exp10.h
index 6a8153cd90ff9..9dee8ef849db8 100644
--- a/libc/src/__support/math/exp10.h
+++ b/libc/src/__support/math/exp10.h
@@ -259,8 +259,10 @@ LIBC_INLINE double exp10_set_exceptional(double x) {
if (xbits.is_nan())
return x;
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
if (fputil::quick_get_round() == FE_UPWARD)
return FPBits::min_subnormal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_UNDERFLOW);
return 0.0;
@@ -272,9 +274,11 @@ LIBC_INLINE double exp10_set_exceptional(double x) {
// x >= log10(2^1024) or +inf/nan
// x is finite
if (x_u < 0x7ff0'0000'0000'0000ULL) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
diff --git a/libc/src/__support/math/exp10f.h b/libc/src/__support/math/exp10f.h
index 41f9a73a09935..d4761adf2962d 100644
--- a/libc/src/__support/math/exp10f.h
+++ b/libc/src/__support/math/exp10f.h
@@ -47,9 +47,11 @@ LIBC_INLINE float exp10f(float x) {
if (xbits.is_pos() && (x_u >= 0x421a'209bU)) {
// x is finite
if (x_u < 0x7f80'0000U) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
@@ -62,8 +64,12 @@ LIBC_INLINE float exp10f(float x) {
// When |x| <= log10(2)*2^-6
if (LIBC_UNLIKELY(x_abs <= 0x3b9a'209bU)) {
if (LIBC_UNLIKELY(x_u == 0xb25e'5bd9U)) { // x = -0x1.bcb7b2p-27f
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ return 0x1.fffffep-1f;
+#else
if (fputil::fenv_is_round_to_nearest())
return 0x1.fffffep-1f;
+#endif
}
// |x| < 2^-25
// 10^x ~ 1 + log(10) * x
diff --git a/libc/src/__support/math/exp10f16.h b/libc/src/__support/math/exp10f16.h
index 935a301544edf..2d0f20ed15cf4 100644
--- a/libc/src/__support/math/exp10f16.h
+++ b/libc/src/__support/math/exp10f16.h
@@ -82,6 +82,11 @@ LIBC_INLINE constexpr float16 exp10f16(float16 x) {
if (x_bits.is_inf())
return FPBits::inf().get_val();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ return FPBits::inf().get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_UPWARD:
@@ -91,6 +96,7 @@ LIBC_INLINE constexpr float16 exp10f16(float16 x) {
default:
return FPBits::max_normal().get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
// When x <= -8.
diff --git a/libc/src/__support/math/exp10m1f.h b/libc/src/__support/math/exp10m1f.h
index c45d3752d6f6c..cb8bb4fb05225 100644
--- a/libc/src/__support/math/exp10m1f.h
+++ b/libc/src/__support/math/exp10m1f.h
@@ -115,9 +115,11 @@ LIBC_INLINE float exp10m1f(float x) {
// When x >= log10(2^128), or x is nan
if (LIBC_UNLIKELY(xbits.is_pos() && x_u >= 0x421a'209bU)) {
if (xbits.is_finite()) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
@@ -164,10 +166,15 @@ LIBC_INLINE float exp10m1f(float x) {
if (xbits.is_nan())
return x;
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
+ if (x_u == 0xc0f0d2f1) // x = log10(2^-25)
+ return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f
+#else
int rounding = fputil::quick_get_round();
if (rounding == FE_UPWARD || rounding == FE_TOWARDZERO ||
(rounding == FE_TONEAREST && x_u == 0xc0f0d2f1))
return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_UNDERFLOW);
@@ -193,21 +200,27 @@ LIBC_INLINE float exp10m1f(float x) {
case 0x40e00000U: // x = 7.0f
return 9'999'999.0f;
case 0x41000000U: { // x = 8.0f
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_UPWARD || rounding == FE_TONEAREST)
return 100'000'000.0f;
+#endif
return 99'999'992.0f;
}
case 0x41100000U: { // x = 9.0f
int rounding = fputil::quick_get_round();
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
if (rounding == FE_UPWARD || rounding == FE_TONEAREST)
return 1'000'000'000.0f;
+#endif
return 999'999'936.0f;
}
case 0x41200000U: { // x = 10.0f
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_UPWARD || rounding == FE_TONEAREST)
return 10'000'000'000.0f;
+#endif
return 9'999'998'976.0f;
}
}
diff --git a/libc/src/__support/math/exp10m1f16.h b/libc/src/__support/math/exp10m1f16.h
index 7ee43a0092476..f9044f333676b 100644
--- a/libc/src/__support/math/exp10m1f16.h
+++ b/libc/src/__support/math/exp10m1f16.h
@@ -93,6 +93,11 @@ LIBC_INLINE constexpr float16 exp10m1f16(float16 x) {
if (x_bits.is_inf())
return FPBits::inf().get_val();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ return FPBits::inf().get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_UPWARD:
@@ -102,6 +107,7 @@ LIBC_INLINE constexpr float16 exp10m1f16(float16 x) {
default:
return FPBits::max_normal().get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
}
// When x < -11 * log10(2).
@@ -117,6 +123,9 @@ LIBC_INLINE constexpr float16 exp10m1f16(float16 x) {
}
// When x < -0x1.ce4p+1, round(10^x - 1, HP, RN) = -1.
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ return FPBits::one(Sign::NEG).get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_DOWNWARD:
@@ -124,6 +133,7 @@ LIBC_INLINE constexpr float16 exp10m1f16(float16 x) {
default:
return fputil::cast<float16>(-0x1.ffcp-1);
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
}
// When |x| <= 2^(-3).
@@ -150,8 +160,8 @@ LIBC_INLINE constexpr float16 exp10m1f16(float16 x) {
}
// When x is 1, 2, or 3. These are hard-to-round cases with exact results.
- // 10^4 - 1 = 9'999 is not exactly representable as a float16, but luckily the
- // polynomial approximation gives the correct result for x = 4 in all
+ // 10^4 - 1 = 9'999 is not exactly representable as a float16, but luckily
+ // the polynomial approximation gives the correct result for x = 4 in all
// rounding modes.
if (LIBC_UNLIKELY((x_u & ~(0x3c00U | 0x4000U | 0x4200U | 0x4400U)) == 0)) {
switch (x_u) {
diff --git a/libc/src/__support/math/exp2.h b/libc/src/__support/math/exp2.h
index 8f025bfa83227..eadcb423dcc13 100644
--- a/libc/src/__support/math/exp2.h
+++ b/libc/src/__support/math/exp2.h
@@ -238,8 +238,10 @@ LIBC_INLINE double set_exceptional(double x) {
if (xbits.is_nan())
return x;
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
if (fputil::quick_get_round() == FE_UPWARD)
return FPBits::min_subnormal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_UNDERFLOW);
return 0.0;
@@ -251,9 +253,11 @@ LIBC_INLINE double set_exceptional(double x) {
// x >= 1024 or +inf/nan
// x is finite
if (x_u < 0x7ff0'0000'0000'0000ULL) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
diff --git a/libc/src/__support/math/exp2f.h b/libc/src/__support/math/exp2f.h
index 439c752969b4d..b79341518c631 100644
--- a/libc/src/__support/math/exp2f.h
+++ b/libc/src/__support/math/exp2f.h
@@ -76,9 +76,11 @@ LIBC_INLINE float exp2f(float x) {
if (xbits.is_pos()) {
// x is finite
if (x_u < 0x7f80'0000U) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
diff --git a/libc/src/__support/math/exp2f16.h b/libc/src/__support/math/exp2f16.h
index 8ff918b665a6f..878d291c87a5d 100644
--- a/libc/src/__support/math/exp2f16.h
+++ b/libc/src/__support/math/exp2f16.h
@@ -66,6 +66,11 @@ LIBC_INLINE constexpr float16 exp2f16(float16 x) {
if (x_bits.is_inf())
return FPBits::inf().get_val();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ return FPBits::inf().get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_UPWARD:
@@ -75,6 +80,7 @@ LIBC_INLINE constexpr float16 exp2f16(float16 x) {
default:
return FPBits::max_normal().get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
// When x <= -25.
diff --git a/libc/src/__support/math/exp2m1f.h b/libc/src/__support/math/exp2m1f.h
index 49cd4589ace85..63b360ef70705 100644
--- a/libc/src/__support/math/exp2m1f.h
+++ b/libc/src/__support/math/exp2m1f.h
@@ -96,9 +96,11 @@ LIBC_INLINE float exp2m1f(float x) {
// x >= 128, or x is nan
if (xbits.is_pos()) {
if (xbits.is_finite()) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
@@ -117,9 +119,11 @@ LIBC_INLINE float exp2m1f(float x) {
if (xbits.is_nan())
return x;
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_UPWARD || rounding == FE_TOWARDZERO)
return -0x1.ffff'fep-1f; // -1.0f + 0x1.0p-24f
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_UNDERFLOW);
diff --git a/libc/src/__support/math/exp2m1f16.h b/libc/src/__support/math/exp2m1f16.h
index 196fc52158cfc..e1499735b0a3a 100644
--- a/libc/src/__support/math/exp2m1f16.h
+++ b/libc/src/__support/math/exp2m1f16.h
@@ -106,6 +106,11 @@ LIBC_INLINE constexpr float16 exp2m1f16(float16 x) {
if (x_bits.is_inf())
return FPBits::inf().get_val();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ return FPBits::inf().get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_UPWARD:
@@ -115,6 +120,7 @@ LIBC_INLINE constexpr float16 exp2m1f16(float16 x) {
default:
return FPBits::max_normal().get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
// When x < -11.
@@ -129,6 +135,9 @@ LIBC_INLINE constexpr float16 exp2m1f16(float16 x) {
fputil::cast<float16>(-0x1.ffcp-1));
// When x <= -12, round(2^x - 1, HP, RN) = -1.
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ return FPBits::one(Sign::NEG).get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_DOWNWARD:
@@ -136,6 +145,7 @@ LIBC_INLINE constexpr float16 exp2m1f16(float16 x) {
default:
return fputil::cast<float16>(-0x1.ffcp-1);
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
// When |x| <= 2^(-3).
diff --git a/libc/src/__support/math/expf.h b/libc/src/__support/math/expf.h
index a635c01d22772..711c0bc51eb4e 100644
--- a/libc/src/__support/math/expf.h
+++ b/libc/src/__support/math/expf.h
@@ -63,9 +63,11 @@ LIBC_INLINE float expf(float x) {
if (xbits.is_pos() && (xbits.uintval() >= 0x42b2'0000)) {
// x is finite
if (xbits.uintval() < 0x7f80'0000U) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
diff --git a/libc/src/__support/math/expf16.h b/libc/src/__support/math/expf16.h
index 2ca7c021e36b9..a6e848cd093ef 100644
--- a/libc/src/__support/math/expf16.h
+++ b/libc/src/__support/math/expf16.h
@@ -75,6 +75,11 @@ LIBC_INLINE constexpr float16 expf16(float16 x) {
if (x_bits.is_inf())
return FPBits::inf().get_val();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ return FPBits::inf().get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_UPWARD:
@@ -84,6 +89,7 @@ LIBC_INLINE constexpr float16 expf16(float16 x) {
default:
return FPBits::max_normal().get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
// When x <= -18.
@@ -95,12 +101,16 @@ LIBC_INLINE constexpr float16 expf16(float16 x) {
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
+ return FPBits::zero().get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_UPWARD:
return FPBits::min_subnormal().get_val();
default:
return FPBits::zero().get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
}
// When 0 < |x| <= 2^(-5).
diff --git a/libc/src/__support/math/expm1.h b/libc/src/__support/math/expm1.h
index 5902fc94578f4..4b82e08b27c51 100644
--- a/libc/src/__support/math/expm1.h
+++ b/libc/src/__support/math/expm1.h
@@ -265,9 +265,11 @@ LIBC_INLINE double set_exceptional(double x) {
// x >= round(log(MAX_NORMAL), D, RU) = 0x1.62e42fefa39fp+9 or +inf/nan
// x is finite
if (x_u < 0x7ff0'0000'0000'0000ULL) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
diff --git a/libc/src/__support/math/expm1f.h b/libc/src/__support/math/expm1f.h
index 60dd62bf7be8d..601fd133f804a 100644
--- a/libc/src/__support/math/expm1f.h
+++ b/libc/src/__support/math/expm1f.h
@@ -38,9 +38,13 @@ LIBC_INLINE float expm1f(float x) {
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Exceptional value
if (LIBC_UNLIKELY(x_u == 0x3e35'bec5U)) { // x = 0x1.6b7d8ap-3f
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+ return 0x1.8dbe64p-3f;
+#else
int round_mode = fputil::quick_get_round();
if (round_mode == FE_TONEAREST || round_mode == FE_UPWARD)
return 0x1.8dbe64p-3f;
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
return 0x1.8dbe62p-3f;
}
#if !defined(LIBC_TARGET_CPU_HAS_FMA_DOUBLE)
@@ -63,17 +67,21 @@ LIBC_INLINE float expm1f(float x) {
// exp(nan) = nan
if (xbits.is_nan())
return x;
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
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
+#endif
return -1.0f;
} else {
// x >= 89 or nan
if (xbits.uintval() >= 0x42b2'0000) {
if (xbits.uintval() < 0x7f80'0000U) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
int rounding = fputil::quick_get_round();
if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
return FPBits::max_normal().get_val();
+#endif
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
diff --git a/libc/src/__support/math/expm1f16.h b/libc/src/__support/math/expm1f16.h
index a14a6f051a834..edb8ce56e1a26 100644
--- a/libc/src/__support/math/expm1f16.h
+++ b/libc/src/__support/math/expm1f16.h
@@ -86,6 +86,11 @@ LIBC_INLINE constexpr float16 expm1f16(float16 x) {
if (x_bits.is_inf())
return FPBits::inf().get_val();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ return FPBits::inf().get_val();
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_UPWARD:
@@ -95,6 +100,7 @@ LIBC_INLINE constexpr float16 expm1f16(float16 x) {
default:
return FPBits::max_normal().get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
}
// When x <= -11 * log(2).
diff --git a/libc/src/__support/math/sin.h b/libc/src/__support/math/sin.h
index cd0c15abe62db..4b266b63e004c 100644
--- a/libc/src/__support/math/sin.h
+++ b/libc/src/__support/math/sin.h
@@ -54,11 +54,13 @@ LIBC_INLINE double sin(double x) {
return fputil::multiply_add(x, -0x1.0p-54, x);
#else
if (LIBC_UNLIKELY(x_e < 4)) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding_mode = fputil::quick_get_round();
if (rounding_mode == FE_TOWARDZERO ||
(xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) ||
(xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD))
return FPBits(xbits.uintval() - 1).get_val();
+#endif // !LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
return fputil::multiply_add(x, -0x1.0p-54, x);
#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
diff --git a/libc/src/__support/math/sincos.h b/libc/src/__support/math/sincos.h
index f7dbb0e0adcc6..5cf3fae72b9b8 100644
--- a/libc/src/__support/math/sincos.h
+++ b/libc/src/__support/math/sincos.h
@@ -67,11 +67,13 @@ LIBC_INLINE void sincos(double x, double *sin_x, double *cos_x) {
*cos_x = fputil::round_result_slightly_down(1.0);
if (LIBC_UNLIKELY(x_e < 4)) {
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding_mode = fputil::quick_get_round();
if (rounding_mode == FE_TOWARDZERO ||
(xbits.sign() == Sign::POS && rounding_mode == FE_DOWNWARD) ||
(xbits.sign() == Sign::NEG && rounding_mode == FE_UPWARD))
*sin_x = FPBits(xbits.uintval() - 1).get_val();
+#endif // !LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
*sin_x = fputil::multiply_add(x, -0x1.0p-54, x);
#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
diff --git a/libc/src/__support/math/sincosf.h b/libc/src/__support/math/sincosf.h
index a6b1cd8ae3010..64fc3e01f032a 100644
--- a/libc/src/__support/math/sincosf.h
+++ b/libc/src/__support/math/sincosf.h
@@ -176,6 +176,10 @@ LIBC_INLINE void sincosf(float x, float *sinp, float *cosp) {
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;
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+ s += EXCEPT_OUTPUTS_SIN[i][3];
+ c += EXCEPT_OUTPUTS_COS[i][3];
+#else
switch (fputil::quick_get_round()) {
case FE_UPWARD:
s += x_sign ? EXCEPT_OUTPUTS_SIN[i][2] : EXCEPT_OUTPUTS_SIN[i][1];
@@ -190,6 +194,7 @@ LIBC_INLINE void sincosf(float x, float *sinp, float *cosp) {
c += EXCEPT_OUTPUTS_COS[i][3];
break;
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
*sinp = x_sign ? -FPBits(s).get_val() : FPBits(s).get_val();
*cosp = FPBits(c).get_val();
diff --git a/libc/src/__support/math/sinf.h b/libc/src/__support/math/sinf.h
index c61beed749900..3641b306efdf9 100644
--- a/libc/src/__support/math/sinf.h
+++ b/libc/src/__support/math/sinf.h
@@ -152,10 +152,12 @@ LIBC_INLINE float sinf(float x) {
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
if (LIBC_UNLIKELY(x_abs == 0x4619'9998U)) { // x = 0x1.33333p13
float r = -0x1.63f4bap-2f;
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if ((rounding == FE_DOWNWARD && xbits.is_pos()) ||
(rounding == FE_UPWARD && xbits.is_neg()))
r = -0x1.63f4bcp-2f;
+#endif // !LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
return xbits.is_neg() ? -r : r;
}
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
diff --git a/libc/src/__support/math/sinf16.h b/libc/src/__support/math/sinf16.h
index ec27c4ec61c53..86972c7379371 100644
--- a/libc/src/__support/math/sinf16.h
+++ b/libc/src/__support/math/sinf16.h
@@ -79,7 +79,9 @@ LIBC_INLINE float16 sinf16(float16 x) {
return r.value();
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
+#endif
// Exhaustive tests show that for |x| <= 0x1.f4p-11, 1ULP rounding errors
// occur. To fix this, the following apply:
@@ -88,6 +90,7 @@ LIBC_INLINE float16 sinf16(float16 x) {
if (LIBC_UNLIKELY(x_abs == 0U))
return x;
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
// When x > 0, and rounding upward, sin(x) == x.
// When x < 0, and rounding downward, sin(x) == x.
if ((rounding == FE_UPWARD && xbits.is_pos()) ||
@@ -99,6 +102,10 @@ LIBC_INLINE float16 sinf16(float16 x) {
x_u--;
return FPBits(x_u).get_val();
}
+#endif // !LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+
+ // TODO: what about the case with rounding nearest and we're in here?
+ // Previously, it's an UB
}
if (xbits.is_inf_or_nan()) {
diff --git a/libc/src/__support/math/sinhf.h b/libc/src/__support/math/sinhf.h
index 141d87ba72bd2..258ddd290295f 100644
--- a/libc/src/__support/math/sinhf.h
+++ b/libc/src/__support/math/sinhf.h
@@ -61,6 +61,7 @@ LIBC_INLINE float sinhf(float x) {
if (xbits.is_inf())
return x;
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
if (xbits.is_neg()) {
if (LIBC_UNLIKELY(rounding == FE_UPWARD || rounding == FE_TOWARDZERO))
@@ -69,6 +70,7 @@ LIBC_INLINE float sinhf(float x) {
if (LIBC_UNLIKELY(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO))
return FPBits::max_normal().get_val();
}
+#endif // !LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_OVERFLOW);
diff --git a/libc/src/__support/math/sinhf16.h b/libc/src/__support/math/sinhf16.h
index 19ff5a3e12423..845fb7e736eec 100644
--- a/libc/src/__support/math/sinhf16.h
+++ b/libc/src/__support/math/sinhf16.h
@@ -129,6 +129,11 @@ LIBC_INLINE constexpr float16 sinhf16(float16 x) {
if (x_bits.is_inf())
return FPBits::inf(x_bits.sign()).get_val();
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
+ return FPBits::inf(x_bits.sign()).get_val();
+#else
int rounding_mode = fputil::quick_get_round();
if (rounding_mode == FE_TONEAREST ||
(x_bits.is_pos() && rounding_mode == FE_UPWARD) ||
@@ -137,6 +142,7 @@ LIBC_INLINE constexpr float16 sinhf16(float16 x) {
fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
return FPBits::inf(x_bits.sign()).get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_TO_NEAREST
return FPBits::max_normal(x_bits.sign()).get_val();
}
diff --git a/libc/src/__support/math/tan.h b/libc/src/__support/math/tan.h
index d1ebc035566f7..5f825ce10951a 100644
--- a/libc/src/__support/math/tan.h
+++ b/libc/src/__support/math/tan.h
@@ -149,10 +149,13 @@ LIBC_INLINE double tan(double x) {
return fputil::multiply_add(x, 0x1.0p-54, x);
#else
if (LIBC_UNLIKELY(x_e < 4)) {
+ // TODO: UB for rounding nearest
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding_mode = fputil::quick_get_round();
if ((xbits.sign() == Sign::POS && rounding_mode == FE_UPWARD) ||
(xbits.sign() == Sign::NEG && rounding_mode == FE_DOWNWARD))
return FPBits(xbits.uintval() + 1).get_val();
+#endif // !LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
return fputil::multiply_add(x, 0x1.0p-54, x);
#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
diff --git a/libc/src/__support/math/tanf16.h b/libc/src/__support/math/tanf16.h
index 6b9b9224fb84d..bbf5115dc0607 100644
--- a/libc/src/__support/math/tanf16.h
+++ b/libc/src/__support/math/tanf16.h
@@ -68,6 +68,7 @@ LIBC_INLINE float16 tanf16(float16 x) {
if (LIBC_UNLIKELY(x_abs == 0))
return x;
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
int rounding = fputil::quick_get_round();
// Exhaustive tests show that, when:
@@ -77,6 +78,7 @@ LIBC_INLINE float16 tanf16(float16 x) {
if ((xbits.is_pos() && rounding == FE_UPWARD) ||
(xbits.is_neg() && rounding == FE_DOWNWARD))
return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
+#endif // !LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
return x;
}
diff --git a/libc/src/__support/math/tanhf16.h b/libc/src/__support/math/tanhf16.h
index 80671691203c9..3119d4d850e8b 100644
--- a/libc/src/__support/math/tanhf16.h
+++ b/libc/src/__support/math/tanhf16.h
@@ -64,6 +64,9 @@ LIBC_INLINE float16 tanhf16(float16 x) {
// When -2^(-14) <= x <= -2^(-9).
if (x_u >= 0x8400U && x_u <= 0x9800U) {
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+ return x;
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
case FE_DOWNWARD:
@@ -71,6 +74,7 @@ LIBC_INLINE float16 tanhf16(float16 x) {
default:
return FPBits(static_cast<uint16_t>(x_u - 1U)).get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
// When |x| <= 0x1.d2p-4.
@@ -98,12 +102,18 @@ LIBC_INLINE float16 tanhf16(float16 x) {
// When |x| >= atanh(1 - 2^(-11)).
fputil::raise_except_if_required(FE_INEXACT);
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+ if (x_abs >= 0x4482U) {
+ return FPBits::one(x_bits.sign()).get_val();
+ }
+#else
int rounding_mode = fputil::quick_get_round();
if ((rounding_mode == FE_TONEAREST && x_abs >= 0x4482U) ||
(rounding_mode == FE_UPWARD && x_bits.is_pos()) ||
(rounding_mode == FE_DOWNWARD && x_bits.is_neg())) {
return FPBits::one(x_bits.sign()).get_val();
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
if (x_bits.is_pos())
return fputil::cast<float16>(0x1.ffcp-1);
return fputil::cast<float16>(-0x1.ffcp-1);
diff --git a/libc/src/__support/str_to_float.h b/libc/src/__support/str_to_float.h
index 830db2485fad8..25d049549d3f8 100644
--- a/libc/src/__support/str_to_float.h
+++ b/libc/src/__support/str_to_float.h
@@ -1146,6 +1146,9 @@ strtofloatingpoint(const CharType *__restrict src) {
}
RoundDirection round_direction = RoundDirection::Nearest;
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+ round_direction = RoundDirection::Nearest;
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
round_direction = RoundDirection::Nearest;
@@ -1160,6 +1163,7 @@ strtofloatingpoint(const CharType *__restrict src) {
round_direction = RoundDirection::Down;
break;
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
StrToNumResult<ExpandedFloat<T>> parse_result({0, 0});
if (base == 16) {
diff --git a/libc/src/stdio/printf_core/float_dec_converter.h b/libc/src/stdio/printf_core/float_dec_converter.h
index 8b99688cb993c..e2a5198752b07 100644
--- a/libc/src/stdio/printf_core/float_dec_converter.h
+++ b/libc/src/stdio/printf_core/float_dec_converter.h
@@ -49,6 +49,14 @@ constexpr uint32_t MAX_BLOCK = 999999999;
LIBC_INLINE RoundDirection get_round_direction(int last_digit, bool truncated,
Sign sign) {
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ if (last_digit != 5) {
+ return last_digit > 5 ? RoundDirection::Up : RoundDirection::Down;
+ } else {
+ return !truncated ? RoundDirection::Even : RoundDirection::Up;
+ }
+#endif
+
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
// Round to nearest, if it's exactly halfway then round to even.
diff --git a/libc/src/stdio/printf_core/float_hex_converter.h b/libc/src/stdio/printf_core/float_hex_converter.h
index 5c6899156ccba..660b067aeba15 100644
--- a/libc/src/stdio/printf_core/float_hex_converter.h
+++ b/libc/src/stdio/printf_core/float_hex_converter.h
@@ -112,6 +112,13 @@ LIBC_INLINE int convert_float_hex_exp(Writer<write_mode> *writer,
mantissa >>= shift_amount;
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ // 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);
+#else
switch (fputil::quick_get_round()) {
case FE_TONEAREST:
// Round to nearest, if it's exactly halfway then round to even.
@@ -131,6 +138,7 @@ LIBC_INLINE int convert_float_hex_exp(Writer<write_mode> *writer,
case FE_TOWARDZERO:
break;
}
+#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
// If the rounding caused an overflow, shift the mantissa and adjust the
// exponent to match.
>From 4ea936deebdf7f6839a99a0ccc287c9767bb5ae3 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Tue, 2 Jun 2026 23:35:59 +0700
Subject: [PATCH 3/4] updated to also handle `fenv_is_round_to_nearest()`'s
dependant
---
libc/src/__support/FPUtil/rounding_mode.h | 3 ---
libc/src/__support/math/sinhf.h | 4 ++++
2 files changed, 4 insertions(+), 3 deletions(-)
diff --git a/libc/src/__support/FPUtil/rounding_mode.h b/libc/src/__support/FPUtil/rounding_mode.h
index 83ad911ed17f2..6ce693d41da50 100644
--- a/libc/src/__support/FPUtil/rounding_mode.h
+++ b/libc/src/__support/FPUtil/rounding_mode.h
@@ -44,9 +44,6 @@ LIBC_INLINE bool fenv_is_round_down() {
// 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() {
-#ifdef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
- return true;
-#endif
static volatile float x = 0x1.0p-24f;
float y = 1.5f + x;
return (y == 1.5f - x);
diff --git a/libc/src/__support/math/sinhf.h b/libc/src/__support/math/sinhf.h
index 258ddd290295f..084b0d6586204 100644
--- a/libc/src/__support/math/sinhf.h
+++ b/libc/src/__support/math/sinhf.h
@@ -32,8 +32,12 @@ LIBC_INLINE float sinhf(float x) {
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// |x| = 0.0005589424981735646724700927734375
if (LIBC_UNLIKELY(x_abs == 0x3a12'85ffU)) {
+#ifdef LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
+ return x;
+#else
if (fputil::fenv_is_round_to_nearest())
return x;
+#endif // LIBC_MATH_HAS_ALWAYS_ROUNDING_NEAREST
}
#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
>From 45582522ae1a75a46c076fec64d293ea7c3e836f Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Wed, 3 Jun 2026 11:09:52 +0700
Subject: [PATCH 4/4] wrong guarding in ldexp
---
libc/src/__support/FPUtil/ManipulationFunctions.h | 11 ++++++-----
1 file changed, 6 insertions(+), 5 deletions(-)
diff --git a/libc/src/__support/FPUtil/ManipulationFunctions.h b/libc/src/__support/FPUtil/ManipulationFunctions.h
index a8c034cfaa118..6dc673daf45e5 100644
--- a/libc/src/__support/FPUtil/ManipulationFunctions.h
+++ b/libc/src/__support/FPUtil/ManipulationFunctions.h
@@ -171,15 +171,16 @@ ldexp(T x, U exp) {
// Make sure that we can safely cast exp to int when not returning early.
static_assert(EXP_LIMIT <= INT_MAX && -EXP_LIMIT >= INT_MIN);
-#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
if (LIBC_UNLIKELY(exp > EXP_LIMIT)) {
- int rounding_mode = quick_get_round();
Sign sign = bits.sign();
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ int rounding_mode = quick_get_round();
if ((sign == Sign::POS && rounding_mode == FE_DOWNWARD) ||
(sign == Sign::NEG && rounding_mode == FE_UPWARD) ||
(rounding_mode == FE_TOWARDZERO))
return FPBits<T>::max_normal(sign).get_val();
+#endif
set_errno_if_required(ERANGE);
raise_except_if_required(FE_OVERFLOW);
@@ -188,18 +189,18 @@ ldexp(T x, U exp) {
// 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();
-
+#ifndef LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
+ int rounding_mode = quick_get_round();
if ((sign == Sign::POS && rounding_mode == FE_UPWARD) ||
(sign == Sign::NEG && rounding_mode == FE_DOWNWARD))
return FPBits<T>::min_subnormal(sign).get_val();
+#endif
set_errno_if_required(ERANGE);
raise_except_if_required(FE_UNDERFLOW);
return FPBits<T>::zero(sign).get_val();
}
-#endif // LIBC_MATH_HAS_ALWAYS_ROUND_NEAREST
// For all other values, NormalFloat to T conversion handles it the right way.
DyadicFloat<FPBits<T>::STORAGE_LEN> normal(bits.get_val());
More information about the libc-commits
mailing list