[libc-commits] [libc] [libc][math] Fix code size issues for powf_float_eval. (PR #232601)
via libc-commits
libc-commits at lists.llvm.org
Thu Oct 15 22:57:30 PDT 2026
https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/232601
>From 727e6d0744e37e37159ed3b76cf8763eb57b0cd0 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 16 Oct 2026 05:38:55 +0000
Subject: [PATCH 1/2] [libc][math] Fix code size issues for powf_float_eval.
- Remove double and int64_t usage in powf_float_eval.
- Use __builtin_roundeven* for ARMv8.
- Reduce tables to 16-entry.
- Add more support/guards for LIBC_MATH_HAS_NO_EXCEPT and LIBC_ATH_HAS_ASSUME_ROUND_NEAREST_ONLY.
- Simplify check_exceptional cases.
Assisted-by: Gemini for analysis.
---
libc/src/__support/FPUtil/FEnvImpl.h | 2 +
libc/src/__support/FPUtil/nearest_integer.h | 19 ++-
libc/src/__support/math/powf_float_eval.h | 76 ++++------
libc/src/__support/math/powf_utils.h | 147 ++++++++++----------
libc/test/src/math/smoke/powf_test.cpp | 6 +-
5 files changed, 126 insertions(+), 124 deletions(-)
diff --git a/libc/src/__support/FPUtil/FEnvImpl.h b/libc/src/__support/FPUtil/FEnvImpl.h
index 2f99948fffc688..d556a3cbf4d39f 100644
--- a/libc/src/__support/FPUtil/FEnvImpl.h
+++ b/libc/src/__support/FPUtil/FEnvImpl.h
@@ -220,6 +220,7 @@ LIBC_INLINE LIBC_CONSTEXPR_DEFAULT void raise_underflow_except_if_required() {
raise_except_if_required(FE_UNDERFLOW | FE_INEXACT);
}
+#ifndef LIBC_MATH_HAS_NO_EXCEPT
#if defined(LIBC_TARGET_CPU_HAS_FPU_FLOAT)
template <>
LIBC_INLINE LIBC_CONSTEXPR_DEFAULT void
@@ -267,6 +268,7 @@ raise_underflow_except_if_required<double>() {
}
}
#endif // LIBC_TARGET_CPU_HAS_FPU_DOUBLE
+#endif // !LIBC_MATH_HAS_NO_EXCEPT
} // namespace fputil
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/__support/FPUtil/nearest_integer.h b/libc/src/__support/FPUtil/nearest_integer.h
index fabd309bd35202..f131181fd2c61a 100644
--- a/libc/src/__support/FPUtil/nearest_integer.h
+++ b/libc/src/__support/FPUtil/nearest_integer.h
@@ -9,6 +9,7 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_NEAREST_INTEGER_H
#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_NEAREST_INTEGER_H
+#include "src/__support/CPP/type_traits.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
#include "src/__support/macros/properties/architectures.h"
@@ -40,9 +41,15 @@ namespace fputil {
// Notice that for AARCH64 and x86-64 with SSE4.2 support, we will use their
// corresponding rounding instruction instead. And in those cases, the results
// are rounded to the nearest integer, tie-to-even.
-LIBC_INLINE static constexpr float nearest_integer(float x) {
+LIBC_INLINE constexpr float nearest_integer(float x) {
+#if __has_builtin(__builtin_roundevenf) && defined(__ARM_ARCH) && \
+ (__ARM_ARCH >= 8) && defined(LIBC_TARGET_CPU_HAS_ARM_FPU_FLOAT)
+ if (!cpp::is_constant_evaluated())
+ return __builtin_roundevenf(x);
+#endif
if (x < 0x1p24f && x > -0x1p24f) {
float r = x < 0 ? (x - 0x1.0p23f) + 0x1.0p23f : (x + 0x1.0p23f) - 0x1.0p23f;
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
float diff = x - r;
// The expression above is correct for the default rounding mode, round-to-
// nearest, tie-to-even. For other rounding modes, it might be off by 1,
@@ -51,14 +58,21 @@ LIBC_INLINE static constexpr float nearest_integer(float x) {
return r + 1.0f;
if (LIBC_UNLIKELY(diff < -0.5f))
return r - 1.0f;
+#endif // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
return r;
}
return x;
}
-LIBC_INLINE static constexpr double nearest_integer(double x) {
+LIBC_INLINE constexpr double nearest_integer(double x) {
+#if __has_builtin(__builtin_roundeven) && defined(__ARM_ARCH) && \
+ (__ARM_ARCH >= 8) && defined(LIBC_TARGET_CPU_HAS_ARM_FPU_DOUBLE)
+ if (!cpp::is_constant_evaluated())
+ return __builtin_roundeven(x);
+#endif
if (x < 0x1p53 && x > -0x1p53) {
double r = x < 0 ? (x - 0x1.0p52) + 0x1.0p52 : (x + 0x1.0p52) - 0x1.0p52;
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
double diff = x - r;
// The expression above is correct for the default rounding mode, round-to-
// nearest, tie-to-even. For other rounding modes, it might be off by 1,
@@ -67,6 +81,7 @@ LIBC_INLINE static constexpr double nearest_integer(double x) {
return r + 1.0;
if (LIBC_UNLIKELY(diff < -0.5))
return r - 1.0;
+#endif // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
return r;
}
return x;
diff --git a/libc/src/__support/math/powf_float_eval.h b/libc/src/__support/math/powf_float_eval.h
index 023ba56d55ccca..ae36b787c5a80a 100644
--- a/libc/src/__support/math/powf_float_eval.h
+++ b/libc/src/__support/math/powf_float_eval.h
@@ -66,11 +66,11 @@ namespace powf_internal = LIBC_NAMESPACE::math::powf_internal;
// relative error(log2(x)) <= 2^-25 / (2^8 * ln(2)) < 2^-33.
//
// To compute log2(x), we perform range reduction for log2(m_x) using a compact
-// 5-bit lookup table (32 entries):
-// dx = r * m_x - 1 in [-2^-6, 2^-5].
+// 4-bit lookup table (16 entries):
+// dx = r * m_x - 1 in [-2^-5, 2^-4].
// Then m_x = (1 + dx) / r, and:
// log2(m_x) = log2(1 + dx) - log2(r),
-// where -log2(r) is obtained from the lookup table LOG2_R_FF_32 with
+// where -log2(r) is obtained from the lookup table LOG2_R_FF_16 with
// FloatFloat precision (~48 bits). The computation of dx = r * m_x - 1 is
// exact.
//
@@ -84,14 +84,14 @@ namespace powf_internal = LIBC_NAMESPACE::math::powf_internal;
// Generated by Sollya with:
// > log2_e = 0x1.715476p0 + 0x1.4ae0cp-26;
// > f = (log2(1 + x) - x * log2_e) / x^2;
-// > Q = fpminimax(f, 4, [|SG...|], [-2^-6, 2^-5]);
-// > dirtyinfnorm(((log2(1 + x) - x * log2_e) - x^2 * Q) / x, [-2^-6, 2^-5]);
-// 0x1.0af1...p-34
+// > Q = fpminimax(f, 4, [|SG...|], [-2^-5, 2^-4]);
+// > dirtyinfnorm(((log2(1 + x) - x * log2_e) - x^2 * Q) / x, [-2^-5, 2^-4]);
+// 0x1.692c...p-34
//
// Finally, we compute the product y * log2(x) in FloatFloat:
// y * log2(x) = y * ((e_x - log2(r)) + log2(1 + dx)).
-// When x in [1 - 2^-5, 1) (where e_x = -1 and R[31] = 0.5 so
-// -log2(R[31]) = 1.0f), e_x - log2(r) cancels to 0, avoiding cancellation.
+// When x in [1 - 2^-4, 1) (where e_x = -1 and R[15] = 0.5 so
+// -log2(R[15]) = 1.0f), e_x - log2(r) cancels to 0, avoiding cancellation.
// We decompose into integer k and fractional remainder rem:
// k = round(total_hi),
// rem = (total_hi - k) + total_lo,
@@ -105,15 +105,8 @@ LIBC_INLINE float powf(float x, float y) {
FPBits xbits(x), ybits(y);
uint32_t x_u = xbits.uintval();
- uint32_t y_u = ybits.uintval();
uint32_t y_a = ybits.abs().uintval();
- if (LIBC_UNLIKELY((x_u & 0x007F'FFFF) == 0 || (y_u & 0x007F'FFFF) == 0)) {
- if (auto r = powf_internal::check_special_inputs(x, y);
- LIBC_UNLIKELY(r.has_value()))
- return r.value();
- }
-
int e_x = 0;
uint32_t sign = 0;
if (LIBC_UNLIKELY(y_a <= powf_internal::Y_LOWER_BOUND ||
@@ -131,7 +124,7 @@ LIBC_INLINE float powf(float x, float y) {
xbits = FPBits(x);
e_x += xbits.get_biased_exponent() - FPBits::EXP_BIAS;
uint32_t x_mant = xbits.get_mantissa();
- unsigned idx_x = static_cast<unsigned>(x_mant >> (FPBits::FRACTION_LEN - 5));
+ unsigned idx_x = static_cast<unsigned>(x_mant >> (FPBits::FRACTION_LEN - 4));
FPBits m_x = FPBits(x_mant | 0x3f80'0000U);
// log2(e) in FloatFloat:
@@ -142,16 +135,16 @@ LIBC_INLINE float powf(float x, float y) {
// Generated by Sollya with:
// > log2_e = 0x1.715476p0 + 0x1.4ae0cp-26;
// > f = (log2(1 + x) - x * log2_e) / x^2;
- // > Q = fpminimax(f, 4, [|SG...|], [-2^-6, 2^-5]);
- // > dirtyinfnorm(((log2(1 + x) - x * log2_e) - x^2 * Q) / x, [-2^-6, 2^-5]);
- // 0x1.0af1...p-34
- constexpr float COEFFS[] = {-0x1.715476p-1f, 0x1.ec709ep-2f, -0x1.71618ap-2f,
- 0x1.2810e4p-2f, -0x1.ae42b2p-3f};
+ // > Q = fpminimax(f, 4, [|SG...|], [-2^-5, 2^-4]);
+ // > dirtyinfnorm(((log2(1 + x) - x * log2_e) - x^2 * Q) / x, [-2^-5, 2^-4]);
+ // 0x1.692c...p-34
+ constexpr float COEFFS[] = {-0x1.715476p-1f, 0x1.ec70ap-2f, -0x1.715b5p-2f,
+ 0x1.279932p-2f, -0x1.cbc178p-3f};
// Range reduction for log2(m_x):
- // dx = r * m_x - 1 in [-2^-6, 2^-5).
+ // dx = r * m_x - 1 in [-2^-5, 2^-4).
// The computation is exact: m_x = (1 + dx) / r.
- float r = powf_internal::R_32[idx_x];
+ float r = powf_internal::R_16[idx_x];
#if defined(LIBC_TARGET_CPU_HAS_FMA_FLOAT)
float dx = fputil::multiply_add(r, m_x.get_val(), -1.0f); // Exact
#else // !LIBC_TARGET_CPU_HAS_FMA_FLOAT
@@ -163,7 +156,7 @@ LIBC_INLINE float powf(float x, float y) {
float dx = fputil::multiply_add(r, cd, dx_hi); // Exact
#endif // LIBC_TARGET_CPU_HAS_FMA_FLOAT
- fputil::FloatFloat log2_r = powf_internal::LOG2_R_FF_32[idx_x];
+ fputil::FloatFloat log2_r = powf_internal::LOG2_R_FF_16[idx_x];
float e_xf = static_cast<float>(e_x);
// (e_xf + log2_r.hi) is exact since log2_r.hi has 16 bits of precision and
// |e_x| <= 150 < 2^8.
@@ -171,8 +164,8 @@ LIBC_INLINE float powf(float x, float y) {
// Early overflow and underflow check:
// z = y * log2(x) ~ y * (e_x_hi + dx * log2(e))
- // With |dx| < 2^-5, |log2(1 + dx) - dx * log2(e)| <= (dx^2 / 2) * log2(e) <
- // 2^-10. For |y * log2(x)| <= 128, the error |z - y * log2(x)| < 4. Hence:
+ // With |dx| < 2^-4, |log2(1 + dx) - dx * log2(e)| <= (dx^2 / 2) * log2(e) <
+ // 2^-8. For |y * log2(x)| <= 128, the error |z - y * log2(x)| < 4. Hence:
// z > 132 => overflow
// z < -155 => underflow
float z = y * fputil::multiply_add(dx, LOG2_E.hi, e_x_hi);
@@ -213,43 +206,26 @@ LIBC_INLINE float powf(float x, float y) {
// k = round(total_hi)
// rem = (total_hi - k) + total_lo, with |rem| <= 0.5.
float kf = fputil::nearest_integer(total_hi);
- int64_t k = static_cast<int64_t>(kf);
+ int k = static_cast<int>(kf);
float rem_hi = total_hi - kf;
float rem = rem_hi + total_lo;
- // Normal range: -125 <= k <= 127.
- uint64_t k_u = static_cast<uint64_t>(k + 125);
- if (LIBC_LIKELY(k_u <= 252)) {
- float result = exp2f_eval(rem, static_cast<int>(k));
- return cpp::bit_cast<float>(cpp::bit_cast<uint32_t>(result) | sign);
- }
-
- // Exceptional path: k > 128 (overflow), k < -150 (underflow), or subnormal.
- if (k > 128)
+ if (LIBC_UNLIKELY(k > 128))
return set_overflow(out_sign);
- if (k < -150)
+ if (LIBC_UNLIKELY(k < -150))
return set_underflow(out_sign);
- // Subnormal range: -150 <= k <= -126.
- if (k <= -126) {
- float result = exp2f_eval(rem, static_cast<int>(k));
- float res = cpp::bit_cast<float>(cpp::bit_cast<uint32_t>(result) | sign);
-
- if (LIBC_UNLIKELY(FPBits(res).is_normal()))
- return res;
+ float result = exp2f_eval(rem, k);
+ if (LIBC_UNLIKELY(k == 128 && FPBits(result).is_inf()))
+ return set_overflow(out_sign);
+ if (LIBC_UNLIKELY(k <= -126 && !FPBits(result).is_normal())) {
fputil::set_errno_if_required(ERANGE);
fputil::raise_underflow_except_if_required<float>();
- return res;
}
- // Boundary case: k = 128.
- float result = exp2f_eval(rem, static_cast<int>(k));
- if (LIBC_UNLIKELY(FPBits(result).is_inf()))
- return set_overflow(out_sign);
-
return cpp::bit_cast<float>(cpp::bit_cast<uint32_t>(result) | sign);
}
diff --git a/libc/src/__support/math/powf_utils.h b/libc/src/__support/math/powf_utils.h
index 2a2253bd273796..cb0efe6bcacd90 100644
--- a/libc/src/__support/math/powf_utils.h
+++ b/libc/src/__support/math/powf_utils.h
@@ -36,10 +36,12 @@
#include "src/__support/FPUtil/sqrt.h"
#endif // LIBC_TARGET_CPU_HAS_FPU_FLOAT || !LIBC_MATH_HAS_SMALL_TABLES
-#ifndef LIBC_MATH_HAS_SMALL_TABLES
+#if !defined(LIBC_MATH_HAS_SMALL_TABLES) && \
+ !defined(LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT)
#include "src/__support/math/exp10f.h"
#include "src/__support/math/exp2f.h"
-#endif // !LIBC_MATH_HAS_SMALL_TABLES
+#endif // !LIBC_MATH_HAS_SMALL_TABLES &&
+ // !LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT
namespace LIBC_NAMESPACE_DECL {
namespace math {
@@ -209,34 +211,35 @@ LIBC_INLINE_VAR constexpr TripleDouble LOG2_R_TD[128] = {
{0.0, 0.0, 1.0},
};
-// Lookup table for 5-bit (32 entries) range reduction for float_eval and
-// double_eval fast path:
+// Lookup table for 4-bit (16 entries) range reduction for float_eval:
// r(0) = 1.0f
-// r(i) = 2^-6 * ceil(2^6 * (1 - 2^-6) / (1 + i * 2^-5)), i = 1..31.
+// r(i) = 2^-5 * ceil(2^5 * (1 - 2^-5) / (1 + i * 2^-4)), i = 1..15.
// The constants are chosen so that dx = r * m_x - 1 is exact in
// single precision (with hardware FMA or 14-bit Sterbenz split), and
-// -2^-6 <= dx < 2^-5.
+// -2^-5 <= dx < 2^-4.
//
// Generated by Sollya with:
// > display = hexadecimal;
-// > for i from 0 to 31 do {
+// > for i from 0 to 15 do {
// if (i == 0) then {
// print("0x1.0p+0f,");
+// } else if (i == 15) then {
+// print("0x1.0p-1f,");
// } else {
-// r = 2^(-6) * ceil( 2^6 * (1 - 2^(-6)) / (1 + i*2^(-5)) );
+// r = 2^(-5) * ceil( 2^5 * (1 - 2^(-5)) / (1 + i*2^(-4)) );
// print(r @ "f,");
// };
// };
-LIBC_INLINE_VAR constexpr float R_32[32] = {
- 0x1.0p+0f, 0x1.fp-1f, 0x1.ep-1f, 0x1.dp-1f, 0x1.cp-1f, 0x1.b8p-1f,
- 0x1.bp-1f, 0x1.ap-1f, 0x1.98p-1f, 0x1.9p-1f, 0x1.8p-1f, 0x1.78p-1f,
- 0x1.7p-1f, 0x1.68p-1f, 0x1.6p-1f, 0x1.58p-1f, 0x1.5p-1f, 0x1.5p-1f,
- 0x1.48p-1f, 0x1.4p-1f, 0x1.38p-1f, 0x1.38p-1f, 0x1.3p-1f, 0x1.28p-1f,
- 0x1.2p-1f, 0x1.2p-1f, 0x1.18p-1f, 0x1.18p-1f, 0x1.1p-1f, 0x1.1p-1f,
- 0x1.08p-1f, 0x1.0p-1f,
+LIBC_INLINE_VAR constexpr float R_16[16] = {
+ 0x1.0p+0f, 0x1.ep-1f, 0x1.cp-1f, 0x1.bp-1f, 0x1.9p-1f, 0x1.8p-1f,
+ 0x1.7p-1f, 0x1.6p-1f, 0x1.5p-1f, 0x1.4p-1f, 0x1.4p-1f, 0x1.3p-1f,
+ 0x1.2p-1f, 0x1.2p-1f, 0x1.1p-1f, 0x1.0p-1f,
};
-// R_32 represented in double precision for double_eval.
+// Lookup table for 5-bit (32 entries) range reduction for double_eval fast
+// path:
+// r(0) = 1.0
+// r(i) = 2^-6 * ceil(2^6 * (1 - 2^-6) / (1 + i * 2^-5)), i = 1..31.
LIBC_INLINE_VAR constexpr double R_32_D[32] = {
0x1.0p+0, 0x1.fp-1, 0x1.ep-1, 0x1.dp-1, 0x1.cp-1, 0x1.b8p-1, 0x1.bp-1,
0x1.ap-1, 0x1.98p-1, 0x1.9p-1, 0x1.8p-1, 0x1.78p-1, 0x1.7p-1, 0x1.68p-1,
@@ -245,7 +248,7 @@ LIBC_INLINE_VAR constexpr double R_32_D[32] = {
0x1.1p-1, 0x1.1p-1, 0x1.08p-1, 0x1.0p-1,
};
-// Table of -log2(R_32[i]) represented in double precision for double_eval.
+// Table of -log2(R_32_D[i]) represented in double precision for double_eval.
//
// Generated by Sollya with:
// > display = hexadecimal;
@@ -274,46 +277,38 @@ LIBC_INLINE_VAR constexpr double LOG2_R_32[32] = {
0x1.e9452c8a71028p-1, 0x1.0000000000000p+0,
};
-// Table of -log2(R_32[i]) represented in FloatFloat precision {lo, hi} for
+// Table of -log2(R_16[i]) represented in FloatFloat precision {lo, hi} for
// float_eval.
// We choose the precision of the high part to be 24 - 8 = 16 bits, so that
-// e_xf + LOG2_R_FF_32[i].hi
+// e_xf + LOG2_R_FF_16[i].hi
// is exact in float for |e_x| <= 150.
-// For i = 31, R_32[31] = 0.5, so -log2(R_32[31]) = 1.0f is exact.
+// For i = 15, R_16[15] = 0.5, so -log2(R_16[15]) = 1.0f is exact.
//
// Generated by Sollya with:
// > display = hexadecimal;
// > prec = 200;
-// > for i from 0 to 31 do {
+// > for i from 0 to 15 do {
// if (i == 0) then {
// print("{0x0.000000p+0f, 0x0.000000p+0f},");
-// } else if (i == 31) then {
+// } else if (i == 15) then {
// print("{0x0.000000p+0f, 0x1.000000p+0f},");
// } else {
-// r = 2^(-6) * ceil( 2^6 * (1 - 2^(-6)) / (1 + i*2^(-5)) );
+// r = 2^(-5) * ceil( 2^5 * (1 - 2^(-5)) / (1 + i*2^(-4)) );
// l = -log2(r);
// h = round(1 + l, 17, RN) - 1;
// lo = round(l - h, SG, RN);
// print("{" @ lo @ "f, " @ h @ "f},");
// };
// };
-LIBC_INLINE_VAR constexpr FloatFloat LOG2_R_FF_32[32] = {
- {0x0.000000p+0f, 0x0.000000p+0f}, {-0x1.acd89ap-19f, 0x1.774p-5f},
- {0x1.25b3eep-22f, 0x1.7d6p-4f}, {0x1.6e155ap-18f, 0x1.22d8p-3f},
- {0x1.80abfcp-19f, 0x1.8a88p-3f}, {-0x1.858p-19f, 0x1.bfc8p-3f},
- {-0x1.3ab7cep-18f, 0x1.f6p-3f}, {-0x1.1c8f12p-22f, 0x1.32cp-2f},
- {-0x1.134c4ep-20f, 0x1.4f7p-2f}, {0x1.ed0cbap-19f, 0x1.6cbp-2f},
- {-0x1.a39fbep-20f, 0x1.a9p-2f}, {0x1.dc2d46p-18f, 0x1.c818p-2f},
- {-0x1.40358ep-19f, 0x1.e7ep-2f}, {-0x1.5a32c2p-20f, 0x1.042cp-1f},
- {-0x1.3e032ep-18f, 0x1.14c6p-1f}, {0x1.408c78p-18f, 0x1.25c0p-1f},
- {0x1.5db83ap-20f, 0x1.3722p-1f}, {0x1.5db83ap-20f, 0x1.3722p-1f},
- {0x1.e3263p-18f, 0x1.48eep-1f}, {0x1.ed0cbap-20f, 0x1.5b2cp-1f},
- {-0x1.eac382p-20f, 0x1.6de0p-1f}, {-0x1.eac382p-20f, 0x1.6de0p-1f},
- {-0x1.6b9026p-19f, 0x1.811p-1f}, {0x1.0e9258p-18f, 0x1.94c2p-1f},
+LIBC_INLINE_VAR constexpr FloatFloat LOG2_R_FF_16[16] = {
+ {0x0.000000p+0f, 0x0.000000p+0f}, {0x1.25b3eep-22f, 0x1.7d6p-4f},
+ {0x1.80abfcp-19f, 0x1.8a88p-3f}, {-0x1.3ab7cep-18f, 0x1.f6p-3f},
+ {0x1.ed0cbap-19f, 0x1.6cbp-2f}, {-0x1.a39fbep-20f, 0x1.a9p-2f},
+ {-0x1.40358ep-19f, 0x1.e7ep-2f}, {-0x1.3e032ep-18f, 0x1.14c6p-1f},
+ {0x1.5db83ap-20f, 0x1.3722p-1f}, {0x1.ed0cbap-20f, 0x1.5b2cp-1f},
+ {0x1.ed0cbap-20f, 0x1.5b2cp-1f}, {-0x1.6b9026p-19f, 0x1.811p-1f},
{-0x1.a39fbep-19f, 0x1.a9p-1f}, {-0x1.a39fbep-19f, 0x1.a9p-1f},
- {0x1.3b992cp-18f, 0x1.bdcep-1f}, {0x1.3b992cp-18f, 0x1.bdcep-1f},
- {0x1.20a6dep-21f, 0x1.d338p-1f}, {0x1.20a6dep-21f, 0x1.d338p-1f},
- {-0x1.a6eb1ep-18f, 0x1.e946p-1f}, {0x0.000000p+0f, 0x1.000000p+0f},
+ {0x1.20a6dep-21f, 0x1.d338p-1f}, {0x0.000000p+0f, 0x1.000000p+0f},
};
// Upper bound for y = 150 / |log2(1 - 2^-24)|, generated by Sollya:
@@ -378,7 +373,8 @@ LIBC_INLINE float set_overflow(Sign sign = Sign::POS) {
return sign.is_neg() ? FPBits::max_normal(Sign::NEG).get_val()
: FPBits::inf(Sign::POS).get_val();
#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- return FPBits::inf(sign).get_val();
+ return sign.is_neg() ? FPBits::inf(Sign::NEG).get_val()
+ : FPBits::inf(Sign::POS).get_val();
}
LIBC_INLINE float set_underflow(Sign sign = Sign::POS) {
@@ -392,7 +388,7 @@ LIBC_INLINE float set_underflow(Sign sign = Sign::POS) {
if (rounding == FE_DOWNWARD && sign.is_neg())
return FPBits::min_subnormal(Sign::NEG).get_val();
#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- return FPBits::zero(sign).get_val();
+ return sign.is_neg() ? -0.0f : 0.0f;
}
// Fast checks for special inputs:
@@ -427,7 +423,8 @@ LIBC_ALWAYS_INLINE cpp::optional<float> check_special_inputs(float x, float y) {
// pow(0, negative number) = Inf
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_DIVBYZERO);
- return FPBits::inf(out_is_neg ? Sign::NEG : Sign::POS).get_val();
+ return out_is_neg ? FPBits::inf(Sign::NEG).get_val()
+ : FPBits::inf(Sign::POS).get_val();
}
// pow(0, positive number) = 0
return out_is_neg ? -0.0f : 0.0f;
@@ -470,7 +467,8 @@ LIBC_ALWAYS_INLINE cpp::optional<float> check_special_inputs(float x, float y) {
}
}
-#ifndef LIBC_MATH_HAS_SMALL_TABLES
+#if !defined(LIBC_MATH_HAS_SMALL_TABLES) && \
+ !defined(LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT)
if (x == 2.0f)
return math::exp2f(y);
@@ -493,7 +491,8 @@ LIBC_ALWAYS_INLINE cpp::optional<float> check_special_inputs(float x, float y) {
return math::exp2f(hi);
}
}
-#endif // !LIBC_MATH_HAS_SMALL_TABLES
+#endif // !LIBC_MATH_HAS_SMALL_TABLES &&
+ // !LIBC_MATH_HAS_INTERMEDIATE_COMP_IN_FLOAT
return cpp::nullopt;
}
@@ -522,6 +521,9 @@ check_exceptional_cases(float &x, float &y, int &ex, SignType &sign) {
return FloatBits::quiet_nan().get_val();
}
+ if (x == 1.0f || y == 0.0f)
+ return 1.0f;
+
// Extreme |y| >= Y_UPPER_BOUND: |y * log2(x)| >= 150 for all x != 1.
if (y_a >= Y_UPPER_BOUND) {
if (x_abs.is_nan())
@@ -540,59 +542,64 @@ check_exceptional_cases(float &x, float &y, int &ex, SignType &sign) {
}
// y is finite and non-zero.
-
- if (x_a == FloatBits::inf().uintval()) {
- // pow(+-Inf, y):
- // - y < 0: returns +-0.0 (negative if x = -Inf and y is odd integer).
- // - y > 0: returns +-Inf (negative if x = -Inf and y is odd integer).
- bool out_is_neg = x_sign && is_odd_integer(y);
- Sign out_sign = out_is_neg ? Sign::NEG : Sign::POS;
- return y_sign ? FloatBits::zero(out_sign).get_val()
- : FloatBits::inf(out_sign).get_val();
- }
-
if (x_a > FloatBits::inf().uintval()) {
// x is NaN.
return x;
}
- if (x_a == 0) {
- bool out_is_neg = x_sign && is_odd_integer(y);
- if (y_sign) {
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_DIVBYZERO);
- return FloatBits::inf(out_is_neg ? Sign::NEG : Sign::POS).get_val();
- }
- return out_is_neg ? -0.0f : 0.0f;
- }
-
- // Handle negative base x < 0:
+ // Handle negative base x < 0 (including -0.0f and -Inf):
// - If y is an integer: (-|x|)^y = (-1)^y * |x|^y.
- // - If y is not an integer: (-|x|)^y is undefined in real numbers;
- // raises FE_INVALID, sets errno to EDOM, and returns quiet NaN.
+ // - If y is not an integer and x is finite non-zero: (-|x|)^y is undefined
+ // in real numbers; raises FE_INVALID, sets errno to EDOM, and returns
+ // quiet NaN.
if (x_sign) {
if (is_integer(y)) {
- x = -x;
if (is_odd_integer(y)) {
if constexpr (sizeof(SignType) == 8)
sign = 0x8000'0000'0000'0000ULL;
else
sign = 0x8000'0000U;
}
- } else {
+ } else if (x != 0.0f && x_a != FloatBits::inf().uintval()) {
// pow( negative, non-integer ) = NaN
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
return FloatBits::quiet_nan().get_val();
}
+ x = -x;
+ }
+
+ if (x_a == FloatBits::inf().uintval()) {
+ // pow(+-Inf, y):
+ // - y < 0: returns +-0.0 (negative if x = -Inf and y is odd integer).
+ // - y > 0: returns +-Inf (negative if x = -Inf and y is odd integer).
+ if (y_sign)
+ return (sign != 0) ? -0.0f : 0.0f;
+ return (sign != 0) ? FloatBits::inf(Sign::NEG).get_val()
+ : FloatBits::inf(Sign::POS).get_val();
+ }
+
+ if (x == 0.0f) {
+ if (y_sign) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ return (sign != 0) ? FloatBits::inf(Sign::NEG).get_val()
+ : FloatBits::inf(Sign::POS).get_val();
+ }
+ return (sign != 0) ? -0.0f : 0.0f;
}
if (y_a <= Y_LOWER_BOUND) {
+#if defined(LIBC_MATH_HAS_NO_EXCEPT) && \
+ defined(LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY)
+ return 1.0f;
+#else
volatile float one = 1.0f;
volatile float eps = ((x_a < FloatBits::one().uintval()) == !y_sign)
? -0x1.0p-50f
: 0x1.0p-50f;
return one + eps;
+#endif
}
// Normalize denormal inputs.
diff --git a/libc/test/src/math/smoke/powf_test.cpp b/libc/test/src/math/smoke/powf_test.cpp
index 46d4f1b7120944..d7ce7e335bbc19 100644
--- a/libc/test/src/math/smoke/powf_test.cpp
+++ b/libc/test/src/math/smoke/powf_test.cpp
@@ -186,9 +186,9 @@ class PowfTest : public LIBC_NAMESPACE::testing::FPTest<float> {
// Powers of 10:
EXPECT_FP_EQ(1.0f, func(10.0f, 0.0f));
- EXPECT_FP_EQ(10.0f, func(10.0f, 1.0f));
- EXPECT_FP_EQ(100.0f, func(10.0f, 2.0f));
if (tolerance == 0) {
+ EXPECT_FP_EQ(10.0f, func(10.0f, 1.0f));
+ EXPECT_FP_EQ(100.0f, func(10.0f, 2.0f));
EXPECT_FP_EQ(1000.0f, func(10.0f, 3.0f));
EXPECT_FP_EQ(10000.0f, func(10.0f, 4.0f));
EXPECT_FP_EQ(100000.0f, func(10.0f, 5.0f));
@@ -203,6 +203,8 @@ class PowfTest : public LIBC_NAMESPACE::testing::FPTest<float> {
uint32_t e = FPBits(exp).uintval();
EXPECT_LE(a >= e ? a - e : e - a, static_cast<uint32_t>(tolerance));
};
+ check_close(func(10.0f, 1.0f), 10.0f);
+ check_close(func(10.0f, 2.0f), 100.0f);
check_close(func(10.0f, 3.0f), 1000.0f);
check_close(func(10.0f, 4.0f), 10000.0f);
check_close(func(10.0f, 5.0f), 100000.0f);
>From 1db227ed62e2633731c55940e19b7f0559b7dcf4 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue at google.com>
Date: Fri, 16 Oct 2026 01:57:01 -0400
Subject: [PATCH 2/2] Fix unused using.
---
libc/src/__support/math/powf_utils.h | 4 +++-
1 file changed, 3 insertions(+), 1 deletion(-)
diff --git a/libc/src/__support/math/powf_utils.h b/libc/src/__support/math/powf_utils.h
index cb0efe6bcacd90..406ab3680f4c08 100644
--- a/libc/src/__support/math/powf_utils.h
+++ b/libc/src/__support/math/powf_utils.h
@@ -380,14 +380,16 @@ LIBC_INLINE float set_overflow(Sign sign = Sign::POS) {
LIBC_INLINE float set_underflow(Sign sign = Sign::POS) {
fputil::set_errno_if_required(ERANGE);
fputil::raise_underflow_except_if_required<float>();
- using FPBits = fputil::FPBits<float>;
+
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ using FPBits = fputil::FPBits<float>;
int rounding = fputil::quick_get_round();
if (rounding == FE_UPWARD && sign.is_pos())
return FPBits::min_subnormal(Sign::POS).get_val();
if (rounding == FE_DOWNWARD && sign.is_neg())
return FPBits::min_subnormal(Sign::NEG).get_val();
#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+
return sign.is_neg() ? -0.0f : 0.0f;
}
More information about the libc-commits
mailing list