[libc-commits] [libc] [libc][math] Correct expm1f16 subnormal rounding and range reporting (PR #236310)

[via libc-commits]

https://github.com/sriramshastry updated https://github.com/llvm/llvm-project/pull/236310

>From 087bd2de91990ccb2800d6aa6594481960def4fd Mon Sep 17 00:00:00 2001
From: sriramshastry <sriramshastry at gmail.com>
Date: Wed, 28 Oct 2026 09:16:36 +0530
Subject: [PATCH] [libc][math] Correct expm1f16 subnormal rounding and range
 reporting

For nonzero binary16 subnormal inputs, expm1(x)-x is positive and less than
half a subnormal ULP. Select the correctly rounded result from x and its
adjacent value, preserving signed zero and directed rounding.

Report underflow and ERANGE according to the rounded result, including the
transition at negative minimum normal. Add exhaustive signed-subnormal and
negative-minimum-normal rounding-mode regressions.

MPFR's mpfr_subnormalize underflow flag differs from this tininess-after-
rounding policy at 0x03ff under upward rounding; see the PR description and
the associated validation results.

Signed-off-by: sriramshastry <sriramshastry at gmail.com>
---
 libc/src/__support/math/expm1f16.h         | 20 +++++++
 libc/test/src/math/smoke/expm1f16_test.cpp | 61 ++++++++++++++++++++++
 2 files changed, 81 insertions(+)

diff --git a/libc/src/__support/math/expm1f16.h b/libc/src/__support/math/expm1f16.h
index 08234a0ad4ace9..072ead28957f71 100644
--- a/libc/src/__support/math/expm1f16.h
+++ b/libc/src/__support/math/expm1f16.h
@@ -120,6 +120,26 @@ LIBC_INLINE constexpr float16 expm1f16(float16 x) {
 
     // When 0 < |x| <= 2^(-3).
     if (x_abs <= 0x3000U && !x_bits.is_zero()) {
+      // The correction to x is nonzero but less than half a subnormal ulp.
+      // Select the adjacent half value directly, preserving signed zero.
+      if (LIBC_UNLIKELY(x_abs < 0x0400U || x_u == 0x8400U)) {
+        uint16_t out = x_u;
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+        int mode = fputil::quick_get_round();
+        bool increase =
+            (x_bits.is_pos() && mode == FE_UPWARD) ||
+            (x_bits.is_neg() && (mode == FE_UPWARD || mode == FE_TOWARDZERO));
+        if (increase)
+          out = static_cast<uint16_t>(x_u + (x_bits.is_neg() ? -1 : 1));
+#endif
+        int excepts = FE_INEXACT;
+        if ((out & 0x7fffU) < 0x0400U) {
+          fputil::set_errno_if_required(ERANGE);
+          excepts |= FE_UNDERFLOW;
+        }
+        fputil::raise_except_if_required(excepts);
+        return FPBits(out).get_val();
+      }
 
 #ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
       if (auto r = EXPM1F16_EXCEPTS_LO.lookup(x_u);
diff --git a/libc/test/src/math/smoke/expm1f16_test.cpp b/libc/test/src/math/smoke/expm1f16_test.cpp
index 545304ced63ed8..2a136bf6b22282 100644
--- a/libc/test/src/math/smoke/expm1f16_test.cpp
+++ b/libc/test/src/math/smoke/expm1f16_test.cpp
@@ -15,6 +15,67 @@
 
 using LlvmLibcExpm1f16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
 
+TEST_F(LlvmLibcExpm1f16Test, AllSubnormalRounding) {
+  using Bits = LIBC_NAMESPACE::fputil::FPBits<float16>;
+  using namespace LIBC_NAMESPACE::fputil::testing;
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+  constexpr RoundingMode MODES[] = {RoundingMode::Nearest};
+#else
+  constexpr RoundingMode MODES[] = {RoundingMode::Nearest, RoundingMode::Upward,
+                                    RoundingMode::Downward,
+                                    RoundingMode::TowardZero};
+#endif
+  for (RoundingMode mode : MODES) {
+    ForceRoundingMode rounding(mode);
+    if (!rounding.success)
+      continue;
+    for (uint16_t magnitude = 1; magnitude < 0x0400; ++magnitude) {
+      constexpr uint16_t SIGNS[] = {0, 0x8000};
+      for (uint16_t sign : SIGNS) {
+        bool negative = sign != 0;
+        uint16_t expected = static_cast<uint16_t>(magnitude | sign);
+        if ((!negative && mode == RoundingMode::Upward) ||
+            (negative && (mode == RoundingMode::Upward ||
+                          mode == RoundingMode::TowardZero)))
+          expected = static_cast<uint16_t>(expected + (negative ? -1 : 1));
+        bool tiny = (expected & 0x7fff) < 0x0400;
+        EXPECT_FP_EQ_WITH_EXCEPTION(
+            Bits(expected).get_val(),
+            LIBC_NAMESPACE::expm1f16(
+                Bits(uint16_t(magnitude | sign)).get_val()),
+            FE_INEXACT | (tiny ? FE_UNDERFLOW : 0));
+        EXPECT_MATH_ERRNO(tiny ? ERANGE : 0);
+      }
+    }
+  }
+}
+
+TEST_F(LlvmLibcExpm1f16Test, NegativeMinNormalTininessAfterRounding) {
+  using Bits = LIBC_NAMESPACE::fputil::FPBits<float16>;
+  using namespace LIBC_NAMESPACE::fputil::testing;
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+  constexpr RoundingMode MODES[] = {RoundingMode::Nearest};
+#else
+  constexpr RoundingMode MODES[] = {RoundingMode::Nearest, RoundingMode::Upward,
+                                    RoundingMode::Downward,
+                                    RoundingMode::TowardZero};
+#endif
+  for (RoundingMode mode : MODES) {
+    ForceRoundingMode rounding(mode);
+    if (!rounding.success)
+      continue;
+    uint16_t expected = 0x8400;
+    if (mode == RoundingMode::Upward || mode == RoundingMode::TowardZero)
+      expected = 0x83ff;
+    bool tiny = (expected & 0x7fff) < 0x0400;
+    EXPECT_FP_EQ_WITH_EXCEPTION(
+        Bits(expected).get_val(),
+        LIBC_NAMESPACE::expm1f16(Bits(uint16_t(0x8400)).get_val()),
+        FE_INEXACT | (tiny ? FE_UNDERFLOW : 0));
+    EXPECT_MATH_ERRNO(tiny ? ERANGE : 0);
+  }
+}
+
 TEST_F(LlvmLibcExpm1f16Test, SpecialNumbers) {
   EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::expm1f16(aNaN));
   EXPECT_MATH_ERRNO(0);