[libc-commits] [libc] [libc][math] Integer-only implementation of expf with 1-ULP errors (PR #209406)
via libc-commits
libc-commits at lists.llvm.org
Fri Jul 24 02:19:39 PDT 2026
=?utf-8?q?Hoàng_Minh_Thiên?=,=?utf-8?q?Hoàng_Minh_Thiên?=,
=?utf-8?q?Hoàng_Minh_Thiên?=,=?utf-8?q?Hoàng_Minh_Thiên?Message-ID:
In-Reply-To: <llvm.org/llvm/llvm-project/pull/209406 at github.com>
llvmorg-github-actions[bot] wrote:
<!--LLVM PR SUMMARY COMMENT-->
@llvm/pr-subscribers-libc
Author: Hoàng Minh Thiên (hmthien050209)
<details>
<summary>Changes</summary>
Integer-only implementation of expf with 1-ULP errors, using 1 single `Frac64`-based pipeline
# Accuracy
At most 1-ULP errors in the unit tests (with `LIBC_MATH_FAST`). All smoke tests passed.
# Code size:
> [!NOTE]
> The report below are from the`armv8m.main-none-eabi` triple.
## Before the patch
```sh
$ ls -lB libc/src/math/generic/CMakeFiles/libc.src.math.generic.expf.dir/
total 16
-rw-r-----@ 1 davidhoang staff 5520 Jul 24 15:59 expf.cpp.obj
$ llvm-nm --radix=d --print-size --size-sort --reverse-sort libc/src/math/generic/CMakeFiles/libc.src.math.generic.expf.dir/expf.cpp.obj
00000000 00001560 V _ZN22__llvm_libc_24_0_0_git4math6EXP_M1E
00000000 00001024 V _ZN22__llvm_libc_24_0_0_git4math6EXP_M2E
00000000 00000628 W _ZN22__llvm_libc_24_0_0_git4math4expfEf
00000000 00000012 T expf
00000000 00000012 T _ZN22__llvm_libc_24_0_0_git4expfEf
00000000 00000004 V _ZZN22__llvm_libc_24_0_0_git6fputil7generic16fenv_is_round_upEvE1x
00000000 00000004 V _ZZN22__llvm_libc_24_0_0_git6fputil7generic15quick_get_roundEvE1x
U __aeabi_f2d
U __aeabi_dmul
U __aeabi_dadd
U __aeabi_d2f
```
## After the patch
```sh
$ ls -lB libc/src/math/generic/CMakeFiles/libc.src.math.generic.expf.dir/
total 8
-rw-r-----@ 1 davidhoang staff 2724 Jul 24 15:54 expf.cpp.obj
$ llvm-nm --radix=d --print-size --size-sort --reverse-sort libc/src/math/generic/CMakeFiles/libc.src.math.generic.expf.dir/expf.cpp.obj
00000000 00001412 W _ZN22__llvm_libc_24_0_0_git4math12integer_only4expfEf
00000000 00000012 T expf
00000000 00000012 T _ZN22__llvm_libc_24_0_0_git4expfEf
```
---
Full diff: https://github.com/llvm/llvm-project/pull/209406.diff
7 Files Affected:
- (modified) libc/src/__support/CMakeLists.txt (+9)
- (added) libc/src/__support/frac64.h (+54)
- (modified) libc/src/__support/math/CMakeLists.txt (+27)
- (added) libc/src/__support/math/exp_integer_utils.h (+62)
- (added) libc/src/__support/math/expf_integer_eval.h (+174)
- (modified) libc/src/math/generic/CMakeLists.txt (+1)
- (modified) libc/src/math/generic/expf.cpp (+10-1)
``````````diff
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/frac64.h b/libc/src/__support/frac64.h
new file mode 100644
index 0000000000000..89296b67fa61a
--- /dev/null
+++ b/libc/src/__support/frac64.h
@@ -0,0 +1,54 @@
+//===-- 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 "src/__support/big_int.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+// Q0.64
+struct Frac64 : public UInt<64> {
+ using UInt<64>::UInt;
+
+ LIBC_INLINE constexpr Frac64 operator~() const { return Frac64(~val[0]); }
+
+ LIBC_INLINE constexpr Frac64 operator+(Frac64 other) const {
+ return Frac64(val[0] + other.val[0]);
+ }
+
+ LIBC_INLINE constexpr Frac64 operator-(Frac64 other) const {
+ return Frac64(val[0] - other.val[0]);
+ }
+
+ 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+=(Frac64 other) {
+ *this = *this + other;
+ return *this;
+ }
+
+ LIBC_INLINE constexpr Frac64 &operator-=(Frac64 other) {
+ *this = *this - other;
+ return *this;
+ }
+
+ LIBC_INLINE constexpr Frac64 &operator*=(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/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 9f4624682eafe..25ea9ba82e41e 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3060,6 +3060,20 @@ 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
@@ -4099,6 +4113,19 @@ 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
new file mode 100644
index 0000000000000..3fda41c6350c4
--- /dev/null
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -0,0 +1,62 @@
+//===-- 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
+//
+//===----------------------------------------------------------------------===//
+
+#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/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 {
+
+// 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
+
+} // 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..b4ad9ada54fa1
--- /dev/null
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -0,0 +1,174 @@
+//===-- 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
+
+#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 {
+
+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);
+}
+
+} // 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
``````````
</details>
https://github.com/llvm/llvm-project/pull/209406
More information about the libc-commits
mailing list