[libc-commits] [libc] [libc][math] Short-circuit certain powf range errors (PR #223719)

[via libc-commits]

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

>From 28685b2da514cbd7a788c3cc2482c2f80b0a91cc Mon Sep 17 00:00:00 2001
From: sriramshastry <sriramshastry at gmail.com>
Date: Fri, 9 Oct 2026 16:21:13 +0000
Subject: [PATCH] [libc][math] Short-circuit provable powf subnormal underflow

For a nonzero subnormal base, reuse the normalization shift to bound
its magnitude below 2^(-125-shift). For positive y, y*(125+shift) >= 150
proves the result is strictly below half the minimum binary32 subnormal.
Return through the existing signed underflow helper after special-case
and domain handling.

Use an exact double product when hardware double is available, and an
exact unsigned significand comparison otherwise. Add rounding-mode and
threshold-neighbor smoke coverage for both evaluators.

The redesigned implementation already handles extreme exponents, so do
not retain the old large-exponent shortcut. Performance remains workload
and layout dependent; this is not a uniform no-regression claim.

Assisted-by: GitHub Copilot
Signed-off-by: sriramshastry <sriramshastry at gmail.com>
---
 libc/src/__support/math/powf_utils.h   | 13 +++++++++
 libc/test/src/math/smoke/powf_test.cpp | 40 ++++++++++++++++++++++++++
 2 files changed, 53 insertions(+)

diff --git a/libc/src/__support/math/powf_utils.h b/libc/src/__support/math/powf_utils.h
index 2a2253bd273796..ac74b851d44ea7 100644
--- a/libc/src/__support/math/powf_utils.h
+++ b/libc/src/__support/math/powf_utils.h
@@ -598,6 +598,19 @@ check_exceptional_cases(float &x, float &y, int &ex, SignType &sign) {
   // Normalize denormal inputs.
   if (x_a < FloatBits::min_normal().uintval()) {
     int shift = cpp::countl_zero(x_a) - 8;
+#ifdef LIBC_TARGET_CPU_HAS_FPU_DOUBLE
+    if (!y_sign &&
+        static_cast<double>(y) * static_cast<double>(125 + shift) >= 150.0)
+      return set_underflow(sign == 0 ? Sign::POS : Sign::NEG);
+#else
+    uint32_t log2_bound = static_cast<uint32_t>(125 + shift);
+    if (!y_sign &&
+        (y_a >= 0x4000'0000U || (y_a >= FloatBits::one().uintval() &&
+                                 (y_a - 0x3f00'0000U) * log2_bound >=
+                                     (150U << FloatBits::FRACTION_LEN))))
+      return set_underflow(sign == 0 ? Sign::POS : Sign::NEG);
+#endif
+
     ex -= shift;
     x = cpp::bit_cast<float>(x_a << shift);
   }
diff --git a/libc/test/src/math/smoke/powf_test.cpp b/libc/test/src/math/smoke/powf_test.cpp
index 46d4f1b7120944..0b07abe78aebba 100644
--- a/libc/test/src/math/smoke/powf_test.cpp
+++ b/libc/test/src/math/smoke/powf_test.cpp
@@ -247,6 +247,43 @@ class PowfTest : public LIBC_NAMESPACE::testing::FPTest<float> {
     }
   }
 
+  void test_range_underflow(float (*func)(float, float)) {
+    constexpr float BASES[] = {0x1.8p-2f, 0x1.fffffep-2f, 0x1.8p-130f};
+    constexpr float EXPONENTS[] = {150.0f, 150.0f, 1.25f};
+    for (int mode = 0; mode < N_ROUNDING_MODES; ++mode) {
+      ForceRoundingMode rounding(ROUNDING_MODES[mode]);
+      if (!rounding.success)
+        continue;
+      float positive =
+          ROUNDING_MODES[mode] == RoundingMode::Upward ? min_denormal : zero;
+      float negative = ROUNDING_MODES[mode] == RoundingMode::Downward
+                           ? -min_denormal
+                           : neg_zero;
+      for (unsigned index = 0; index < 3; ++index)
+        EXPECT_FP_EQ_WITH_EXCEPTION(positive,
+                                    func(BASES[index], EXPONENTS[index]),
+                                    FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_FP_EQ_WITH_EXCEPTION(negative, func(-0x1.8p-2f, 151.0f),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_FP_EQ_WITH_EXCEPTION(negative, func(-0x1.8p-130f, 3.0f),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      float boundary = ROUNDING_MODES[mode] == RoundingMode::Upward ||
+                               ROUNDING_MODES[mode] == RoundingMode::Nearest
+                           ? min_denormal
+                           : zero;
+      EXPECT_FP_EQ_WITH_EXCEPTION(boundary, func(0x1.fffffep-2f, 0x1.2bfffep7f),
+                                  FE_UNDERFLOW | FE_INEXACT);
+      constexpr float SUBNORMAL_BOUND = 150.0f / 126.0f;
+      uint32_t bound_bits = FPBits(SUBNORMAL_BOUND).uintval();
+      EXPECT_FP_EQ_WITH_EXCEPTION(
+          boundary, func(0x1.fffffcp-127f, FPBits(bound_bits - 1).get_val()),
+          FE_UNDERFLOW | FE_INEXACT);
+      EXPECT_FP_EQ_WITH_EXCEPTION(
+          positive, func(0x1.fffffcp-127f, FPBits(bound_bits + 1).get_val()),
+          FE_UNDERFLOW | FE_INEXACT);
+    }
+  }
+
 #ifdef LIBC_TEST_FTZ_DAZ
   void test_ftz(float (*func)(float, float)) {
     LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::FTZ);
@@ -298,6 +335,9 @@ class PowfTest : public LIBC_NAMESPACE::testing::FPTest<float> {
   TEST_F(LlvmLibcPowfTest##suffix, SubnormalBase) {                            \
     test_subnormal_base(&func, tolerance);                                     \
   }                                                                            \
+  TEST_F(LlvmLibcPowfTest##suffix, RangeUnderflow) {                           \
+    test_range_underflow(&func);                                               \
+  }                                                                            \
   LIST_POWF_FTZ_DAZ_TESTS(suffix, func)                                        \
   static_assert(true, "Require semicolon.")