[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 ad8355f54a0e25fbc8537f7ed44b326bfa888ec6 Mon Sep 17 00:00:00 2001
From: Sriram Shastry <sriramshastry at gmail.com>
Date: Tue, 15 Sep 2026 11:28:12 +0530
Subject: [PATCH] [libc][math] Short-circuit certain powf range errors

Return overflow and underflow results directly when float bounds
prove that the exact powf result is outside the representable range.
Preserve rounding-mode results, errno, and floating-point exceptions
while avoiding the log2 and exp2 approximations for these inputs.
---
 libc/src/__support/math/CMakeLists.txt |  1 +
 libc/src/__support/math/powf.h         | 36 ++++++++++++++++++++++++++
 2 files changed, 37 insertions(+)

diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 0871a5bf229b6..d774cc7c26936 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -5301,6 +5301,7 @@ add_header_library(
     libc.src.__support.FPUtil.multiply_add
     libc.src.__support.FPUtil.nearest_integer
     libc.src.__support.FPUtil.polyeval
+    libc.src.__support.FPUtil.rounding_mode
     libc.src.__support.FPUtil.sqrt
     libc.src.__support.FPUtil.triple_double
     libc.src.__support.macros.optimization
diff --git a/libc/src/__support/math/powf.h b/libc/src/__support/math/powf.h
index 944b477ac4461..c66f36dd8a837 100644
--- a/libc/src/__support/math/powf.h
+++ b/libc/src/__support/math/powf.h
@@ -37,6 +37,7 @@
 #include "src/__support/FPUtil/double_double.h"
 #include "src/__support/FPUtil/multiply_add.h"
 #include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/FPUtil/rounding_mode.h"
 #include "src/__support/FPUtil/sqrt.h" // Speedup for powf(x, 1/2) = sqrtf(x)
 #include "src/__support/FPUtil/triple_double.h"
 #include "src/__support/common.h"
@@ -693,6 +694,26 @@ LIBC_INLINE float powf(float x, float y) {
       return FloatBits::quiet_nan().get_val();
     }
 
+    if (y_abs > 0x4f17'0000U && y_abs < 0x7f80'0000U && x_u != 0 &&
+        x_u < 0x7f80'0000U && x_u != 0x3f80'0000U && x_u != 0x4000'0000U &&
+        x_u != 0x4120'0000U) {
+      const bool overflow = (x_u > 0x3f80'0000U) != (y_u > 0x8000'0000U);
+      if (overflow) {
+        fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
+        static volatile float rounding_test = 0x1.0p-25f;
+        if (1.0f - rounding_test != 1.0f)
+          return FloatBits::max_normal().get_val();
+        fputil::set_errno_if_required(ERANGE);
+        return FloatBits::inf().get_val();
+      }
+
+      fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
+      if (fputil::fenv_is_round_up())
+        return FloatBits::min_subnormal().get_val();
+      fputil::set_errno_if_required(ERANGE);
+      return 0.0f;
+    }
+
     // Exceptional exponents.
     if (y == 0.0f)
       return 1.0f;
@@ -781,6 +802,21 @@ LIBC_INLINE float powf(float x, float y) {
     }
   }
 
+  if (LIBC_UNLIKELY(y_u >= 0x3f80'0000U && y_u < 0x7f80'0000U && x_u > 0 &&
+                    x_u < 0x3f00'0000U)) {
+    // If x is in [2^(e - 127), 2^(e - 126)), then
+    // log2(x) < -(126 - e).  Use this bound to identify results below
+    // 2^-150 without evaluating log2(x).
+    int log2_bound = 126 - static_cast<int>(x_u >> FloatBits::FRACTION_LEN);
+    if (LIBC_UNLIKELY(static_cast<double>(y) * log2_bound >= 150.0)) {
+      fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
+      if (fputil::fenv_is_round_up())
+        return FloatBits::min_subnormal().get_val();
+      fputil::set_errno_if_required(ERANGE);
+      return 0.0f;
+    }
+  }
+
   int ex = -FloatBits::EXP_BIAS;
   uint64_t sign = 0;