[libc-commits] [libc] [libc][math] Integer-only, statically rounded implementation of expf (PR #209406)
Hoàng Minh Thiên via libc-commits
libc-commits at lists.llvm.org
Mon Aug 10 06:05:13 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 01/17] 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 02/17] 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 03/17] 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 04/17] 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 05/17] 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
>From eb9a3ddf41c67c43a66971cbe7e6148df6e4b581 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Mon, 27 Jul 2026 23:39:05 +0700
Subject: [PATCH 06/17] feat: Frac32-based pipeline
but with 2-ULP errors, fixes pending
---
libc/src/__support/CMakeLists.txt | 9 +
libc/src/__support/frac32.h | 54 ++++++
libc/src/__support/math/exp_integer_utils.h | 19 ++
libc/src/__support/math/expf_integer_eval.h | 200 +++++++++++++++++---
4 files changed, 254 insertions(+), 28 deletions(-)
create mode 100644 libc/src/__support/frac32.h
diff --git a/libc/src/__support/CMakeLists.txt b/libc/src/__support/CMakeLists.txt
index 3b3ca4eabd9c4..28be2be4c29c0 100644
--- a/libc/src/__support/CMakeLists.txt
+++ b/libc/src/__support/CMakeLists.txt
@@ -384,6 +384,15 @@ add_header_library(
libc.src.__support.macros.config
)
+add_header_library(
+ frac32
+ HDRS
+ frac32.h
+ DEPENDS
+ .big_int
+ libc.src.__support.macros.config
+)
+
add_header_library(
uint128
HDRS
diff --git a/libc/src/__support/frac32.h b/libc/src/__support/frac32.h
new file mode 100644
index 0000000000000..176a29a83bd45
--- /dev/null
+++ b/libc/src/__support/frac32.h
@@ -0,0 +1,54 @@
+//===-- 32-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_FRAC32_H
+#define LLVM_LIBC_SRC___SUPPORT_FRAC32_H
+
+#include "src/__support/big_int.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+// Q0.32
+struct Frac32 : public UInt<32> {
+ using UInt<32>::UInt;
+
+ LIBC_INLINE constexpr Frac32 operator~() const { return Frac32(~val[0]); }
+
+ LIBC_INLINE constexpr Frac32 operator+(Frac32 other) const {
+ return Frac32(val[0] + other.val[0]);
+ }
+
+ LIBC_INLINE constexpr Frac32 operator-(Frac32 other) const {
+ return Frac32(val[0] - other.val[0]);
+ }
+
+ LIBC_INLINE constexpr Frac32 operator*(Frac32 other) const {
+ UInt<32> r = UInt<32>::quick_mul_hi(UInt<32>(other));
+ return Frac32(r.val[0]);
+ }
+
+ LIBC_INLINE constexpr Frac32 &operator+=(Frac32 other) {
+ *this = *this + other;
+ return *this;
+ }
+
+ LIBC_INLINE constexpr Frac32 &operator-=(Frac32 other) {
+ *this = *this - other;
+ return *this;
+ }
+
+ LIBC_INLINE constexpr Frac32 &operator*=(Frac32 other) {
+ *this = *this * other;
+ return *this;
+ }
+};
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_FRAC32_H
diff --git a/libc/src/__support/math/exp_integer_utils.h b/libc/src/__support/math/exp_integer_utils.h
index 3fda41c6350c4..87d2ccbd1eab5 100644
--- a/libc/src/__support/math/exp_integer_utils.h
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -14,7 +14,9 @@
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/big_int.h"
+#include "src/__support/frac32.h"
#include "src/__support/frac64.h"
+#include "src/__support/macros/attributes.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/math_extras.h"
@@ -53,6 +55,23 @@ LIBC_INLINE_VAR constexpr Frac64 EXPF_COEFFS[] = {
Frac64(0x000e'4a74'a170'46e8), // x^6
};
+// print(2+round(1/log(2), 32, RN));
+// LSB(INV_LN2_FRAC32) = 2^-31
+LIBC_INLINE_VAR constexpr Frac32 INV_LN2_FRAC32 = Frac32(0xb8aa'3b29);
+
+// 1-ULP
+// Degree-6 still works fine (p-29), degree-5 yields p-23 accuracy
+// This is still the case for Frac32!
+// 32-bit polynomial approximation of 2^x coefficients generated with Sollya:
+// > P = fpminimax(2^x, 6, [|1, 32...|], [0, 1], absolute, fixed);
+// Store the fractional part of the coefficients below
+// > dirtyinfnorm(2^x - P(x), [0, 1]);
+// 0x1.9ded...p-29
+// ULPs of coeffs = 2^-32
+LIBC_INLINE_VAR constexpr Frac32 EXPF_COEFFS_FRAC32[] = {
+ Frac32(0xb172'14e8), Frac32(0x3d7f'b5f8), Frac32(0x0e34'1554),
+ Frac32(0x027a'7f04), Frac32(0x0051'55fe), Frac32(0x000e'4abc)};
+
} // 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 b4ad9ada54fa1..28d7b934fe251 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -12,7 +12,7 @@
#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
-#include "src/__support/frac64.h"
+#include "src/__support/frac32.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/math/check/exp_exceptions.h"
@@ -26,6 +26,152 @@ 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_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_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 (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();
+// }
+
+// // quiet NaN
+// return x;
+// }
+
+// // e^-inf = 0
+// // e^+inf = +inf
+// if (xbits.is_inf()) {
+// return is_neg ? 0.0f : FPBits::inf().get_val();
+// }
+
+// // 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;
+// }
+// }
+
+// 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();
+// 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 (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);
+// } 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, 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
+// // Can be rounded up
+// // + e_y_unbiased
+
+// uint32_t result = (static_cast<uint32_t>(p.val[0] >> 40) + 1);
+// result >>= 1;
+
+// result += e_y_unbiased;
+
+// return cpp::bit_cast<float>(result);
+// }
+
+// Round-nearest, no except implementation of expf using integer-only
+// arithmetic.
+// Frac32 implementation
LIBC_INLINE float expf(float x) {
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
@@ -81,18 +227,16 @@ LIBC_INLINE float expf(float x) {
// Main calculations
uint16_t x_e = xbits.get_biased_exponent();
- uint64_t x_u = xbits.get_mantissa();
+ uint32_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.
+ // See EXPF_COEFFS_FRAC32 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;
+ x_u |= uint32_t(1) << FPBits::FRACTION_LEN;
int x_e_unbiased = static_cast<int>(x_e) - FPBits::EXP_BIAS;
@@ -103,37 +247,38 @@ LIBC_INLINE float expf(float x) {
x_u >>= -x_e_unbiased;
}
- // LSB(x_u_frac) = 2^-55
- Frac64 x_u_frac(x_u);
+ // LSB(x_u_frac) = 2^-23
+ Frac32 x_u_frac(x_u);
- // LSB(x_ln2) = 2^-54
- Frac64 x_ln2 = x_u_frac * INV_LN2;
+ // LSB(x_ln2) = 2^-22
+ Frac32 x_ln2 = x_u_frac * INV_LN2_FRAC32;
- constexpr uint64_t FRAC_MASK = (uint64_t(1) << 54) - 1;
- uint64_t x_ln2_bit = x_ln2.val[0];
+ constexpr uint32_t FRAC_MASK = (uint32_t(1) << 22) - 1;
+ uint32_t x_ln2_bit = x_ln2.val[0];
- uint64_t e_y, l2y_r;
+ uint32_t e_y, l2y_r;
uint32_t e_y_unbiased;
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);
+ e_y_unbiased = (FPBits::EXP_BIAS << 23) - static_cast<uint32_t>(e_y << 1);
} 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);
+ e_y_unbiased = (FPBits::EXP_BIAS << 23) + static_cast<uint32_t>(e_y << 1);
}
- // LSB(l2y_r_frac) = LSB(l2y_r) * 2^-10 = 2^-64
- Frac64 l2y_r_frac(l2y_r << 10);
+ // LSB(l2y_r_frac) = LSB(l2y_r) * 2^-10 = 2^-32
+ Frac32 l2y_r_frac(l2y_r << 10);
// p = 2^l2y_r_frac - 1
- 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]);
+ Frac32 p = fputil::polyeval(l2y_r_frac, Frac32(0), EXPF_COEFFS_FRAC32[0],
+ EXPF_COEFFS_FRAC32[1], EXPF_COEFFS_FRAC32[2],
+ EXPF_COEFFS_FRAC32[3], EXPF_COEFFS_FRAC32[4],
+ EXPF_COEFFS_FRAC32[5]);
- uint32_t k = static_cast<uint32_t>(e_y >> 54);
+ uint32_t k = static_cast<uint32_t>(e_y >> 22);
int d = static_cast<int>(k) - FPBits::EXP_BIAS;
if (LIBC_UNLIKELY(is_neg && d >= 23)) { // underflow
@@ -141,23 +286,22 @@ LIBC_INLINE float expf(float x) {
}
if (LIBC_UNLIKELY(is_neg && d >= 0)) { // subnormal
- uint64_t full_val = (uint64_t(1) << 63) | (p.val[0] >> 1);
+ uint64_t full_val = (uint64_t(1) << 32) | p.val[0]; // 1 + p
// add rounding bit
- full_val += (uint64_t(1) << (40 + d));
+ full_val += (uint64_t(1) << (9 + d));
// shift back to align to 32-bit float representation
- uint32_t result = static_cast<uint32_t>(full_val >> (41 + d));
+ uint32_t result = static_cast<uint32_t>(full_val >> (10 + d));
return cpp::bit_cast<float>(result);
}
- // RN 23 bits --> shift for the LSB to be 2^-24 --> +1, shift for another
- // bit
+ // 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);
+ uint32_t result = (static_cast<uint32_t>(p.val[0] >> 8) + 1);
result >>= 1;
result += e_y_unbiased;
>From 230e8acf10f294ae1425da7c6e43611ee37889fa Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Mon, 27 Jul 2026 23:51:05 +0700
Subject: [PATCH 07/17] chore: fix formatting
---
libc/src/__support/math/expf_integer_eval.h | 3 ++-
1 file changed, 2 insertions(+), 1 deletion(-)
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 28d7b934fe251..f1e9c1143e26b 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -286,7 +286,8 @@ LIBC_INLINE float expf(float x) {
}
if (LIBC_UNLIKELY(is_neg && d >= 0)) { // subnormal
- uint64_t full_val = (uint64_t(1) << 32) | p.val[0]; // 1 + p
+ // 1 + p
+ uint64_t full_val = (uint64_t(1) << 32) | p.val[0];
// add rounding bit
full_val += (uint64_t(1) << (9 + d));
>From 96ff74e6ce3be5eddfa970baac788a96d4edeb26 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Tue, 28 Jul 2026 18:43:20 +0700
Subject: [PATCH 08/17] chore: update file headers to conform with the new
template
---
libc/src/__support/frac32.h | 7 ++++++-
libc/src/__support/frac64.h | 8 ++++++--
libc/src/__support/math/exp_integer_utils.h | 11 ++++++++++-
libc/src/__support/math/expf_integer_eval.h | 8 +++++++-
4 files changed, 29 insertions(+), 5 deletions(-)
diff --git a/libc/src/__support/frac32.h b/libc/src/__support/frac32.h
index 176a29a83bd45..9e224e6a5008f 100644
--- a/libc/src/__support/frac32.h
+++ b/libc/src/__support/frac32.h
@@ -1,10 +1,15 @@
-//===-- 32-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
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the declaration of 32-bit unsigned fractional type
+///
+//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_FRAC32_H
#define LLVM_LIBC_SRC___SUPPORT_FRAC32_H
diff --git a/libc/src/__support/frac64.h b/libc/src/__support/frac64.h
index 89296b67fa61a..bbb0d3bb30903 100644
--- a/libc/src/__support/frac64.h
+++ b/libc/src/__support/frac64.h
@@ -1,10 +1,15 @@
-//===-- 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
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the declaration of 64-bit unsigned fractional type
+///
+//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_FRAC64_H
#define LLVM_LIBC_SRC___SUPPORT_FRAC64_H
@@ -14,7 +19,6 @@
namespace LIBC_NAMESPACE_DECL {
-// Q0.64
struct Frac64 : public UInt<64> {
using UInt<64>::UInt;
diff --git a/libc/src/__support/math/exp_integer_utils.h b/libc/src/__support/math/exp_integer_utils.h
index 87d2ccbd1eab5..aa3be0875ca24 100644
--- a/libc/src/__support/math/exp_integer_utils.h
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -1,10 +1,19 @@
-//===-- 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.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains utilities for implementing integer-only implementation
+/// of expf(x)
+///
+//===----------------------------------------------------------------------===//
+
+// TODO: this file is expected to be used for exp*(x) functions, not only
+// limited to expf(x)
#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 f1e9c1143e26b..77a56aa8fec64 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -1,10 +1,16 @@
-//===-- 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
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the integer-only implementation of expf(x)
+///
+//===----------------------------------------------------------------------===//
+
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
>From 9ee2810f72452b797413b44b4c7726f73e33a982 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Tue, 28 Jul 2026 18:45:34 +0700
Subject: [PATCH 09/17] chore: fix formatting
---
libc/src/__support/math/expf_integer_eval.h | 1 -
1 file changed, 1 deletion(-)
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 77a56aa8fec64..6e966f2203e58 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -11,7 +11,6 @@
///
//===----------------------------------------------------------------------===//
-
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
>From c3585195de8f66d866a731f4e97adcecc8e5b0c4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Sun, 9 Aug 2026 16:23:16 +0700
Subject: [PATCH 10/17] fix: ~0.5-ULP 64-bit
alongside cleanups
---
libc/src/__support/math/exp_integer_utils.h | 25 ++-
libc/src/__support/math/expf_integer_eval.h | 229 ++++----------------
2 files changed, 64 insertions(+), 190 deletions(-)
diff --git a/libc/src/__support/math/exp_integer_utils.h b/libc/src/__support/math/exp_integer_utils.h
index aa3be0875ca24..005e977f97ce3 100644
--- a/libc/src/__support/math/exp_integer_utils.h
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -55,13 +55,26 @@ LIBC_INLINE_VAR constexpr Frac64 INV_LN2 = Frac64(0xb8aa'3b29'5c17'f0bc);
// > 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
+// };
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
+ Frac64(0xb172'17f7'd1cf'b7cf), // x
+ Frac64(0x3d7f'7bff'057d'4a5e), // x^2
+ Frac64(0x0e35'846b'8363'9484), // x^3
+ Frac64(0x0276'556d'ec97'dcd4), // x^4
+ Frac64(0x0057'61ff'dc04'c7ff), // x^5
+ Frac64(0x000a'1847'b6e7'92ec), // x^6
+ Frac64(0x0000'ffe8'14e5'7033), // x^7
+ Frac64(0x0000'1628'b6e9'70c8), // x^8
+ Frac64(0x0000'01b8'8ce7'4088), // x^9
+ Frac64(0x0000'001c'18d5'cb29), // x^10
+ Frac64(0x0000'0002'b43f'4490), // x^11
};
// print(2+round(1/log(2), 32, RN));
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 6e966f2203e58..98dc16a6c5836 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -22,6 +22,8 @@
#include "src/__support/macros/optimization.h"
#include "src/__support/math/check/exp_exceptions.h"
#include "src/__support/math/exp_integer_utils.h"
+#include "src/__support/uint128.h"
+#include "src/string/memory_utils/aarch64/inline_memset.h"
namespace LIBC_NAMESPACE_DECL {
@@ -31,152 +33,10 @@ 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_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_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 (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();
-// }
-
-// // quiet NaN
-// return x;
-// }
-
-// // e^-inf = 0
-// // e^+inf = +inf
-// if (xbits.is_inf()) {
-// return is_neg ? 0.0f : FPBits::inf().get_val();
-// }
-
-// // 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;
-// }
-// }
-
-// 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();
-// 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 (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);
-// } 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, 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
-// // Can be rounded up
-// // + e_y_unbiased
-
-// uint32_t result = (static_cast<uint32_t>(p.val[0] >> 40) + 1);
-// result >>= 1;
-
-// result += e_y_unbiased;
-
-// return cpp::bit_cast<float>(result);
-// }
-
-// Round-nearest, no except implementation of expf using integer-only
-// arithmetic.
-// Frac32 implementation
+// Failing exhaustive tests:
+// Negative:
+// Test failed for 90853 inputs in range: 3267362816 to 3268411392 [0xc2c00000,
+// 0xc2d00000), [-0x1.8p+6, -0x1.ap+6)
LIBC_INLINE float expf(float x) {
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
@@ -185,16 +45,13 @@ LIBC_INLINE float expf(float x) {
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_val_abs >= 0x42b2'0000U || x_val_abs <= 0x3380'0000U)) {
- if (x_val_abs < 0x3300'0000U) { // |x| < 2^-25
+ // When |x| >= 89, |x| <= 2^-25, or x is NaN
+ if (LIBC_UNLIKELY(x_val_abs >= 0x42b1'7218 || x_val_abs <= 0x3300'0000U)) {
+ // |x| <= 2^-25
+ if (xbits.get_biased_exponent() <= 102) {
return 1.0f;
}
- if (x_val_abs < 0x3380'0000U) { // |x| < 2^-24
- return 1.0f + x;
- }
-
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
@@ -224,24 +81,21 @@ 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();
- uint32_t x_u = xbits.get_mantissa();
+ 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_FRAC32 for more details on the approximation polynomial
- // used.
+ // See EXPF_COEFFS for more details on the approximation polynomial used.
// add leading bit = 1
- x_u |= uint32_t(1) << FPBits::FRACTION_LEN;
+ 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;
@@ -252,62 +106,69 @@ LIBC_INLINE float expf(float x) {
x_u >>= -x_e_unbiased;
}
- // LSB(x_u_frac) = 2^-23
- Frac32 x_u_frac(x_u);
+ // LSB(x_u_frac) = 2^-55
+ Frac64 x_u_frac(x_u);
- // LSB(x_ln2) = 2^-22
- Frac32 x_ln2 = x_u_frac * INV_LN2_FRAC32;
+ // LSB(x_ln2) = 2^-54
+ Frac64 x_ln2 = x_u_frac * INV_LN2;
- constexpr uint32_t FRAC_MASK = (uint32_t(1) << 22) - 1;
- uint32_t x_ln2_bit = x_ln2.val[0];
+ constexpr uint64_t FRAC_MASK = (uint64_t(1) << 54) - 1;
+ uint64_t x_ln2_bit = x_ln2.val[0];
- uint32_t e_y, l2y_r;
+ uint64_t e_y, l2y_r;
uint32_t e_y_unbiased;
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 << 1);
+ 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 << 1);
+ e_y_unbiased = (FPBits::EXP_BIAS << 23) + static_cast<uint32_t>(e_y >> 31);
}
- // LSB(l2y_r_frac) = LSB(l2y_r) * 2^-10 = 2^-32
- Frac32 l2y_r_frac(l2y_r << 10);
+ // LSB(l2y_r_frac) = LSB(l2y_r) * 2^-10 = 2^-64
+ Frac64 l2y_r_frac(l2y_r << 10);
// p = 2^l2y_r_frac - 1
- Frac32 p = fputil::polyeval(l2y_r_frac, Frac32(0), EXPF_COEFFS_FRAC32[0],
- EXPF_COEFFS_FRAC32[1], EXPF_COEFFS_FRAC32[2],
- EXPF_COEFFS_FRAC32[3], EXPF_COEFFS_FRAC32[4],
- EXPF_COEFFS_FRAC32[5]);
+ Frac64 p = l2y_r_frac *
+ fputil::polyeval(l2y_r_frac, EXPF_COEFFS[0], EXPF_COEFFS[1],
+ EXPF_COEFFS[2], EXPF_COEFFS[3], EXPF_COEFFS[4],
+ EXPF_COEFFS[5], EXPF_COEFFS[6], EXPF_COEFFS[7],
+ EXPF_COEFFS[8], EXPF_COEFFS[9], EXPF_COEFFS[10]);
- uint32_t k = static_cast<uint32_t>(e_y >> 22);
+ 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;
}
+ // TODO: uint64_t will cause overflow, dialing to UInt128 works flawlessly,
+ // but wondering if there is any lighter solution without touching the 128-bit
+ // types?
if (LIBC_UNLIKELY(is_neg && d >= 0)) { // subnormal
- // 1 + p
- uint64_t full_val = (uint64_t(1) << 32) | p.val[0];
+ UInt128 full_val = (UInt128(1) << 63) | (p.val[0] >> 1); // 1 + p
// add rounding bit
- full_val += (uint64_t(1) << (9 + d));
+ // TODO: skip for R0, RD
+ full_val += (UInt128(1) << (40 + d));
+
+ // TODO: RU --> mask
// shift back to align to 32-bit float representation
- uint32_t result = static_cast<uint32_t>(full_val >> (10 + d));
+ 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
+ // 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] >> 8) + 1);
+ uint32_t result = (static_cast<uint32_t>(p.val[0] >> 40) + 1);
result >>= 1;
result += e_y_unbiased;
>From 2ebb466c9c3e700980020ab8364be68a115acae1 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Sun, 9 Aug 2026 16:34:34 +0700
Subject: [PATCH 11/17] chore: fix formatting
---
libc/src/__support/math/expf_integer_eval.h | 3 ++-
1 file changed, 2 insertions(+), 1 deletion(-)
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 98dc16a6c5836..6904622b6264c 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -149,7 +149,8 @@ LIBC_INLINE float expf(float x) {
// but wondering if there is any lighter solution without touching the 128-bit
// types?
if (LIBC_UNLIKELY(is_neg && d >= 0)) { // subnormal
- UInt128 full_val = (UInt128(1) << 63) | (p.val[0] >> 1); // 1 + p
+ // 1 + p
+ UInt128 full_val = (UInt128(1) << 63) | (p.val[0] >> 1);
// add rounding bit
// TODO: skip for R0, RD
>From d8dbc145cdc1846403153a3f4d174b33c6112c66 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Sun, 9 Aug 2026 18:13:38 +0700
Subject: [PATCH 12/17] fix: comment inaccuracies, and clean Frac32 left-overs
---
libc/src/__support/CMakeLists.txt | 9 ----
libc/src/__support/frac32.h | 59 ---------------------
libc/src/__support/math/exp_integer_utils.h | 41 ++------------
libc/src/__support/math/expf_integer_eval.h | 9 ++--
4 files changed, 8 insertions(+), 110 deletions(-)
delete mode 100644 libc/src/__support/frac32.h
diff --git a/libc/src/__support/CMakeLists.txt b/libc/src/__support/CMakeLists.txt
index 28be2be4c29c0..3b3ca4eabd9c4 100644
--- a/libc/src/__support/CMakeLists.txt
+++ b/libc/src/__support/CMakeLists.txt
@@ -384,15 +384,6 @@ add_header_library(
libc.src.__support.macros.config
)
-add_header_library(
- frac32
- HDRS
- frac32.h
- DEPENDS
- .big_int
- libc.src.__support.macros.config
-)
-
add_header_library(
uint128
HDRS
diff --git a/libc/src/__support/frac32.h b/libc/src/__support/frac32.h
deleted file mode 100644
index 9e224e6a5008f..0000000000000
--- a/libc/src/__support/frac32.h
+++ /dev/null
@@ -1,59 +0,0 @@
-//===----------------------------------------------------------------------===//
-//
-// 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
-//
-//===----------------------------------------------------------------------===//
-///
-/// \file
-/// This file contains the declaration of 32-bit unsigned fractional type
-///
-//===----------------------------------------------------------------------===//
-
-#ifndef LLVM_LIBC_SRC___SUPPORT_FRAC32_H
-#define LLVM_LIBC_SRC___SUPPORT_FRAC32_H
-
-#include "src/__support/big_int.h"
-#include "src/__support/macros/config.h"
-
-namespace LIBC_NAMESPACE_DECL {
-
-// Q0.32
-struct Frac32 : public UInt<32> {
- using UInt<32>::UInt;
-
- LIBC_INLINE constexpr Frac32 operator~() const { return Frac32(~val[0]); }
-
- LIBC_INLINE constexpr Frac32 operator+(Frac32 other) const {
- return Frac32(val[0] + other.val[0]);
- }
-
- LIBC_INLINE constexpr Frac32 operator-(Frac32 other) const {
- return Frac32(val[0] - other.val[0]);
- }
-
- LIBC_INLINE constexpr Frac32 operator*(Frac32 other) const {
- UInt<32> r = UInt<32>::quick_mul_hi(UInt<32>(other));
- return Frac32(r.val[0]);
- }
-
- LIBC_INLINE constexpr Frac32 &operator+=(Frac32 other) {
- *this = *this + other;
- return *this;
- }
-
- LIBC_INLINE constexpr Frac32 &operator-=(Frac32 other) {
- *this = *this - other;
- return *this;
- }
-
- LIBC_INLINE constexpr Frac32 &operator*=(Frac32 other) {
- *this = *this * other;
- return *this;
- }
-};
-
-} // namespace LIBC_NAMESPACE_DECL
-
-#endif // LLVM_LIBC_SRC___SUPPORT_FRAC32_H
diff --git a/libc/src/__support/math/exp_integer_utils.h b/libc/src/__support/math/exp_integer_utils.h
index 005e977f97ce3..0f9acf88cf28d 100644
--- a/libc/src/__support/math/exp_integer_utils.h
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -23,21 +23,12 @@
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/big_int.h"
-#include "src/__support/frac32.h"
#include "src/__support/frac64.h"
#include "src/__support/macros/attributes.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 {
@@ -48,21 +39,12 @@ namespace integer_only {
// 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);
+// > P = fpminimax(2^x, 11, [|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
-// };
+// 0x1.6238...p-58
+// LSB(EXPF_COEFFS[i]) = 2^-64
LIBC_INLINE_VAR constexpr Frac64 EXPF_COEFFS[] = {
Frac64(0xb172'17f7'd1cf'b7cf), // x
Frac64(0x3d7f'7bff'057d'4a5e), // x^2
@@ -77,23 +59,6 @@ LIBC_INLINE_VAR constexpr Frac64 EXPF_COEFFS[] = {
Frac64(0x0000'0002'b43f'4490), // x^11
};
-// print(2+round(1/log(2), 32, RN));
-// LSB(INV_LN2_FRAC32) = 2^-31
-LIBC_INLINE_VAR constexpr Frac32 INV_LN2_FRAC32 = Frac32(0xb8aa'3b29);
-
-// 1-ULP
-// Degree-6 still works fine (p-29), degree-5 yields p-23 accuracy
-// This is still the case for Frac32!
-// 32-bit polynomial approximation of 2^x coefficients generated with Sollya:
-// > P = fpminimax(2^x, 6, [|1, 32...|], [0, 1], absolute, fixed);
-// Store the fractional part of the coefficients below
-// > dirtyinfnorm(2^x - P(x), [0, 1]);
-// 0x1.9ded...p-29
-// ULPs of coeffs = 2^-32
-LIBC_INLINE_VAR constexpr Frac32 EXPF_COEFFS_FRAC32[] = {
- Frac32(0xb172'14e8), Frac32(0x3d7f'b5f8), Frac32(0x0e34'1554),
- Frac32(0x027a'7f04), Frac32(0x0051'55fe), Frac32(0x000e'4abc)};
-
} // 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 6904622b6264c..8317aae4a0192 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -17,13 +17,11 @@
#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
-#include "src/__support/frac32.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"
#include "src/__support/uint128.h"
-#include "src/string/memory_utils/aarch64/inline_memset.h"
namespace LIBC_NAMESPACE_DECL {
@@ -45,8 +43,11 @@ LIBC_INLINE float expf(float x) {
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_val_abs >= 0x42b1'7218 || x_val_abs <= 0x3300'0000U)) {
+ // When |x| >= smallest value that will cause overflow, |x| <= 2^-25, or x is
+ // NaN
+ if (LIBC_UNLIKELY(x_val_abs >=
+ check::exp_internal::Bounds<float>::UPPER_BITS ||
+ x_val_abs <= 0x3300'0000U)) {
// |x| <= 2^-25
if (xbits.get_biased_exponent() <= 102) {
return 1.0f;
>From 5277560e8a0871730479c4c64723f1cb7edcae1b Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Mon, 10 Aug 2026 10:49:25 +0700
Subject: [PATCH 13/17] feat: dropping LSBs to allow processing only on
uint64_t at max
---
libc/src/__support/math/expf_integer_eval.h | 14 +++++++-------
1 file changed, 7 insertions(+), 7 deletions(-)
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 8317aae4a0192..c5126aef5e130 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -21,7 +21,6 @@
#include "src/__support/macros/optimization.h"
#include "src/__support/math/check/exp_exceptions.h"
#include "src/__support/math/exp_integer_utils.h"
-#include "src/__support/uint128.h"
namespace LIBC_NAMESPACE_DECL {
@@ -146,21 +145,22 @@ LIBC_INLINE float expf(float x) {
return 0.0f;
}
- // TODO: uint64_t will cause overflow, dialing to UInt128 works flawlessly,
- // but wondering if there is any lighter solution without touching the 128-bit
- // types?
+ // Dropping some last bits (won't need them as we're casting into 32-bit float
+ // anyway)
+ constexpr uint32_t DROP_BITS = 2;
if (LIBC_UNLIKELY(is_neg && d >= 0)) { // subnormal
// 1 + p
- UInt128 full_val = (UInt128(1) << 63) | (p.val[0] >> 1);
+ uint64_t full_val =
+ (uint64_t(1) << (64 - DROP_BITS)) | (p.val[0] >> DROP_BITS);
// add rounding bit
// TODO: skip for R0, RD
- full_val += (UInt128(1) << (40 + d));
+ full_val += (uint64_t(1) << ((41 - DROP_BITS) + d));
// TODO: RU --> mask
// shift back to align to 32-bit float representation
- uint32_t result = static_cast<uint32_t>(full_val >> (41 + d));
+ uint32_t result = static_cast<uint32_t>(full_val >> ((42 - DROP_BITS) + d));
return cpp::bit_cast<float>(result);
}
>From 957a2e9a72db562a87e733f0f014c1cfe5a5faf4 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Mon, 10 Aug 2026 12:15:07 +0700
Subject: [PATCH 14/17] feat: migrate from shift-compare to a single compare op
---
libc/src/__support/math/expf_integer_eval.h | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index c5126aef5e130..35f783c6b2f8a 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -48,7 +48,7 @@ LIBC_INLINE float expf(float x) {
check::exp_internal::Bounds<float>::UPPER_BITS ||
x_val_abs <= 0x3300'0000U)) {
// |x| <= 2^-25
- if (xbits.get_biased_exponent() <= 102) {
+ if (x_val_abs <= 0x3300'0000U) {
return 1.0f;
}
>From 4f27e9ee02e6dcc6dad70f6594fe685e307c5456 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Mon, 10 Aug 2026 16:03:02 +0700
Subject: [PATCH 15/17] feat: no more failed exhaustive tests
---
libc/src/__support/math/expf_integer_eval.h | 20 +++++++++-----------
1 file changed, 9 insertions(+), 11 deletions(-)
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 35f783c6b2f8a..7757edf1283fc 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -30,10 +30,6 @@ namespace integer_only {
// Round-nearest, no except implementation of expf using integer-only
// arithmetic.
-// Failing exhaustive tests:
-// Negative:
-// Test failed for 90853 inputs in range: 3267362816 to 3268411392 [0xc2c00000,
-// 0xc2d00000), [-0x1.8p+6, -0x1.ap+6)
LIBC_INLINE float expf(float x) {
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
@@ -128,6 +124,15 @@ LIBC_INLINE float expf(float x) {
e_y_unbiased = (FPBits::EXP_BIAS << 23) + static_cast<uint32_t>(e_y >> 31);
}
+ uint32_t k = static_cast<uint32_t>(e_y >> 54);
+ int d = static_cast<int>(k) - FPBits::EXP_BIAS;
+
+ // d >= 24 --> k >= 151
+ // --> guaranteed to below 2^-150
+ if (LIBC_UNLIKELY(is_neg && d >= 24)) { // underflow
+ return 0.0f;
+ }
+
// LSB(l2y_r_frac) = LSB(l2y_r) * 2^-10 = 2^-64
Frac64 l2y_r_frac(l2y_r << 10);
@@ -138,13 +143,6 @@ LIBC_INLINE float expf(float x) {
EXPF_COEFFS[5], EXPF_COEFFS[6], EXPF_COEFFS[7],
EXPF_COEFFS[8], EXPF_COEFFS[9], EXPF_COEFFS[10]);
- 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;
- }
-
// Dropping some last bits (won't need them as we're casting into 32-bit float
// anyway)
constexpr uint32_t DROP_BITS = 2;
>From 5219fc2615940306629adddceecd2041f158d3a2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Mon, 10 Aug 2026 19:57:13 +0700
Subject: [PATCH 16/17] feat: statically rounded symbol, separate unit and
exhaustive tests
---
.../math/static_rounding/expf.h} | 72 ++++--
libc/shared/math/static_rounding_math.h | 30 +++
libc/src/__support/math/CMakeLists.txt | 27 --
libc/src/__support/math/exp_integer_utils.h | 68 -----
libc/src/math/generic/CMakeLists.txt | 1 -
libc/src/math/generic/expf.cpp | 4 +-
libc/test/src/math/CMakeLists.txt | 18 +-
libc/test/src/math/exhaustive/CMakeLists.txt | 25 ++
.../exhaustive_test_static_rounding.h | 244 ++++++++++++++++++
.../exhaustive/expf_static_rounding_test.cpp | 32 +++
.../src/math/expf_static_rounding_test.cpp | 148 +++++++++++
11 files changed, 551 insertions(+), 118 deletions(-)
rename libc/{src/__support/math/expf_integer_eval.h => shared/math/static_rounding/expf.h} (70%)
create mode 100644 libc/shared/math/static_rounding_math.h
delete mode 100644 libc/src/__support/math/exp_integer_utils.h
create mode 100644 libc/test/src/math/exhaustive/exhaustive_test_static_rounding.h
create mode 100644 libc/test/src/math/exhaustive/expf_static_rounding_test.cpp
create mode 100644 libc/test/src/math/expf_static_rounding_test.cpp
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/shared/math/static_rounding/expf.h
similarity index 70%
rename from libc/src/__support/math/expf_integer_eval.h
rename to libc/shared/math/static_rounding/expf.h
index 7757edf1283fc..7d69747cc37b9 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/shared/math/static_rounding/expf.h
@@ -7,30 +7,58 @@
//===----------------------------------------------------------------------===//
///
/// \file
-/// This file contains the integer-only implementation of expf(x)
+/// This file contains the statically rounded implementation of expf(x)
///
//===----------------------------------------------------------------------===//
-#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
-#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
+#ifndef LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_EXPF_H
+#define LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_EXPF_H
#include "src/__support/CPP/bit.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 {
-namespace math {
+using math::check::exp_internal::Bounds;
+
+namespace shared {
-namespace integer_only {
+namespace math {
-// Round-nearest, no except implementation of expf using integer-only
+namespace static_rounding {
+
+// 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);
+
+// 64-bit polynomial approximation of 2^x coefficients generated with Sollya:
+// > P = fpminimax(2^x, 11, [|1, 64...|], [0, 1], absolute, fixed);
+// Store the fractional part of the coefficients below
+// > dirtyinfnorm(2^x - P(x), [0, 1]);
+// 0x1.6238...p-58
+// LSB(EXPF_COEFFS[i]) = 2^-64
+LIBC_INLINE_VAR constexpr Frac64 EXPF_COEFFS[] = {
+ Frac64(0xb172'17f7'd1cf'b7cf), // x
+ Frac64(0x3d7f'7bff'057d'4a5e), // x^2
+ Frac64(0x0e35'846b'8363'9484), // x^3
+ Frac64(0x0276'556d'ec97'dcd4), // x^4
+ Frac64(0x0057'61ff'dc04'c7ff), // x^5
+ Frac64(0x000a'1847'b6e7'92ec), // x^6
+ Frac64(0x0000'ffe8'14e5'7033), // x^7
+ Frac64(0x0000'1628'b6e9'70c8), // x^8
+ Frac64(0x0000'01b8'8ce7'4088), // x^9
+ Frac64(0x0000'001c'18d5'cb29), // x^10
+ Frac64(0x0000'0002'b43f'4490), // x^11
+};
+
+// Statically rounded, no except implementation of expf using integer-only
// arithmetic.
-LIBC_INLINE float expf(float x) {
+LIBC_INLINE float expf(float x, int rounding) {
using FPBits = typename fputil::FPBits<float>;
FPBits xbits(x);
@@ -40,8 +68,7 @@ LIBC_INLINE float expf(float x) {
// When |x| >= smallest value that will cause overflow, |x| <= 2^-25, or x is
// NaN
- if (LIBC_UNLIKELY(x_val_abs >=
- check::exp_internal::Bounds<float>::UPPER_BITS ||
+ if (LIBC_UNLIKELY(x_val_abs >= Bounds<float>::UPPER_BITS ||
x_val_abs <= 0x3300'0000U)) {
// |x| <= 2^-25
if (x_val_abs <= 0x3300'0000U) {
@@ -151,11 +178,20 @@ LIBC_INLINE float expf(float x) {
uint64_t full_val =
(uint64_t(1) << (64 - DROP_BITS)) | (p.val[0] >> DROP_BITS);
- // add rounding bit
- // TODO: skip for R0, RD
- full_val += (uint64_t(1) << ((41 - DROP_BITS) + d));
-
- // TODO: RU --> mask
+// add rounding bit
+// skip for R0, RD
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (LIBC_UNLIKELY(rounding != FE_TONEAREST))
+#endif
+ full_val += (uint64_t(1) << ((41 - DROP_BITS) + d));
+
+ // RU
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (LIBC_UNLIKELY(rounding == FE_UPWARD)) {
+ constexpr uint64_t ROUND_UP_MASK = (uint64_t(1) << DROP_BITS) - 1;
+ full_val += ROUND_UP_MASK;
+ }
+#endif
// shift back to align to 32-bit float representation
uint32_t result = static_cast<uint32_t>(full_val >> ((42 - DROP_BITS) + d));
@@ -176,10 +212,12 @@ LIBC_INLINE float expf(float x) {
return cpp::bit_cast<float>(result);
}
-} // namespace integer_only
+} // namespace static_rounding
} // namespace math
+} // namespace shared
+
} // namespace LIBC_NAMESPACE_DECL
-#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
+#endif // LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_EXPF_H
diff --git a/libc/shared/math/static_rounding_math.h b/libc/shared/math/static_rounding_math.h
new file mode 100644
index 0000000000000..cb4b8a1544479
--- /dev/null
+++ b/libc/shared/math/static_rounding_math.h
@@ -0,0 +1,30 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the statically rounded implementations of math functions
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_MATH_H
+#define LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_MATH_H
+
+#include "shared/libc_common.h"
+#include "static_rounding/expf.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace shared {
+
+using math::static_rounding::expf;
+
+} // namespace shared
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_MATH_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 25ea9ba82e41e..9f4624682eafe 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3060,20 +3060,6 @@ add_header_library(
libc.src.__support.macros.config
)
-add_header_library(
- expf_integer_eval
- HDRS
- expf_integer_eval.h
- DEPENDS
- .exp_constants
- .exp_integer_utils
- libc.src.__support.CPP.bit
- libc.src.__support.FPUtil.fp_bits
- libc.src.__support.macros.config
- libc.src.__support.macros.optimization
- libc.src.__support.frac64
-)
-
add_header_library(
expf
HDRS
@@ -4113,19 +4099,6 @@ add_header_library(
libc.src.__support.macros.optimization
)
-add_header_library(
- exp_integer_utils
- HDRS
- exp_integer_utils.h
- DEPENDS
- .exp_constants
- libc.src.__support.CPP.bit
- libc.src.__support.FPUtil.fp_bits
- libc.src.__support.macros.config
- libc.src.__support.macros.optimization
- libc.src.__support.frac64
-)
-
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
deleted file mode 100644
index 0f9acf88cf28d..0000000000000
--- a/libc/src/__support/math/exp_integer_utils.h
+++ /dev/null
@@ -1,68 +0,0 @@
-//===----------------------------------------------------------------------===//
-//
-// 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
-//
-//===----------------------------------------------------------------------===//
-///
-/// \file
-/// This file contains utilities for implementing integer-only implementation
-/// of expf(x)
-///
-//===----------------------------------------------------------------------===//
-
-// TODO: this file is expected to be used for exp*(x) functions, not only
-// limited to expf(x)
-
-#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/frac64.h"
-#include "src/__support/macros/attributes.h"
-#include "src/__support/macros/config.h"
-#include "src/__support/macros/optimization.h"
-#include "src/__support/math_extras.h"
-
-namespace LIBC_NAMESPACE_DECL {
-
-namespace math {
-
-namespace integer_only {
-
-// 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);
-
-// 64-bit polynomial approximation of 2^x coefficients generated with Sollya:
-// > P = fpminimax(2^x, 11, [|1, 64...|], [0, 1], absolute, fixed);
-// Store the fractional part of the coefficients below
-// > dirtyinfnorm(2^x - P(x), [0, 1]);
-// 0x1.6238...p-58
-// LSB(EXPF_COEFFS[i]) = 2^-64
-LIBC_INLINE_VAR constexpr Frac64 EXPF_COEFFS[] = {
- Frac64(0xb172'17f7'd1cf'b7cf), // x
- Frac64(0x3d7f'7bff'057d'4a5e), // x^2
- Frac64(0x0e35'846b'8363'9484), // x^3
- Frac64(0x0276'556d'ec97'dcd4), // x^4
- Frac64(0x0057'61ff'dc04'c7ff), // x^5
- Frac64(0x000a'1847'b6e7'92ec), // x^6
- Frac64(0x0000'ffe8'14e5'7033), // x^7
- Frac64(0x0000'1628'b6e9'70c8), // x^8
- Frac64(0x0000'01b8'8ce7'4088), // x^9
- Frac64(0x0000'001c'18d5'cb29), // x^10
- Frac64(0x0000'0002'b43f'4490), // x^11
-};
-
-} // namespace integer_only
-
-} // namespace math
-
-} // namespace LIBC_NAMESPACE_DECL
-
-#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_UTILS_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 8fff101fa96a6..d5661f6f55e91 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1227,7 +1227,6 @@ 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 de8a97a46f1af..44c47d8ac42b6 100644
--- a/libc/src/math/generic/expf.cpp
+++ b/libc/src/math/generic/expf.cpp
@@ -7,8 +7,8 @@
//===----------------------------------------------------------------------===//
#include "src/math/expf.h"
+#include "shared/math/static_rounding/expf.h"
#include "src/__support/math/expf.h"
-#include "src/__support/math/expf_integer_eval.h"
namespace LIBC_NAMESPACE_DECL {
@@ -16,7 +16,7 @@ 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);
+ return shared::math::static_rounding::expf(x, FE_TONEAREST);
#else
return math::expf(x);
#endif
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 240689c0de02f..2f82df610bf80 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1193,6 +1193,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ expf_static_rounding_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ expf_static_rounding_test.cpp
+ DEPENDS
+ libc.src.errno.errno
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
expf16_test
NEED_MPFR
@@ -2569,7 +2581,7 @@ add_fp_unittest(
SRCS
asinf16_test.cpp
DEPENDS
- libc.src.math.asinf16
+ libc.src.math.asinf16
)
add_fp_unittest(
@@ -2650,7 +2662,7 @@ add_fp_unittest(
SRCS
acosf16_test.cpp
DEPENDS
- libc.src.math.acosf16
+ libc.src.math.acosf16
)
add_fp_unittest(
@@ -2672,7 +2684,7 @@ add_fp_unittest(
SRCS
acospif16_test.cpp
DEPENDS
- libc.src.math.acospif16
+ libc.src.math.acospif16
)
add_fp_unittest(
diff --git a/libc/test/src/math/exhaustive/CMakeLists.txt b/libc/test/src/math/exhaustive/CMakeLists.txt
index 2d1301d3a1e66..4a0feb89cb344 100644
--- a/libc/test/src/math/exhaustive/CMakeLists.txt
+++ b/libc/test/src/math/exhaustive/CMakeLists.txt
@@ -10,6 +10,16 @@ add_header_library(
libc.src.__support.macros.properties.types
)
+add_header_library(
+ exhaustive_test_static_rounding
+ HDRS
+ exhaustive_test_static_rounding.h
+ DEPENDS
+ libc.src.__support.CPP.type_traits
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.macros.properties.types
+)
+
add_fp_unittest(
sqrtf_test
NO_RUN_POSTBUILD
@@ -215,6 +225,21 @@ add_fp_unittest(
-lpthread
)
+add_fp_unittest(
+ expf_static_rounding_test
+ NO_RUN_POSTBUILD
+ SUITE
+ libc_math_exhaustive_tests
+ SRCS
+ expf_static_rounding_test.cpp
+ DEPENDS
+ .exhaustive_test_static_rounding
+ libc.src.__support.math.expf
+ libc.src.__support.FPUtil.fp_bits
+ LINK_LIBRARIES
+ -lpthread
+)
+
add_fp_unittest(
exp2f_test
NO_RUN_POSTBUILD
diff --git a/libc/test/src/math/exhaustive/exhaustive_test_static_rounding.h b/libc/test/src/math/exhaustive/exhaustive_test_static_rounding.h
new file mode 100644
index 0000000000000..5d248f5399018
--- /dev/null
+++ b/libc/test/src/math/exhaustive/exhaustive_test_static_rounding.h
@@ -0,0 +1,244 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains exhaustive test template for statically rounded math
+/// functions
+///
+//===----------------------------------------------------------------------===//
+
+// This file is modeled after exhaustive_test.h, modified for testing statically
+// rounded math functions.
+
+#include "exhaustive_test.h"
+#include "src/__support/CPP/type_traits.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/macros/properties/types.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/RoundingModeUtils.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#include <atomic>
+#include <functional>
+#include <iostream>
+#include <mutex>
+#include <sstream>
+#include <thread>
+#include <vector>
+
+template <typename OutType, typename InType = OutType>
+using StaticallyRoundedUnaryOp = OutType(InType, int);
+
+template <typename OutType, typename InType,
+ UnaryOp<OutType, InType> BaselineFunc,
+ StaticallyRoundedUnaryOp<OutType, InType> Func>
+struct StaticallyRoundedUnaryOpChecker
+ : public virtual LIBC_NAMESPACE::testing::Test {
+ using FloatType = InType;
+ using FPBits = LIBC_NAMESPACE::fputil::FPBits<FloatType>;
+ using StorageType = typename FPBits::StorageType;
+ using RoundingMode = LIBC_NAMESPACE::fputil::testing::RoundingMode;
+ using ForceRoundingMode = LIBC_NAMESPACE::fputil::testing::ForceRoundingMode;
+
+ // TODO: dedupe
+ LIBC_INLINE int rounding_mode_to_fenv_rounding_mode(RoundingMode rounding) {
+ switch (rounding) {
+ case RoundingMode::Nearest:
+ return FE_TONEAREST;
+ case RoundingMode::Downward:
+ return FE_DOWNWARD;
+ case RoundingMode::Upward:
+ return FE_UPWARD;
+ case RoundingMode::TowardZero:
+ return FE_TOWARDZERO;
+ }
+ return FE_TONEAREST; // Default case, should not happen
+ };
+
+ // Check in a range, return the number of failures.
+ uint64_t check(StorageType start, StorageType stop, RoundingMode rounding) {
+ ForceRoundingMode r(rounding);
+ if (!r.success)
+ return (stop > start);
+
+ const int fenv_rounding = rounding_mode_to_fenv_rounding_mode(rounding);
+
+ StorageType bits = start;
+ uint64_t failed = 0;
+
+ do {
+ FPBits xbits(bits);
+ FloatType x = xbits.get_val();
+ bool correct = TEST_FP_EQ(BaselineFunc(x), Func(x, fenv_rounding));
+ failed += (!correct);
+ // Uncomment to print out failed values.
+ if (!correct) {
+ EXPECT_FP_EQ_ROUNDING_MODE(BaselineFunc(x), Func(x, fenv_rounding),
+ rounding);
+ }
+ } while (bits++ < stop);
+
+ return failed;
+ }
+};
+
+// Modeled after `LlvmLibcExhaustiveMathTest`
+//
+// Checker class needs inherit from LIBC_NAMESPACE::testing::Test and provide
+// StorageType and check method.
+template <typename Checker, size_t Increment = 1 << 20>
+struct LlvmLibcExhaustiveStaticallyRoundedMathTest
+ : public virtual LIBC_NAMESPACE::testing::Test,
+ public Checker {
+ using FloatType = typename Checker::FloatType;
+ using FPBits = typename Checker::FPBits;
+ using StorageType = typename Checker::StorageType;
+ using RoundingMode = typename LIBC_NAMESPACE::fputil::testing::RoundingMode;
+
+ void explain_failed_range(std::stringstream &msg, StorageType x_begin,
+ StorageType x_end) {
+#ifdef LIBC_TYPES_HAS_FLOAT16
+ using T = LIBC_NAMESPACE::cpp::conditional_t<
+ LIBC_NAMESPACE::cpp::is_same_v<FloatType, float16>, float, FloatType>;
+#else
+ using T = FloatType;
+#endif
+
+ msg << x_begin << " to " << x_end << " [0x" << std::hex << x_begin << ", 0x"
+ << x_end << "), [" << std::hexfloat
+ << static_cast<T>(FPBits(x_begin).get_val()) << ", "
+ << static_cast<T>(FPBits(x_end).get_val()) << ")";
+ }
+
+ void explain_failed_range(std::stringstream &msg, StorageType x_begin,
+ StorageType x_end, StorageType y_begin,
+ StorageType y_end) {
+ msg << "x ";
+ explain_failed_range(msg, x_begin, x_end);
+ msg << ", y ";
+ explain_failed_range(msg, y_begin, y_end);
+ }
+
+ // Break [start, stop) into `nthreads` subintervals and apply *check to each
+ // subinterval in parallel.
+ template <typename... T>
+ void test_full_range(RoundingMode rounding, StorageType start,
+ StorageType stop, T... extra_range_bounds) {
+ int n_threads = std::thread::hardware_concurrency();
+ std::vector<std::thread> thread_list;
+ std::mutex mx_cur_val;
+ int current_percent = -1;
+ StorageType current_value = start;
+ std::atomic<uint64_t> failed(0);
+
+ for (int i = 0; i < n_threads; ++i) {
+ thread_list.emplace_back([&, this]() {
+ while (true) {
+ StorageType range_begin, range_end;
+ int new_percent = -1;
+ {
+ std::lock_guard<std::mutex> lock(mx_cur_val);
+ if (current_value == stop)
+ return;
+
+ range_begin = current_value;
+ if (stop >= Increment && stop - Increment >= current_value) {
+ range_end = static_cast<StorageType>(current_value + Increment);
+ } else {
+ range_end = stop;
+ }
+ current_value = range_end;
+ int pc =
+ static_cast<int>(100.0 * (range_end - start) / (stop - start));
+ if (current_percent != pc) {
+ new_percent = pc;
+ current_percent = pc;
+ }
+ }
+ if (new_percent >= 0) {
+ std::stringstream msg;
+ msg << new_percent << "% is in process \r";
+ std::cout << msg.str() << std::flush;
+ }
+
+ uint64_t failed_in_range = Checker::check(
+ range_begin, range_end, extra_range_bounds..., rounding);
+ if (failed_in_range > 0) {
+ std::stringstream msg;
+ msg << "Test failed for " << std::dec << failed_in_range
+ << " inputs in range: ";
+ explain_failed_range(msg, range_begin, range_end,
+ extra_range_bounds...);
+ msg << "\n";
+ std::cerr << msg.str() << std::flush;
+
+ failed.fetch_add(failed_in_range);
+ }
+ }
+ });
+ }
+
+ for (auto &thread : thread_list) {
+ if (thread.joinable()) {
+ thread.join();
+ }
+ }
+
+ std::cout << std::endl;
+ std::cout << "Test " << ((failed > 0) ? "FAILED" : "PASSED") << std::endl;
+ ASSERT_EQ(failed.load(), uint64_t(0));
+ }
+
+ void test_full_range_all_roundings(StorageType start, StorageType stop) {
+ std::cout << "-- Testing for FE_TONEAREST in range [0x" << std::hex << start
+ << ", 0x" << stop << ") --" << std::dec << std::endl;
+ test_full_range(RoundingMode::Nearest, start, stop);
+
+ std::cout << "-- Testing for FE_UPWARD in range [0x" << std::hex << start
+ << ", 0x" << stop << ") --" << std::dec << std::endl;
+ test_full_range(RoundingMode::Upward, start, stop);
+
+ std::cout << "-- Testing for FE_DOWNWARD in range [0x" << std::hex << start
+ << ", 0x" << stop << ") --" << std::dec << std::endl;
+ test_full_range(RoundingMode::Downward, start, stop);
+
+ std::cout << "-- Testing for FE_TOWARDZERO in range [0x" << std::hex
+ << start << ", 0x" << stop << ") --" << std::dec << std::endl;
+ test_full_range(RoundingMode::TowardZero, start, stop);
+ }
+
+ void test_full_range_all_roundings(StorageType x_start, StorageType x_stop,
+ StorageType y_start, StorageType y_stop) {
+ std::cout << "-- Testing for FE_TONEAREST in x range [0x" << std::hex
+ << x_start << ", 0x" << x_stop << "), y range [0x" << y_start
+ << ", 0x" << y_stop << ") --" << std::dec << std::endl;
+ test_full_range(RoundingMode::Nearest, x_start, x_stop, y_start, y_stop);
+
+ std::cout << "-- Testing for FE_UPWARD in x range [0x" << std::hex
+ << x_start << ", 0x" << x_stop << "), y range [0x" << y_start
+ << ", 0x" << y_stop << ") --" << std::dec << std::endl;
+ test_full_range(RoundingMode::Upward, x_start, x_stop, y_start, y_stop);
+
+ std::cout << "-- Testing for FE_DOWNWARD in x range [0x" << std::hex
+ << x_start << ", 0x" << x_stop << "), y range [0x" << y_start
+ << ", 0x" << y_stop << ") --" << std::dec << std::endl;
+ test_full_range(RoundingMode::Downward, x_start, x_stop, y_start, y_stop);
+
+ std::cout << "-- Testing for FE_TOWARDZERO in x range [0x" << std::hex
+ << x_start << ", 0x" << x_stop << "), y range [0x" << y_start
+ << ", 0x" << y_stop << ") --" << std::dec << std::endl;
+ test_full_range(RoundingMode::TowardZero, x_start, x_stop, y_start, y_stop);
+ }
+};
+
+template <typename FloatType, UnaryOp<FloatType> BaselineFunc,
+ StaticallyRoundedUnaryOp<FloatType> Func>
+using LlvmLibcStaticallyRoundedUnaryOpExhaustiveMathTest =
+ LlvmLibcExhaustiveStaticallyRoundedMathTest<StaticallyRoundedUnaryOpChecker<
+ FloatType, FloatType, BaselineFunc, Func>>;
diff --git a/libc/test/src/math/exhaustive/expf_static_rounding_test.cpp b/libc/test/src/math/exhaustive/expf_static_rounding_test.cpp
new file mode 100644
index 0000000000000..2c4dbeeebef78
--- /dev/null
+++ b/libc/test/src/math/exhaustive/expf_static_rounding_test.cpp
@@ -0,0 +1,32 @@
+//===-- Exhaustive test for expf ------------------------------------------===//
+//
+// 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 "exhaustive_test_static_rounding.h"
+#include "shared/math/static_rounding/expf.h"
+#include "src/__support/math/expf.h"
+
+using LlvmLibcStaticallyRoundedExpfExhaustiveTest =
+ LlvmLibcStaticallyRoundedUnaryOpExhaustiveMathTest<
+ float, LIBC_NAMESPACE::math::expf,
+ LIBC_NAMESPACE::shared::math::static_rounding::expf>;
+
+// Range: [0, Inf];
+static constexpr uint32_t POS_START = 0x0000'0000U;
+static constexpr uint32_t POS_STOP = 0x7f80'0000U;
+
+TEST_F(LlvmLibcStaticallyRoundedExpfExhaustiveTest, PostiveRange) {
+ test_full_range_all_roundings(POS_START, POS_STOP);
+}
+
+// Range: [-Inf, 0];
+static constexpr uint32_t NEG_START = 0xb000'0000U;
+static constexpr uint32_t NEG_STOP = 0xff80'0000U;
+
+TEST_F(LlvmLibcStaticallyRoundedExpfExhaustiveTest, NegativeRange) {
+ test_full_range_all_roundings(NEG_START, NEG_STOP);
+}
diff --git a/libc/test/src/math/expf_static_rounding_test.cpp b/libc/test/src/math/expf_static_rounding_test.cpp
new file mode 100644
index 0000000000000..806f779661912
--- /dev/null
+++ b/libc/test/src/math/expf_static_rounding_test.cpp
@@ -0,0 +1,148 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// This file contains the unit tests for statically-rounded implementation of
+/// static_rounding::expf(x)
+///
+//===----------------------------------------------------------------------===//
+
+#include "hdr/math_macros.h"
+#include "hdr/stdint_proxy.h"
+#include "shared/math/static_rounding/expf.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/libc_errno.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/expf.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/RoundingModeUtils.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcExpfStaticRoundingTest = LIBC_NAMESPACE::testing::FPTest<float>;
+using RoundingMode = LIBC_NAMESPACE::fputil::testing::RoundingMode;
+
+namespace static_rounding = LIBC_NAMESPACE::shared::math::static_rounding;
+namespace math = LIBC_NAMESPACE::math;
+
+LIBC_INLINE constexpr int
+rounding_mode_to_fenv_rounding_mode(RoundingMode rounding) {
+ switch (rounding) {
+ case RoundingMode::Nearest:
+ return FE_TONEAREST;
+ case RoundingMode::Downward:
+ return FE_DOWNWARD;
+ case RoundingMode::Upward:
+ return FE_UPWARD;
+ case RoundingMode::TowardZero:
+ return FE_TOWARDZERO;
+ }
+ return FE_TONEAREST; // Default case, should not happen
+};
+
+TEST_F(LlvmLibcExpfStaticRoundingTest, SpecialNumbers) {
+ for (auto rounding : ROUNDING_MODES) {
+ const int fenv_rounding = rounding_mode_to_fenv_rounding_mode(rounding);
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::expf(aNaN), static_rounding::expf(aNaN, fenv_rounding), rounding);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::expf(inf), static_rounding::expf(inf, fenv_rounding), rounding);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ROUNDING_MODE(math::expf(neg_inf),
+ static_rounding::expf(neg_inf, fenv_rounding),
+ rounding);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::expf(0.0f), static_rounding::expf(0.0f, fenv_rounding), rounding);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ROUNDING_MODE(math::expf(-0.0f),
+ static_rounding::expf(-0.0f, fenv_rounding),
+ rounding);
+ EXPECT_MATH_ERRNO(0);
+ }
+}
+
+TEST_F(LlvmLibcExpfStaticRoundingTest, Overflow) {
+ constexpr float VALUES[] = {FPBits(0x7f7fffffU).get_val(),
+ FPBits(0x42cffff8U).get_val(),
+ FPBits(0x42d00008U).get_val()};
+
+ for (auto rounding : ROUNDING_MODES) {
+ const int fenv_rounding = rounding_mode_to_fenv_rounding_mode(rounding);
+
+ // Statically rounded expf doesn't raise exceptions
+
+ for (auto x : VALUES) {
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::expf(x), static_rounding::expf(x, fenv_rounding), rounding);
+ }
+ }
+}
+
+TEST_F(LlvmLibcExpfStaticRoundingTest, Underflow) {
+ constexpr float VALUES[] = {FPBits(0xff7fffffU).get_val(),
+ FPBits(0xc2cffff8U).get_val(),
+ FPBits(0xc2d00008U).get_val()};
+
+ for (auto rounding : ROUNDING_MODES) {
+ const int fenv_rounding = rounding_mode_to_fenv_rounding_mode(rounding);
+
+ // Statically rounded expf doesn't raise exceptions
+ for (auto x : VALUES) {
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::expf(x), static_rounding::expf(x, fenv_rounding), rounding);
+ }
+ }
+}
+
+// Test with inputs which are the borders of underflow/overflow but still
+// produce valid results without setting errno.
+TEST_F(LlvmLibcExpfStaticRoundingTest, Borderline) {
+ constexpr float VALUES[] = {
+ FPBits(0x42affff8U).get_val(), FPBits(0x42b00008U).get_val(),
+ FPBits(0xc2affff8U).get_val(), FPBits(0xc2b00008U).get_val(),
+ FPBits(0xc236bd8cU).get_val()};
+
+ for (auto rounding : ROUNDING_MODES) {
+ const int fenv_rounding = rounding_mode_to_fenv_rounding_mode(rounding);
+
+ for (auto x : VALUES) {
+ // printf(
+ // "mode = %d, 0x%X, math::expf = 0x%X, static_rounding::expf =
+ // 0x%X\n", fenv_rounding, FPBits(x).uintval(),
+ // FPBits(math::expf(x)).uintval(), FPBits(static_rounding::expf(x,
+ // fenv_rounding)).uintval());
+ // TODO: fix failing case:
+ // 0xC2B00008, math::expf = 0x41ECBD, static_rounding::expf = 0x41ECBC
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::expf(x), static_rounding::expf(x, fenv_rounding), rounding);
+ }
+ }
+}
+
+TEST_F(LlvmLibcExpfStaticRoundingTest, InFloatRange) {
+ constexpr uint32_t COUNT = 1'231;
+ constexpr uint32_t STEP = UINT32_MAX / COUNT;
+ for (auto rounding : ROUNDING_MODES) {
+ const int fenv_rounding = rounding_mode_to_fenv_rounding_mode(rounding);
+
+ for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) {
+ float x = FPBits(v).get_val();
+ if (FPBits(v).is_nan() || FPBits(v).is_inf())
+ continue;
+ libc_errno = 0;
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::expf(x), static_rounding::expf(x, fenv_rounding), rounding);
+ }
+ }
+}
>From 9de7429ba6e3e80dd0975119d491a340311c092a Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Mon, 10 Aug 2026 20:04:23 +0700
Subject: [PATCH 17/17] fix: wrong calling condition, and cleaning up includes
---
libc/src/math/generic/expf.cpp | 6 +++---
.../src/math/exhaustive/exhaustive_test_static_rounding.h | 2 --
libc/test/src/math/expf_static_rounding_test.cpp | 1 -
3 files changed, 3 insertions(+), 6 deletions(-)
diff --git a/libc/src/math/generic/expf.cpp b/libc/src/math/generic/expf.cpp
index 44c47d8ac42b6..1ec621b22677f 100644
--- a/libc/src/math/generic/expf.cpp
+++ b/libc/src/math/generic/expf.cpp
@@ -13,9 +13,9 @@
namespace LIBC_NAMESPACE_DECL {
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)
+#if !defined(LIBC_TARGET_CPU_HAS_FPU_DOUBLE) && \
+ defined(LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY) && \
+ defined(LIBC_MATH_HAS_NO_EXCEPT) && defined(LIBC_MATH_HAS_NO_ERRNO)
return shared::math::static_rounding::expf(x, FE_TONEAREST);
#else
return math::expf(x);
diff --git a/libc/test/src/math/exhaustive/exhaustive_test_static_rounding.h b/libc/test/src/math/exhaustive/exhaustive_test_static_rounding.h
index 5d248f5399018..6ab34250fe341 100644
--- a/libc/test/src/math/exhaustive/exhaustive_test_static_rounding.h
+++ b/libc/test/src/math/exhaustive/exhaustive_test_static_rounding.h
@@ -22,10 +22,8 @@
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/RoundingModeUtils.h"
#include "test/UnitTest/Test.h"
-#include "utils/MPFRWrapper/MPFRUtils.h"
#include <atomic>
-#include <functional>
#include <iostream>
#include <mutex>
#include <sstream>
diff --git a/libc/test/src/math/expf_static_rounding_test.cpp b/libc/test/src/math/expf_static_rounding_test.cpp
index 806f779661912..66688c64f23e1 100644
--- a/libc/test/src/math/expf_static_rounding_test.cpp
+++ b/libc/test/src/math/expf_static_rounding_test.cpp
@@ -22,7 +22,6 @@
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/RoundingModeUtils.h"
#include "test/UnitTest/Test.h"
-#include "utils/MPFRWrapper/MPFRUtils.h"
using LlvmLibcExpfStaticRoundingTest = LIBC_NAMESPACE::testing::FPTest<float>;
using RoundingMode = LIBC_NAMESPACE::fputil::testing::RoundingMode;
More information about the libc-commits
mailing list