[libc-commits] [libc] [libc][math] Integer-only implementation of expf with 1-ULP errors (PR #209406)
Hoàng Minh Thiên via libc-commits
libc-commits at lists.llvm.org
Fri Jul 24 02:08:04 PDT 2026
https://github.com/hmthien050209 updated https://github.com/llvm/llvm-project/pull/209406
>From a84bf583dcdfd32ff76cddd703f7f5b80e66936d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Tue, 14 Jul 2026 21:35:37 +0700
Subject: [PATCH 1/5] feat: initialize compilable stubs for expf
---
libc/src/__support/math/CMakeLists.txt | 31 +++++++++
libc/src/__support/math/exp_integer_utils.h | 42 ++++++++++++
libc/src/__support/math/expf_integer_eval.h | 75 +++++++++++++++++++++
libc/src/math/generic/CMakeLists.txt | 1 +
libc/src/math/generic/expf.cpp | 11 ++-
5 files changed, 159 insertions(+), 1 deletion(-)
create mode 100644 libc/src/__support/math/exp_integer_utils.h
create mode 100644 libc/src/__support/math/expf_integer_eval.h
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 9f4624682eafe..b36cc1f537c69 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3060,6 +3060,22 @@ add_header_library(
libc.src.__support.macros.config
)
+add_header_library(
+ expf_integer_eval
+ HDRS
+ expf_integer_eval.h
+ DEPENDS
+ # TODO: strip the below deps down if able to
+ .exp_constants
+ .exp_integer_utils
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.frac128
+)
+
add_header_library(
expf
HDRS
@@ -4099,6 +4115,21 @@ add_header_library(
libc.src.__support.macros.optimization
)
+add_header_library(
+ exp_integer_utils
+ HDRS
+ exp_integer_utils.h
+ DEPENDS
+ # TODO: strip the below deps down if able to
+ .exp_constants
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.frac128
+)
+
add_header_library(
exp2
HDRS
diff --git a/libc/src/__support/math/exp_integer_utils.h b/libc/src/__support/math/exp_integer_utils.h
new file mode 100644
index 0000000000000..d6216b6ee183f
--- /dev/null
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -0,0 +1,42 @@
+//===-- e^x range reduction and evaluation using integer-only --*- C++ -*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+// TODO: update the description of this file above
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_UTILS_H
+
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/big_int.h"
+#include "src/__support/frac128.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math_extras.h"
+
+#undef LIBC_TARGET_IS_BIG_ENDIAN
+#if !defined(__BYTE_ORDER__) || !defined(__ORDER_LITTLE_ENDIAN__) || \
+ !defined(__ORDER_BIG_ENDIAN__)
+#define LIBC_TARGET_IS_BIG_ENDIAN 0
+#else
+#define LIBC_TARGET_IS_BIG_ENDIAN (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
+#endif // /LIBC_TARGET_IS_BIG_ENDIAN
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace integer_only {} // namespace integer_only
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_UTILS_H
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
new file mode 100644
index 0000000000000..1bab06719d0c8
--- /dev/null
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -0,0 +1,75 @@
+//===-- Implementation header for expf using integer-only --------*- C++-*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
+
+// TODO: clean up includes
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/frac128.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/exp_integer_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace integer_only {
+
+LIBC_INLINE float expf(float x) {
+ using FPBits = typename fputil::FPBits<float>;
+ FPBits xbits(x);
+
+ bool is_neg = xbits.is_neg();
+ uint16_t x_e = xbits.get_biased_exponent();
+ // uint64_t x_u = xbits.get_mantissa();
+
+ // Exceptional values
+ // TODO: optimize the branching order
+ if (xbits.is_zero()) {
+ // e^0 = 1 (exact), for both +/-0
+ return 1.0f;
+ } else {
+ // x is inf or NaN
+ if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
+ // e^NaN = NaN
+ if (xbits.is_signaling_nan()) {
+ // Per conversation with lntue, we don't need to raise exception here,
+ // as we're assuming no FPUs/fenv in this kind of environment
+
+ // silencing
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // e^-inf = 0
+ // e^+inf = +inf
+ if (xbits.is_inf()) {
+ return is_neg ? 0.0f : FPBits::inf().get_val();
+ }
+
+ // x is a quiet NaN
+ return x;
+ } else {
+ // TODO: out-of-range (overflow/underflow), exeecute range reduction
+ }
+ }
+
+ // TODO: execute the normal exp here
+ return 0.0f;
+}
+
+} // namespace integer_only
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index d5661f6f55e91..8fff101fa96a6 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1227,6 +1227,7 @@ add_entrypoint_object(
../expf.h
DEPENDS
libc.src.__support.math.expf
+ libc.src.__support.math.expf_integer_eval
libc.src.errno.errno
)
diff --git a/libc/src/math/generic/expf.cpp b/libc/src/math/generic/expf.cpp
index de11f51ac64a0..de8a97a46f1af 100644
--- a/libc/src/math/generic/expf.cpp
+++ b/libc/src/math/generic/expf.cpp
@@ -8,9 +8,18 @@
#include "src/math/expf.h"
#include "src/__support/math/expf.h"
+#include "src/__support/math/expf_integer_eval.h"
namespace LIBC_NAMESPACE_DECL {
-LLVM_LIBC_FUNCTION(float, expf, (float x)) { return math::expf(x); }
+LLVM_LIBC_FUNCTION(float, expf, (float x)) {
+#if defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) && \
+ defined(LIBC_MATH_SMALL_TABLES) && \
+ !defined(LIBC_TARGET_CPU_HAS_FPU_DOUBLE)
+ return math::integer_only::expf(x);
+#else
+ return math::expf(x);
+#endif
+}
} // namespace LIBC_NAMESPACE_DECL
>From 8b5848cf73f2a35162c14e91dca3b071b24693d9 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Wed, 15 Jul 2026 23:37:57 +0700
Subject: [PATCH 2/5] more stubs
---
libc/src/__support/frac64.h | 59 +++++++++++++++++++++
libc/src/__support/math/exp_integer_utils.h | 26 +++++++--
libc/src/__support/math/expf_integer_eval.h | 53 ++++++++++--------
3 files changed, 112 insertions(+), 26 deletions(-)
create mode 100644 libc/src/__support/frac64.h
diff --git a/libc/src/__support/frac64.h b/libc/src/__support/frac64.h
new file mode 100644
index 0000000000000..2172f4cbaeeb8
--- /dev/null
+++ b/libc/src/__support/frac64.h
@@ -0,0 +1,59 @@
+//===-- 64-bit unsigned fractional type -------------------------*- C++ -*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_FRAC64_H
+#define LLVM_LIBC_SRC___SUPPORT_FRAC64_H
+
+#include "hdr/stdint_proxy.h" // uint64_t
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_INLINE
+#include "src/__support/uint128.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+// Q0.64
+struct Frac64 {
+ uint64_t val;
+
+ LIBC_INLINE constexpr Frac64() : val(0) {}
+ LIBC_INLINE constexpr explicit Frac64(uint64_t value) : val(value) {}
+
+ LIBC_INLINE constexpr Frac64 operator~() const { return Frac64(~val); }
+
+ LIBC_INLINE constexpr Frac64 operator+(const Frac64 &other) const {
+ return Frac64(val + other.val);
+ }
+
+ LIBC_INLINE constexpr Frac64 operator-(const Frac64 &other) const {
+ return Frac64(val - other.val);
+ }
+
+ LIBC_INLINE constexpr Frac64 operator*(const Frac64 &other) const {
+ return Frac64(
+ static_cast<uint64_t>((static_cast<UInt128>(val) * other.val) >> 64));
+ }
+
+ LIBC_INLINE constexpr Frac64 &operator+=(const Frac64 &other) {
+ *this = *this + other;
+ return *this;
+ }
+
+ LIBC_INLINE constexpr Frac64 &operator-=(const Frac64 &other) {
+ *this = *this - other;
+ return *this;
+ }
+
+ LIBC_INLINE constexpr Frac64 &operator*=(const Frac64 &other) {
+ *this = *this * other;
+ return *this;
+ }
+};
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_FRAC64_H
diff --git a/libc/src/__support/math/exp_integer_utils.h b/libc/src/__support/math/exp_integer_utils.h
index d6216b6ee183f..4cd1d32511106 100644
--- a/libc/src/__support/math/exp_integer_utils.h
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -16,7 +16,7 @@
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/big_int.h"
-#include "src/__support/frac128.h"
+#include "src/__support/frac64.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/math_extras.h"
@@ -27,13 +27,33 @@
#define LIBC_TARGET_IS_BIG_ENDIAN 0
#else
#define LIBC_TARGET_IS_BIG_ENDIAN (__BYTE_ORDER__ == __ORDER_BIG_ENDIAN__)
-#endif // /LIBC_TARGET_IS_BIG_ENDIAN
+#endif // LIBC_TARGET_IS_BIG_ENDIAN
namespace LIBC_NAMESPACE_DECL {
namespace math {
-namespace integer_only {} // namespace integer_only
+namespace integer_only {
+
+// round(1/log(2), D, RN)
+LIBC_INLINE_VAR constexpr double INV_LN2 = 0x1.71547652b82fep+0;
+
+// TODO: Remez eval
+// TODO: test around and see what at what degree the accuracy is "acceptable"
+LIBC_INLINE_VAR constexpr Frac64 EXPF_COEFFS[6] = {};
+
+LIBC_INLINE constexpr float
+expf_range_reduction([[maybe_unused]] const float &x) {
+ // TODO: declare the range reduction params and this function
+ return 0.0f;
+}
+
+LIBC_INLINE constexpr float expf_eval() {
+ // TODO: declare params and this function
+ return 0.0f;
+}
+
+} // namespace integer_only
} // namespace math
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 1bab06719d0c8..b7ee31cf6856c 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -34,35 +34,42 @@ LIBC_INLINE float expf(float x) {
// Exceptional values
// TODO: optimize the branching order
- if (xbits.is_zero()) {
+ if (LIBC_UNLIKELY(xbits.is_zero())) {
// e^0 = 1 (exact), for both +/-0
return 1.0f;
- } else {
- // x is inf or NaN
- if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
- // e^NaN = NaN
- if (xbits.is_signaling_nan()) {
- // Per conversation with lntue, we don't need to raise exception here,
- // as we're assuming no FPUs/fenv in this kind of environment
-
- // silencing
- return FPBits::quiet_nan().get_val();
- }
-
- // e^-inf = 0
- // e^+inf = +inf
- if (xbits.is_inf()) {
- return is_neg ? 0.0f : FPBits::inf().get_val();
- }
-
- // x is a quiet NaN
- return x;
- } else {
- // TODO: out-of-range (overflow/underflow), exeecute range reduction
+ }
+
+ // x is inf or NaN
+ if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
+ // e^NaN = NaN
+ if (xbits.is_signaling_nan()) {
+ // Per conversation with lntue, we don't need to raise exception here,
+ // as we're assuming no FPUs/fenv in this kind of environment
+
+ // silencing
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // e^-inf = 0
+ // e^+inf = +inf
+ if (xbits.is_inf()) {
+ return is_neg ? 0.0f : FPBits::inf().get_val();
}
+
+ // x is a quiet NaN
+ return x;
}
+ // strats:
+ // 1. u32 fast path + u64 accurate path
+ // 2. dropping u32 altogether, just single u64 accurate path
+ // aim for minimal code size
+ // (speed/acc still important)
+
// TODO: execute the normal exp here
+ // TODO: evaluate ULPs, relative errors, etc.
+ // We're doing round-to-nearest only in this pass
+
return 0.0f;
}
>From ef44e3809c42f0038172bc2bdd67dc5f0b996590 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Thu, 23 Jul 2026 19:43:31 +0700
Subject: [PATCH 3/5] feat: basic algo
still failing on edge cases with too large errors
---
libc/src/__support/frac64.h | 30 +++--
libc/src/__support/math/exp_integer_utils.h | 36 +++---
libc/src/__support/math/expf_integer_eval.h | 120 +++++++++++++++-----
3 files changed, 126 insertions(+), 60 deletions(-)
diff --git a/libc/src/__support/frac64.h b/libc/src/__support/frac64.h
index 2172f4cbaeeb8..5763055e9b64c 100644
--- a/libc/src/__support/frac64.h
+++ b/libc/src/__support/frac64.h
@@ -10,6 +10,7 @@
#define LLVM_LIBC_SRC___SUPPORT_FRAC64_H
#include "hdr/stdint_proxy.h" // uint64_t
+#include "src/__support/big_int.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_INLINE
#include "src/__support/uint128.h"
@@ -17,38 +18,35 @@
namespace LIBC_NAMESPACE_DECL {
// Q0.64
-struct Frac64 {
- uint64_t val;
+struct Frac64 : public UInt<64> {
+ using UInt<64>::UInt;
- LIBC_INLINE constexpr Frac64() : val(0) {}
- LIBC_INLINE constexpr explicit Frac64(uint64_t value) : val(value) {}
+ LIBC_INLINE constexpr Frac64 operator~() const { return Frac64(~val[0]); }
- LIBC_INLINE constexpr Frac64 operator~() const { return Frac64(~val); }
-
- LIBC_INLINE constexpr Frac64 operator+(const Frac64 &other) const {
- return Frac64(val + other.val);
+ LIBC_INLINE constexpr Frac64 operator+(Frac64 other) const {
+ return Frac64(val[0] + other.val[0]);
}
- LIBC_INLINE constexpr Frac64 operator-(const Frac64 &other) const {
- return Frac64(val - other.val);
+ LIBC_INLINE constexpr Frac64 operator-(Frac64 other) const {
+ return Frac64(val[0] - other.val[0]);
}
- LIBC_INLINE constexpr Frac64 operator*(const Frac64 &other) const {
- return Frac64(
- static_cast<uint64_t>((static_cast<UInt128>(val) * other.val) >> 64));
+ LIBC_INLINE constexpr Frac64 operator*(Frac64 other) const {
+ UInt<64> r = UInt<64>::quick_mul_hi(UInt<64>(other));
+ return Frac64(r.val[0]);
}
- LIBC_INLINE constexpr Frac64 &operator+=(const Frac64 &other) {
+ LIBC_INLINE constexpr Frac64 &operator+=(Frac64 other) {
*this = *this + other;
return *this;
}
- LIBC_INLINE constexpr Frac64 &operator-=(const Frac64 &other) {
+ LIBC_INLINE constexpr Frac64 &operator-=(Frac64 other) {
*this = *this - other;
return *this;
}
- LIBC_INLINE constexpr Frac64 &operator*=(const Frac64 &other) {
+ LIBC_INLINE constexpr Frac64 &operator*=(Frac64 other) {
*this = *this * other;
return *this;
}
diff --git a/libc/src/__support/math/exp_integer_utils.h b/libc/src/__support/math/exp_integer_utils.h
index 4cd1d32511106..97688c654213f 100644
--- a/libc/src/__support/math/exp_integer_utils.h
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -35,23 +35,25 @@ namespace math {
namespace integer_only {
-// round(1/log(2), D, RN)
-LIBC_INLINE_VAR constexpr double INV_LN2 = 0x1.71547652b82fep+0;
-
-// TODO: Remez eval
-// TODO: test around and see what at what degree the accuracy is "acceptable"
-LIBC_INLINE_VAR constexpr Frac64 EXPF_COEFFS[6] = {};
-
-LIBC_INLINE constexpr float
-expf_range_reduction([[maybe_unused]] const float &x) {
- // TODO: declare the range reduction params and this function
- return 0.0f;
-}
-
-LIBC_INLINE constexpr float expf_eval() {
- // TODO: declare params and this function
- return 0.0f;
-}
+// print(2+round(1/log(2), 64, RN));
+// LSB(INV_LN2) = 2^-63
+LIBC_INLINE_VAR constexpr Frac64 INV_LN2 = Frac64(0xb8aa'3b29'5c17'f0bc);
+
+// 1-ULP
+// 64-bit polynomial approximation of 2^x coefficients generated with Sollya:
+// > P = fpminimax(2^x, 6, [|1, 64...|], [0, 1], absolute, fixed);
+// Store the fractional part of the coefficients below
+// > dirtyinfnorm(2^x - P(x), [0, 1]);
+// 0x1.9d39...p-29
+// ULPs of coeffs = 2^-64
+LIBC_INLINE_VAR constexpr Frac64 EXPF_COEFFS[] = {
+ Frac64(0xb172'14ea'215c'7750), // x
+ Frac64(0x3d7f'b5e7'4e78'9f2b), // x^2
+ Frac64(0x0e34'15ac'7481'5dee), // x^3
+ Frac64(0x027a'7e40'a2eb'6584), // x^4
+ Frac64(0x0051'56c0'9d53'15f3), // x^5
+ Frac64(0x000e'4a74'a170'46e8), // x^6
+};
} // namespace integer_only
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index b7ee31cf6856c..ec9380ede5cec 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -13,7 +13,8 @@
#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/frac128.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/frac64.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/math/exp_integer_utils.h"
@@ -29,25 +30,26 @@ LIBC_INLINE float expf(float x) {
FPBits xbits(x);
bool is_neg = xbits.is_neg();
- uint16_t x_e = xbits.get_biased_exponent();
- // uint64_t x_u = xbits.get_mantissa();
-
- // Exceptional values
- // TODO: optimize the branching order
- if (LIBC_UNLIKELY(xbits.is_zero())) {
- // e^0 = 1 (exact), for both +/-0
- return 1.0f;
- }
+ uint32_t x_bits_abs = xbits.uintval() & 0x7fff'ffffU;
+
+ // When |x| >= 89, |x| < 2^-25, or x is NaN
+ if (LIBC_UNLIKELY(x_bits_abs >= 0x42b2'0000U || x_bits_abs <= 0x3380'0000U)) {
+ // |x| < 2^-24
+ // e^x ~ 1 + x
+ if (xbits.get_biased_exponent() <= 103) {
+ return 1.0f + x;
+ }
- // x is inf or NaN
- if (LIBC_UNLIKELY(x_e > 2 * FPBits::EXP_BIAS)) {
- // e^NaN = NaN
- if (xbits.is_signaling_nan()) {
+ if (LIBC_UNLIKELY(xbits.is_nan())) {
// Per conversation with lntue, we don't need to raise exception here,
// as we're assuming no FPUs/fenv in this kind of environment
+ if (xbits.is_signaling_nan()) {
+ // silencing
+ return FPBits::quiet_nan().get_val();
+ }
- // silencing
- return FPBits::quiet_nan().get_val();
+ // quiet NaN
+ return x;
}
// e^-inf = 0
@@ -56,21 +58,85 @@ LIBC_INLINE float expf(float x) {
return is_neg ? 0.0f : FPBits::inf().get_val();
}
- // x is a quiet NaN
- return x;
+ // Large finite positive --> overflow
+ if (!is_neg) {
+ return FPBits::inf().get_val();
+ }
+
+ // x < log(2^-150) or NaN (NaN is already handled above)
+ if (xbits.uintval() >= 0xc2cf'f1b5U) {
+ return 0.0f;
+ }
}
- // strats:
- // 1. u32 fast path + u64 accurate path
- // 2. dropping u32 altogether, just single u64 accurate path
- // aim for minimal code size
- // (speed/acc still important)
+ // Main calculations
+
+ uint16_t x_e = xbits.get_biased_exponent();
+ uint64_t x_u = xbits.get_mantissa();
+
+ // Range reduction
+ // The algorithm near the end of this function estimates 2^r,
+ // where r is the fractional part of x * log2(e) and is in [0, 1].
+ // See EXPF_COEFFS for more details on the approximation polynomial used.
+
+ // add leading bit = 1
+ x_u |= uint64_t(1) << FPBits::FRACTION_LEN;
+
+ // shift to top 32 bit --> decimal point at hidden bit that we've added
+ x_u <<= 32;
+
+ int x_e_unbiased = static_cast<int>(x_e) - FPBits::EXP_BIAS;
+
+ // shift for the decimal point to be at the hidden bit
+ if (x_e_unbiased > 0) {
+ x_u <<= x_e_unbiased;
+ } else if (x_e_unbiased < 0) {
+ x_u >>= -x_e_unbiased;
+ }
+
+ // LSB(x_u_frac) = 2^-55
+ Frac64 x_u_frac(x_u);
+
+ // LSB(x_ln2) = 2^-54
+ Frac64 x_ln2 = x_u_frac * INV_LN2;
+
+ constexpr uint64_t FRAC_MASK = (uint64_t(1) << 54) - 1;
+ uint64_t x_ln2_bit = x_ln2.val[0];
+
+ uint64_t e_y, l2y_r;
+ uint32_t e_y_unbiased;
+
+ if (is_neg) {
+ e_y = (x_ln2_bit + FRAC_MASK) & ~FRAC_MASK;
+ l2y_r = e_y - x_ln2_bit;
+ e_y_unbiased = (FPBits::EXP_BIAS << 23) - static_cast<uint32_t>(e_y >> 31);
+ } else {
+ e_y = x_ln2_bit & ~FRAC_MASK;
+ l2y_r = x_ln2_bit - e_y;
+ e_y_unbiased = (FPBits::EXP_BIAS << 23) + static_cast<uint32_t>(e_y >> 31);
+ }
+
+ // LSB(l2y_r_frac) = LSB(l2y_r) * 2^-10 = 2^-64
+ Frac64 l2y_r_frac(l2y_r << 10);
+
+ // p = 2^l2y_r_frac - 1
+ Frac64 p = fputil::polyeval(l2y_r_frac, EXPF_COEFFS[0], EXPF_COEFFS[1],
+ EXPF_COEFFS[2], EXPF_COEFFS[3], EXPF_COEFFS[4],
+ EXPF_COEFFS[5]);
+ // TODO: handle denorm, underflow, overflow, boundary cases
+ // We're currently focusing on rounding to nearest only
+
+ // RN 23 bits --> shift for the LSB to be 2^-24 --> +1, shift for another
+ // bit
+ // Can be rounded up
+ // + e_y_unbiased
+
+ uint32_t result = (static_cast<uint32_t>(p.val[0] >> 40) + 1);
+ result >>= 1;
- // TODO: execute the normal exp here
- // TODO: evaluate ULPs, relative errors, etc.
- // We're doing round-to-nearest only in this pass
+ result += e_y_unbiased;
- return 0.0f;
+ return cpp::bit_cast<float>(result);
}
} // namespace integer_only
>From 0f1f82f4669501396f826e9a88f04e6e73d89466 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Thu, 23 Jul 2026 23:13:59 +0700
Subject: [PATCH 4/5] feat: clean up CMakeLists.txt deps
---
libc/src/__support/CMakeLists.txt | 9 +++++++++
libc/src/__support/math/CMakeLists.txt | 8 ++------
2 files changed, 11 insertions(+), 6 deletions(-)
diff --git a/libc/src/__support/CMakeLists.txt b/libc/src/__support/CMakeLists.txt
index 12761e6a524b9..3b3ca4eabd9c4 100644
--- a/libc/src/__support/CMakeLists.txt
+++ b/libc/src/__support/CMakeLists.txt
@@ -375,6 +375,15 @@ add_header_library(
libc.src.__support.macros.config
)
+add_header_library(
+ frac64
+ HDRS
+ frac64.h
+ DEPENDS
+ .big_int
+ libc.src.__support.macros.config
+)
+
add_header_library(
uint128
HDRS
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index b36cc1f537c69..25ea9ba82e41e 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3065,15 +3065,13 @@ add_header_library(
HDRS
expf_integer_eval.h
DEPENDS
- # TODO: strip the below deps down if able to
.exp_constants
.exp_integer_utils
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.fp_bits
- libc.src.__support.FPUtil.fenv_impl
libc.src.__support.macros.config
libc.src.__support.macros.optimization
- libc.src.__support.frac128
+ libc.src.__support.frac64
)
add_header_library(
@@ -4120,14 +4118,12 @@ add_header_library(
HDRS
exp_integer_utils.h
DEPENDS
- # TODO: strip the below deps down if able to
.exp_constants
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.fp_bits
- libc.src.__support.FPUtil.fenv_impl
libc.src.__support.macros.config
libc.src.__support.macros.optimization
- libc.src.__support.frac128
+ libc.src.__support.frac64
)
add_header_library(
>From 8c98338a8ad63e2d6bffe6c65b6da5d17eb867f6 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Fri, 24 Jul 2026 16:06:31 +0700
Subject: [PATCH 5/5] fix: inaccuracies and subnormal results handling
---
libc/src/__support/frac64.h | 3 --
libc/src/__support/math/exp_integer_utils.h | 4 +-
libc/src/__support/math/expf_integer_eval.h | 54 +++++++++++++++------
3 files changed, 41 insertions(+), 20 deletions(-)
diff --git a/libc/src/__support/frac64.h b/libc/src/__support/frac64.h
index 5763055e9b64c..89296b67fa61a 100644
--- a/libc/src/__support/frac64.h
+++ b/libc/src/__support/frac64.h
@@ -9,11 +9,8 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_FRAC64_H
#define LLVM_LIBC_SRC___SUPPORT_FRAC64_H
-#include "hdr/stdint_proxy.h" // uint64_t
#include "src/__support/big_int.h"
#include "src/__support/macros/config.h"
-#include "src/__support/macros/optimization.h" // LIBC_INLINE
-#include "src/__support/uint128.h"
namespace LIBC_NAMESPACE_DECL {
diff --git a/libc/src/__support/math/exp_integer_utils.h b/libc/src/__support/math/exp_integer_utils.h
index 97688c654213f..3fda41c6350c4 100644
--- a/libc/src/__support/math/exp_integer_utils.h
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -1,4 +1,4 @@
-//===-- e^x range reduction and evaluation using integer-only --*- C++ -*-===//
+//===-- e^x integer-only utility functions ----------------------*- C++ -*-===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
@@ -6,8 +6,6 @@
//
//===----------------------------------------------------------------------===//
-// TODO: update the description of this file above
-
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_UTILS_H
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index ec9380ede5cec..b4ad9ada54fa1 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -9,14 +9,13 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
-// TODO: clean up includes
#include "src/__support/CPP/bit.h"
-#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/frac64.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
+#include "src/__support/math/check/exp_exceptions.h"
#include "src/__support/math/exp_integer_utils.h"
namespace LIBC_NAMESPACE_DECL {
@@ -25,22 +24,27 @@ namespace math {
namespace integer_only {
+// Round-nearest, no except implementation of expf using integer-only
+// arithmetic.
LIBC_INLINE float expf(float x) {
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
bool is_neg = xbits.is_neg();
- uint32_t x_bits_abs = xbits.uintval() & 0x7fff'ffffU;
+ uint32_t x_val = xbits.uintval();
+ uint32_t x_val_abs = x_val & 0x7fff'ffffU;
// When |x| >= 89, |x| < 2^-25, or x is NaN
- if (LIBC_UNLIKELY(x_bits_abs >= 0x42b2'0000U || x_bits_abs <= 0x3380'0000U)) {
- // |x| < 2^-24
- // e^x ~ 1 + x
- if (xbits.get_biased_exponent() <= 103) {
+ if (LIBC_UNLIKELY(x_val_abs >= 0x42b2'0000U || x_val_abs <= 0x3380'0000U)) {
+ if (x_val_abs < 0x3300'0000U) { // |x| < 2^-25
+ return 1.0f;
+ }
+
+ if (x_val_abs < 0x3380'0000U) { // |x| < 2^-24
return 1.0f + x;
}
- if (LIBC_UNLIKELY(xbits.is_nan())) {
+ if (xbits.is_nan()) {
// Per conversation with lntue, we don't need to raise exception here,
// as we're assuming no FPUs/fenv in this kind of environment
if (xbits.is_signaling_nan()) {
@@ -69,6 +73,11 @@ LIBC_INLINE float expf(float x) {
}
}
+ if (LIBC_UNLIKELY(x_val >= check::exp_internal::Bounds<float>::UPPER_BITS &&
+ !is_neg)) { // overflow
+ return FPBits::inf().get_val();
+ }
+
// Main calculations
uint16_t x_e = xbits.get_biased_exponent();
@@ -106,7 +115,7 @@ LIBC_INLINE float expf(float x) {
uint64_t e_y, l2y_r;
uint32_t e_y_unbiased;
- if (is_neg) {
+ if (LIBC_UNLIKELY(is_neg)) {
e_y = (x_ln2_bit + FRAC_MASK) & ~FRAC_MASK;
l2y_r = e_y - x_ln2_bit;
e_y_unbiased = (FPBits::EXP_BIAS << 23) - static_cast<uint32_t>(e_y >> 31);
@@ -120,11 +129,28 @@ LIBC_INLINE float expf(float x) {
Frac64 l2y_r_frac(l2y_r << 10);
// p = 2^l2y_r_frac - 1
- Frac64 p = fputil::polyeval(l2y_r_frac, EXPF_COEFFS[0], EXPF_COEFFS[1],
- EXPF_COEFFS[2], EXPF_COEFFS[3], EXPF_COEFFS[4],
- EXPF_COEFFS[5]);
- // TODO: handle denorm, underflow, overflow, boundary cases
- // We're currently focusing on rounding to nearest only
+ Frac64 p = fputil::polyeval(l2y_r_frac, Frac64(0), EXPF_COEFFS[0],
+ EXPF_COEFFS[1], EXPF_COEFFS[2], EXPF_COEFFS[3],
+ EXPF_COEFFS[4], EXPF_COEFFS[5]);
+
+ uint32_t k = static_cast<uint32_t>(e_y >> 54);
+ int d = static_cast<int>(k) - FPBits::EXP_BIAS;
+
+ if (LIBC_UNLIKELY(is_neg && d >= 23)) { // underflow
+ return 0.0f;
+ }
+
+ if (LIBC_UNLIKELY(is_neg && d >= 0)) { // subnormal
+ uint64_t full_val = (uint64_t(1) << 63) | (p.val[0] >> 1);
+
+ // add rounding bit
+ full_val += (uint64_t(1) << (40 + d));
+
+ // shift back to align to 32-bit float representation
+ uint32_t result = static_cast<uint32_t>(full_val >> (41 + d));
+
+ return cpp::bit_cast<float>(result);
+ }
// RN 23 bits --> shift for the LSB to be 2^-24 --> +1, shift for another
// bit
More information about the libc-commits
mailing list