[libc-commits] [libc] [libc][math] Optimize exp2f and exp10f hot paths (PR #224178)
via libc-commits
libc-commits at lists.llvm.org
Fri Sep 25 08:57:13 PDT 2026
https://github.com/sriramshastry updated https://github.com/llvm/llvm-project/pull/224178
>From f00e87f44531456e110f547173e2dff8203a3b9e Mon Sep 17 00:00:00 2001
From: Sriram Shastry <sriramshastry at gmail.com>
Date: Wed, 16 Sep 2026 23:15:42 +0530
Subject: [PATCH] [libc][math] Optimize exp2f and exp10f hot paths
Route common finite inputs directly to the existing degree-5 kernels so
exp2f avoids its exceptional-path frame and exp10f keeps cold handling out
of line. Preserve the proven polynomial arithmetic and special cases.
Add official range-based performance coverage for both functions.
PerfTest medians from alternating baseline/patch runs were:
baseline patch speedup
exp2f hot normal 2.33805 1.92944 17.02%
exp2f close to one 2.28649 1.76252 22.78%
exp10f hot normal 2.83441 2.36588 15.80%
exp10f close to one 2.89851 2.25680 20.67%
Times are ns/op. Each row pools positive and negative ranges, with six
alternating pairs per range. Official exhaustive MPFR tests pass both sign
ranges in all four rounding modes for each function.
Signed-off-by: Sriram Shastry <sriramshastry at gmail.com>
---
libc/src/__support/math/exp10f_double_eval.h | 147 ++++++++++--------
libc/src/__support/math/exp2f_double_eval.h | 96 ++++++------
.../math/performance_testing/CMakeLists.txt | 11 ++
.../math/performance_testing/exp10f_perf.cpp | 52 +++++++
.../math/performance_testing/exp2f_perf.cpp | 40 ++++-
5 files changed, 229 insertions(+), 117 deletions(-)
create mode 100644 libc/test/src/math/performance_testing/exp10f_perf.cpp
diff --git a/libc/src/__support/math/exp10f_double_eval.h b/libc/src/__support/math/exp10f_double_eval.h
index e1a7cc9d68e98..eec615ef40429 100644
--- a/libc/src/__support/math/exp10f_double_eval.h
+++ b/libc/src/__support/math/exp10f_double_eval.h
@@ -27,67 +27,49 @@ namespace LIBC_NAMESPACE_DECL {
namespace math {
namespace double_eval {
-LIBC_INLINE float exp10f(float x) {
+LIBC_INLINE float exp10f_mid(float x) {
+ // Range reduction: 10^x = 2^(mid + hi) * 10^lo
+ auto rr = exp_b_range_reduc<Exp10Base>(x);
+
+ // The low part is approximated by a degree-5 minimax polynomial.
+ using fputil::multiply_add;
+ double lo2 = rr.lo * rr.lo;
+ double c0 = multiply_add(rr.lo, Exp10Base::COEFFS[0], 1.0);
+ double c1 = multiply_add(rr.lo, Exp10Base::COEFFS[2], Exp10Base::COEFFS[1]);
+ double c2 = multiply_add(rr.lo, Exp10Base::COEFFS[4], Exp10Base::COEFFS[3]);
+ double p = multiply_add(lo2, c2, c1);
+ // 10^x = 2^(mid + hi) * 10^lo
+ // ~ mh * (1 + COEFFS[0] * lo + ... + COEFFS[4] * lo^5)
+ return static_cast<float>(multiply_add(p, lo2 * rr.mh, c0 * rr.mh));
+}
+
+[[gnu::cold, gnu::noinline]] static float exp10f_slow(float x) {
using FPBits = fputil::FPBits<float>;
FPBits xbits(x);
uint32_t x_u = xbits.uintval();
uint32_t x_abs = x_u & 0x7fff'ffffU;
- // When |x| >= log10(2^128), or x is nan
- if (LIBC_UNLIKELY(x_abs >= 0x421a'209bU)) {
- // When x < log10(2^-150) or nan
- if (x_u > 0xc234'9e35U) {
- // exp(-Inf) = 0
- if (xbits.is_inf())
- return 0.0f;
- // exp(nan) = nan
- if (xbits.is_nan())
- return x;
+ // x >= log10(2^128), or positive NaN.
+ if (LIBC_UNLIKELY(x_u >= 0x421a'209bU && x_u < 0x8000'0000U)) {
+ // x is finite.
+ if (x_u < 0x7f80'0000U) {
#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- if (fputil::fenv_is_round_up())
- return FPBits::min_subnormal().get_val();
-#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_except_if_required(FE_UNDERFLOW);
- return 0.0f;
- }
- // x >= log10(2^128) or nan
- if (xbits.is_pos() && (x_u >= 0x421a'209bU)) {
- // x is finite
- if (x_u < 0x7f80'0000U) {
-#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- int rounding = fputil::quick_get_round();
- if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
- return FPBits::max_normal().get_val();
+ int rounding = fputil::quick_get_round();
+ if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
+ return FPBits::max_normal().get_val();
#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_except_if_required(FE_OVERFLOW);
- }
- // x is +inf or nan
- return x + FPBits::inf().get_val();
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_OVERFLOW);
}
+ // x is +inf or nan
+ return x + FPBits::inf().get_val();
}
// 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_ASSUME_ROUND_NEAREST_ONLY
- return 0x1.fffffep-1f;
-#else // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- if (fputil::fenv_is_round_to_nearest())
- return 0x1.fffffep-1f;
-#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- }
- // |x| < 2^-25
- // 10^x ~ 1 + log(10) * x
- if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
- return fputil::multiply_add(x, 0x1.26bb1cp+1f, 1.0f);
- }
-
+ if (LIBC_UNLIKELY(x_abs <= 0x3b9a'209bU))
return static_cast<float>(Exp10Base::powb_lo(x));
- }
// Exceptional value.
if (LIBC_UNLIKELY(x_u == 0x3d14'd956U)) { // x = 0x1.29b2acp-5f
@@ -124,28 +106,57 @@ LIBC_INLINE float exp10f(float x) {
}
}
- // Range reduction: 10^x = 2^(mid + hi) * 10^lo
- // rr = (2^(mid + hi), lo)
- auto rr = exp_b_range_reduc<Exp10Base>(x);
+ // When |x| >= log10(2^128), or x is nan
+ if (LIBC_UNLIKELY(x_abs >= 0x421a'209bU)) {
+ // When x < log10(2^-150) or nan
+ if (x_u > 0xc234'9e35U) {
+ // exp(-Inf) = 0
+ if (xbits.is_inf())
+ return 0.0f;
+ // exp(nan) = nan
+ if (xbits.is_nan())
+ return x;
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (fputil::fenv_is_round_up())
+ return FPBits::min_subnormal().get_val();
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_UNDERFLOW);
+ return 0.0f;
+ }
+ }
- // The low part is approximated by a degree-5 minimax polynomial.
- // 10^lo ~ 1 + COEFFS[0] * lo + ... + COEFFS[4] * lo^5
- using fputil::multiply_add;
- double lo2 = rr.lo * rr.lo;
- // c0 = 1 + COEFFS[0] * lo
- double c0 = multiply_add(rr.lo, Exp10Base::COEFFS[0], 1.0);
- // c1 = COEFFS[1] + COEFFS[2] * lo
- double c1 = multiply_add(rr.lo, Exp10Base::COEFFS[2], Exp10Base::COEFFS[1]);
- // c2 = COEFFS[3] + COEFFS[4] * lo
- double c2 = multiply_add(rr.lo, Exp10Base::COEFFS[4], Exp10Base::COEFFS[3]);
- // p = c1 + c2 * lo^2
- // = COEFFS[1] + COEFFS[2] * lo + COEFFS[3] * lo^2 + COEFFS[4] * lo^3
- double p = multiply_add(lo2, c2, c1);
- // 10^lo ~ c0 + p * lo^2
- // 10^x = 2^(mid + hi) * 10^lo
- // ~ mh * (c0 + p * lo^2)
- // = (mh * c0) + p * (mh * lo^2)
- return static_cast<float>(multiply_add(p, lo2 * rr.mh, c0 * rr.mh));
+ return exp10f_mid(x);
+}
+
+LIBC_INLINE float exp10f(float x) {
+ using FPBits = fputil::FPBits<float>;
+ FPBits xbits(x);
+
+ uint32_t x_u = xbits.uintval();
+ uint32_t x_abs = x_u & 0x7fff'ffffU;
+
+ if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
+ if (LIBC_UNLIKELY(x_u == 0xb25e'5bd9U)) { // x = -0x1.bcb7b2p-27f
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ return 0x1.fffffep-1f;
+#else // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (fputil::fenv_is_round_to_nearest())
+ return 0x1.fffffep-1f;
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ }
+ // |x| < 2^-25
+ // 10^x ~ 1 + log(10) * x
+ return fputil::multiply_add(x, 0x1.26bb1cp+1f, 1.0f);
+ }
+
+ if (LIBC_LIKELY(
+ x_abs > 0x3b9a'209bU &&
+ (x_u < 0x421a'209bU || (xbits.is_neg() && x_u <= 0xc234'9e35U)) &&
+ x_u != 0x3d14'd956U && (x_u & 0x800f'ffffU) != 0))
+ return exp10f_mid(x);
+
+ return exp10f_slow(x);
}
} // namespace double_eval
diff --git a/libc/src/__support/math/exp2f_double_eval.h b/libc/src/__support/math/exp2f_double_eval.h
index e0387f32e193f..c9548712ab379 100644
--- a/libc/src/__support/math/exp2f_double_eval.h
+++ b/libc/src/__support/math/exp2f_double_eval.h
@@ -38,15 +38,59 @@ LIBC_INLINE float exp2f(float x) {
uint32_t x_u = xbits.uintval();
uint32_t x_abs = x_u & 0x7fff'ffffU;
+ // |x| < 2^-25
+ if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U))
+ return 1.0f + x;
+
+ auto exp2f_mid = [&]() -> float {
+ // kf = (hi + mid) * 2^5 = round(x * 2^5)
+ float kf = 0;
+ int k = 0;
+#ifdef LIBC_TARGET_CPU_HAS_NEAREST_INT
+ kf = fputil::nearest_integer(x * 32.0f);
+ k = static_cast<int>(kf);
+#else // !LIBC_TARGET_CPU_HAS_NEAREST_INT
+ constexpr float HALF[2] = {0.5f, -0.5f};
+ k = static_cast<int>(fputil::multiply_add(x, 32.0f, HALF[x < 0.0f]));
+ kf = static_cast<float>(k);
+#endif // LIBC_TARGET_CPU_HAS_NEAREST_INT
+
+ // dx = lo = x - (hi + mid) = x - kf * 2^(-5)
+ double dx = fputil::multiply_add(-0x1.0p-5f, kf, x);
+
+ // hi = floor(kf * 2^(-4))
+ // exp_hi = shift hi to the exponent field of double precision.
+ int64_t exp_hi =
+ static_cast<int64_t>(static_cast<uint64_t>(k >> ExpBase::MID_BITS)
+ << fputil::FPBits<double>::FRACTION_LEN);
+ // mh = 2^hi * 2^mid
+ // mh_bits = bit field of mh
+ int64_t mh_bits = ExpBase::EXP_2_MID[k & ExpBase::MID_MASK] + exp_hi;
+ double mh = fputil::FPBits<double>(uint64_t(mh_bits)).get_val();
+
+ // Degree-5 polynomial approximating (2^x - 1)/x generated by Sollya.
+ constexpr double COEFFS[5] = {0x1.62e42fefa39efp-1, 0x1.ebfbdff8131c4p-3,
+ 0x1.c6b08d7061695p-5, 0x1.3b2b1bee74b2ap-7,
+ 0x1.5d88091198529p-10};
+ double dx_sq = dx * dx;
+ double c1 = fputil::multiply_add(dx, COEFFS[0], 1.0);
+ double c2 = fputil::multiply_add(dx, COEFFS[2], COEFFS[1]);
+ double c3 = fputil::multiply_add(dx, COEFFS[4], COEFFS[3]);
+ double p = fputil::multiply_add(dx_sq, c3, c2);
+ // 2^x = 2^(hi + mid + lo)
+ // = 2^(hi + mid) * 2^lo
+ // ~ mh * (1 + lo * P(lo))
+ // = mh + (mh*lo) * P(lo)
+ return static_cast<float>(fputil::multiply_add(p, dx_sq * mh, c1 * mh));
+ };
+
+ if (LIBC_LIKELY(x_abs > 0x3d00'0000U && x_abs < 0x4300'0000U))
+ return exp2f_mid();
+
// When |x| >= 128, or x is nan, or |x| <= 2^-5
if (LIBC_UNLIKELY(x_abs >= 0x4300'0000U || x_abs <= 0x3d00'0000U)) {
// |x| <= 2^-5
if (x_abs <= 0x3d00'0000) {
- // |x| < 2^-25
- if (LIBC_UNLIKELY(x_abs <= 0x3280'0000U)) {
- return 1.0f + x;
- }
-
#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
constexpr uint32_t EXVAL1 = 0x3b42'9d37U;
constexpr uint32_t EXVAL2 = 0xbcf3'a937U;
@@ -128,47 +172,7 @@ LIBC_INLINE float exp2f(float x) {
// generated by Sollya.
// We perform 2^hi * 2^mid by simply add hi to the exponent field
// of 2^mid.
-
- // kf = (hi + mid) * 2^5 = round(x * 2^5)
- float kf = 0;
- int k = 0;
-#ifdef LIBC_TARGET_CPU_HAS_NEAREST_INT
- kf = fputil::nearest_integer(x * 32.0f);
- k = static_cast<int>(kf);
-#else // !LIBC_TARGET_CPU_HAS_NEAREST_INT
- constexpr float HALF[2] = {0.5f, -0.5f};
- k = static_cast<int>(fputil::multiply_add(x, 32.0f, HALF[x < 0.0f]));
- kf = static_cast<float>(k);
-#endif // LIBC_TARGET_CPU_HAS_NEAREST_INT
-
- // dx = lo = x - (hi + mid) = x - kf * 2^(-5)
- double dx = fputil::multiply_add(-0x1.0p-5f, kf, x);
-
- // hi = floor(kf * 2^(-4))
- // exp_hi = shift hi to the exponent field of double precision.
- int64_t exp_hi =
- static_cast<int64_t>(static_cast<uint64_t>(k >> ExpBase::MID_BITS)
- << fputil::FPBits<double>::FRACTION_LEN);
- // mh = 2^hi * 2^mid
- // mh_bits = bit field of mh
- int64_t mh_bits = ExpBase::EXP_2_MID[k & ExpBase::MID_MASK] + exp_hi;
- double mh = fputil::FPBits<double>(uint64_t(mh_bits)).get_val();
-
- // Degree-5 polynomial approximating (2^x - 1)/x generating by Sollya with:
- // > P = fpminimax((2^x - 1)/x, 5, [|D...|], [-1/32. 1/32]);
- constexpr double COEFFS[5] = {0x1.62e42fefa39efp-1, 0x1.ebfbdff8131c4p-3,
- 0x1.c6b08d7061695p-5, 0x1.3b2b1bee74b2ap-7,
- 0x1.5d88091198529p-10};
- double dx_sq = dx * dx;
- double c1 = fputil::multiply_add(dx, COEFFS[0], 1.0);
- double c2 = fputil::multiply_add(dx, COEFFS[2], COEFFS[1]);
- double c3 = fputil::multiply_add(dx, COEFFS[4], COEFFS[3]);
- double p = fputil::multiply_add(dx_sq, c3, c2);
- // 2^x = 2^(hi + mid + lo)
- // = 2^(hi + mid) * 2^lo
- // ~ mh * (1 + lo * P(lo))
- // = mh + (mh*lo) * P(lo)
- return static_cast<float>(fputil::multiply_add(p, dx_sq * mh, c1 * mh));
+ return exp2f_mid();
}
} // namespace double_eval
diff --git a/libc/test/src/math/performance_testing/CMakeLists.txt b/libc/test/src/math/performance_testing/CMakeLists.txt
index 534bd4384e39a..2053c68d5a4a4 100644
--- a/libc/test/src/math/performance_testing/CMakeLists.txt
+++ b/libc/test/src/math/performance_testing/CMakeLists.txt
@@ -155,6 +155,17 @@ add_perf_binary(
-fno-builtin
)
+add_perf_binary(
+ exp10f_perf
+ SRCS
+ exp10f_perf.cpp
+ DEPENDS
+ .perf_test
+ libc.src.math.exp10f
+ COMPILE_OPTIONS
+ -fno-builtin
+)
+
add_perf_binary(
exp2f_perf
SRCS
diff --git a/libc/test/src/math/performance_testing/exp10f_perf.cpp b/libc/test/src/math/performance_testing/exp10f_perf.cpp
new file mode 100644
index 0000000000000..76dab21781111
--- /dev/null
+++ b/libc/test/src/math/performance_testing/exp10f_perf.cpp
@@ -0,0 +1,52 @@
+//===-- Differential test for exp10f --------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+
+#include "PerfTest.h"
+#include "src/math/exp10f.h"
+
+#include <fstream>
+#include <math.h>
+
+int main(int argc, char **argv) {
+ using Perf = LIBC_NAMESPACE::testing::PerfTest<float, float>;
+ using FuncPtr = Perf::UnaryFuncPtr;
+ auto libc_func = static_cast<FuncPtr>(&LIBC_NAMESPACE::exp10f);
+ auto system_func = static_cast<FuncPtr>(&::exp10f);
+
+ constexpr size_t SAMPLE_COUNT = 1'000'001;
+ constexpr size_t ROUND_COUNT = 20;
+ constexpr char LOG_FILE[] = "exp10f_perf.log";
+ char selected_range = argc > 1 ? argv[1][0] : '\0';
+
+ if (selected_range == '\0')
+ Perf::run_perf<false>(libc_func, system_func, 10, "LLVM libc exp10f",
+ "system libc exp10f", LOG_FILE);
+
+ std::ofstream log(LOG_FILE,
+ selected_range == '\0' ? std::ios::app : std::ios::trunc);
+ auto run_range = [&](char id, const char *name, uint32_t start,
+ uint32_t stop) {
+ if (selected_range != '\0' && selected_range != id)
+ return;
+ log << "\n " << name << ":\n";
+ Perf::run_perf_in_range<false>(
+ libc_func, system_func, start, stop, SAMPLE_COUNT, ROUND_COUNT,
+ "LLVM libc exp10f", "system libc exp10f", log);
+ };
+
+ run_range('1', "Normal outputs, nonnegative inputs", 0x0000'0000U,
+ 0x421a'209aU);
+ run_range('2', "Normal outputs, negative inputs", 0x8000'0000U, 0xc217'b818U);
+ run_range('3', "Subnormal outputs", 0xc217'b819U, 0xc234'9e35U);
+ run_range('4', "Close to one, nonnegative inputs", 0x0000'0000U,
+ 0x3b9a'209bU);
+ run_range('5', "Close to one, negative inputs", 0x8000'0000U, 0xbb9a'209bU);
+ run_range('6', "Overflow inputs", 0x421a'209bU, 0x7f7f'ffffU);
+ run_range('7', "Underflow inputs", 0xc234'9e36U, 0xff7f'ffffU);
+ return 0;
+}
diff --git a/libc/test/src/math/performance_testing/exp2f_perf.cpp b/libc/test/src/math/performance_testing/exp2f_perf.cpp
index fa5a6ad452e94..571db1490baa6 100644
--- a/libc/test/src/math/performance_testing/exp2f_perf.cpp
+++ b/libc/test/src/math/performance_testing/exp2f_perf.cpp
@@ -9,10 +9,44 @@
#include "PerfTest.h"
#include "src/math/exp2f.h"
+#include <fstream>
#include <math.h>
-int main() {
- SINGLE_INPUT_SINGLE_OUTPUT_PERF(float, LIBC_NAMESPACE::exp2f, ::exp2f,
- "exp2f_perf.log")
+int main(int argc, char **argv) {
+ using Perf = LIBC_NAMESPACE::testing::PerfTest<float, float>;
+ using FuncPtr = Perf::UnaryFuncPtr;
+ auto libc_func = static_cast<FuncPtr>(&LIBC_NAMESPACE::exp2f);
+ auto system_func = static_cast<FuncPtr>(&::exp2f);
+
+ constexpr size_t SAMPLE_COUNT = 1'000'001;
+ constexpr size_t ROUND_COUNT = 20;
+ constexpr char LOG_FILE[] = "exp2f_perf.log";
+ char selected_range = argc > 1 ? argv[1][0] : '\0';
+
+ if (selected_range == '\0')
+ Perf::run_perf<false>(libc_func, system_func, 10, "LLVM libc exp2f",
+ "system libc exp2f", LOG_FILE);
+
+ std::ofstream log(LOG_FILE,
+ selected_range == '\0' ? std::ios::app : std::ios::trunc);
+ auto run_range = [&](char id, const char *name, uint32_t start,
+ uint32_t stop) {
+ if (selected_range != '\0' && selected_range != id)
+ return;
+ log << "\n " << name << ":\n";
+ Perf::run_perf_in_range<false>(libc_func, system_func, start, stop,
+ SAMPLE_COUNT, ROUND_COUNT, "LLVM libc exp2f",
+ "system libc exp2f", log);
+ };
+
+ run_range('1', "Normal outputs, nonnegative inputs", 0x0000'0000U,
+ 0x42ff'ffffU);
+ run_range('2', "Normal outputs, negative inputs", 0x8000'0000U, 0xc2fc'0000U);
+ run_range('3', "Subnormal outputs", 0xc2fc'0001U, 0xc315'ffffU);
+ run_range('4', "Close to one, nonnegative inputs", 0x0000'0000U,
+ 0x3c80'0000U);
+ run_range('5', "Close to one, negative inputs", 0x8000'0000U, 0xbc80'0000U);
+ run_range('6', "Overflow inputs", 0x4300'0000U, 0x7f7f'ffffU);
+ run_range('7', "Underflow inputs", 0xc316'0000U, 0xff7f'ffffU);
return 0;
}
More information about the libc-commits
mailing list