[libc-commits] [libc] [llvm] [libc][math] Optimize integer-only expf and factor out shared utilities. (PR #233118)
via libc-commits
libc-commits at lists.llvm.org
Sat Oct 17 10:36:24 PDT 2026
llvmorg-github-actions[bot] wrote:
<!--LLVM PR SUMMARY COMMENT-->
@llvm/pr-subscribers-libc
Author: lntue
<details>
<summary>Changes</summary>
- Optimize `expf_integer_eval` for both code size and performance.
- Factor out the shared range reduction and `2^u` fixed point evaluation to be used by `exp2f` and `exp10f` integer-only versions.
Assisted-by: Gemini for analysis and refactoring.
---
Patch is 37.79 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/233118.diff
14 Files Affected:
- (modified) libc/shared/math/static_rounding/expf.h (+2)
- (modified) libc/src/__support/CMakeLists.txt (+1)
- (modified) libc/src/__support/big_int.h (+48-10)
- (modified) libc/src/__support/frac64.h (+16-12)
- (modified) libc/src/__support/math/CMakeLists.txt (+19-2)
- (added) libc/src/__support/math/exp2f_integer_utils.h (+182)
- (modified) libc/src/__support/math/expf_integer_eval.h (+62-238)
- (modified) libc/test/shared/CMakeLists.txt (+1)
- (modified) libc/test/shared/shared_math_test.cpp (+36)
- (modified) libc/test/src/__support/big_int_test.cpp (+1)
- (modified) libc/test/src/math/exhaustive/expf_test.cpp (+1-1)
- (modified) libc/test/src/math/expf_test.cpp (+2-1)
- (modified) libc/test/src/math/smoke/expf_test.cpp (+2-1)
- (modified) utils/bazel/llvm-project-overlay/libc/BUILD.bazel (+18-1)
``````````diff
diff --git a/libc/shared/math/static_rounding/expf.h b/libc/shared/math/static_rounding/expf.h
index 0f4aeef673ba7a..78dd010f7be432 100644
--- a/libc/shared/math/static_rounding/expf.h
+++ b/libc/shared/math/static_rounding/expf.h
@@ -19,9 +19,11 @@
namespace LIBC_NAMESPACE_DECL {
namespace shared {
+namespace static_rounding {
using math::static_rounding::expf;
+} // namespace static_rounding
} // namespace shared
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/__support/CMakeLists.txt b/libc/src/__support/CMakeLists.txt
index cc31ced36cf2ee..dfb1ed0e13ff39 100644
--- a/libc/src/__support/CMakeLists.txt
+++ b/libc/src/__support/CMakeLists.txt
@@ -416,6 +416,7 @@ add_header_library(
frac64.h
DEPENDS
.big_int
+ libc.src.__support.macros.attributes
libc.src.__support.macros.config
)
diff --git a/libc/src/__support/big_int.h b/libc/src/__support/big_int.h
index 3bc9833664c13a..3461cd23728ccb 100644
--- a/libc/src/__support/big_int.h
+++ b/libc/src/__support/big_int.h
@@ -50,7 +50,7 @@ template <typename T> struct DoubleWide final : cpp::array<T, 2> {
};
// Converts an unsigned value into a DoubleWide<half_width_t<T>>.
-template <typename T> LIBC_INLINE constexpr auto split(T value) {
+template <typename T> LIBC_ALWAYS_INLINE constexpr auto split(T value) {
static_assert(cpp::is_unsigned_v<T>);
using half_type = half_width_t<T>;
return DoubleWide<half_type>(
@@ -59,19 +59,21 @@ template <typename T> LIBC_INLINE constexpr auto split(T value) {
}
// The low part of a DoubleWide value.
-template <typename T> LIBC_INLINE constexpr T lo(const DoubleWide<T> &value) {
+template <typename T>
+LIBC_ALWAYS_INLINE constexpr T lo(const DoubleWide<T> &value) {
return value[0];
}
// The high part of a DoubleWide value.
-template <typename T> LIBC_INLINE constexpr T hi(const DoubleWide<T> &value) {
+template <typename T>
+LIBC_ALWAYS_INLINE constexpr T hi(const DoubleWide<T> &value) {
return value[1];
}
// The low part of an unsigned value.
-template <typename T> LIBC_INLINE constexpr half_width_t<T> lo(T value) {
+template <typename T> LIBC_ALWAYS_INLINE constexpr half_width_t<T> lo(T value) {
return lo(split(value));
}
// The high part of an unsigned value.
-template <typename T> LIBC_INLINE constexpr half_width_t<T> hi(T value) {
+template <typename T> LIBC_ALWAYS_INLINE constexpr half_width_t<T> hi(T value) {
return hi(split(value));
}
@@ -239,9 +241,40 @@ LIBC_INLINE constexpr word multiply_with_carry(cpp::array<word, O> &dst,
}
template <typename word, size_t N>
-LIBC_INLINE constexpr void quick_mul_hi(cpp::array<word, N> &dst,
- const cpp::array<word, N> &lhs,
- const cpp::array<word, N> &rhs) {
+LIBC_ALWAYS_INLINE constexpr void quick_mul_hi(cpp::array<word, N> &dst,
+ const cpp::array<word, N> &lhs,
+ const cpp::array<word, N> &rhs) {
+ if constexpr (N == 1) {
+ if constexpr (cpp::is_same_v<word, uint8_t>) {
+ dst[0] = hi<uint16_t>(uint16_t(lhs[0]) * uint16_t(rhs[0]));
+ } else if constexpr (cpp::is_same_v<word, uint16_t>) {
+ dst[0] = hi<uint32_t>(uint32_t(lhs[0]) * uint32_t(rhs[0]));
+ }
+#ifdef LIBC_TYPES_HAS_INT64
+ else if constexpr (cpp::is_same_v<word, uint32_t>) {
+ dst[0] = hi<uint64_t>(uint64_t(lhs[0]) * uint64_t(rhs[0]));
+ }
+#endif // LIBC_TYPES_HAS_INT64
+#ifdef LIBC_TYPES_HAS_INT128
+ else if constexpr (cpp::is_same_v<word, uint64_t>) {
+ dst[0] = hi<__uint128_t>(__uint128_t(lhs[0]) * __uint128_t(rhs[0]));
+ }
+#endif // LIBC_TYPES_HAS_INT128
+ else {
+ // Split the single word into half-words and perform 3 half-word
+ // multiplications instead of 4.
+ word a = lhs[0];
+ word b = rhs[0];
+ word a_lo = lo(a);
+ word a_hi = hi(a);
+ word b_lo = lo(b);
+ word b_hi = hi(b);
+ word t1 = a_hi * b_lo;
+ word t2 = a_lo * b_hi + lo(t1);
+ dst[0] = a_hi * b_hi + hi(t1) + hi(t2);
+ }
+ return;
+ }
Accumulator<word> acc;
word carry = 0;
// First round of accumulation for those at N - 1 in the full product.
@@ -638,7 +671,11 @@ struct BigInt {
// `Bits` least significant bits of the full product, while this function will
// approximate `Bits` most significant bits of the full product with errors
// bounded by:
- // 0 <= (a.full_mul(b) >> Bits) - a.quick_mul_hi(b)) <= WORD_COUNT - 1.
+ // 0 <= (a.full_mul(b) >> Bits) - a.quick_mul_hi(b)) <= max(WORD_COUNT - 1,
+ // 1)
+ // (when WORD_COUNT == 1 and WordType is natively supported by a double-width
+ // integer type, the error is 0; otherwise when split into half-words, the
+ // error is bounded by 1).
//
// An example usage of this is to quickly (but less accurately) compute the
// product of (normalized) mantissas of floating point numbers:
@@ -651,11 +688,12 @@ struct BigInt {
// Performance summary:
// Number of 64-bit x 64-bit -> 128-bit multiplications performed.
// Bits WORD_COUNT ful_mul quick_mul_hi Error bound
+ // 64 1 1 1 0 (1 without int128)
// 128 2 4 3 1
// 196 3 9 6 2
// 256 4 16 10 3
// 512 8 64 36 7
- LIBC_INLINE constexpr BigInt quick_mul_hi(const BigInt &other) const {
+ LIBC_ALWAYS_INLINE constexpr BigInt quick_mul_hi(const BigInt &other) const {
BigInt result{};
multiword::quick_mul_hi(result.val, val, other.val);
return result;
diff --git a/libc/src/__support/frac64.h b/libc/src/__support/frac64.h
index 350b0a8b701b59..14fd5999f119b2 100644
--- a/libc/src/__support/frac64.h
+++ b/libc/src/__support/frac64.h
@@ -7,7 +7,7 @@
//===----------------------------------------------------------------------===//
///
/// \file
-/// This file contains the declaration of 64-bit unsigned fractional type
+/// This file contains the declaration of 64-bit unsigned fractional type.
///
//===----------------------------------------------------------------------===//
@@ -15,6 +15,7 @@
#define LLVM_LIBC_SRC___SUPPORT_FRAC64_H
#include "src/__support/big_int.h"
+#include "src/__support/macros/attributes.h"
#include "src/__support/macros/config.h"
namespace LIBC_NAMESPACE_DECL {
@@ -22,41 +23,44 @@ namespace LIBC_NAMESPACE_DECL {
struct Frac64 : public UInt<64> {
using UInt<64>::UInt;
- LIBC_INLINE constexpr Frac64 operator~() const { return Frac64(~val[0]); }
+ LIBC_ALWAYS_INLINE constexpr Frac64 operator~() const {
+ return Frac64(~val[0]);
+ }
- LIBC_INLINE constexpr Frac64 operator+(Frac64 other) const {
+ LIBC_ALWAYS_INLINE constexpr Frac64 operator+(Frac64 other) const {
return Frac64(val[0] + other.val[0]);
}
- LIBC_INLINE constexpr Frac64 operator-(Frac64 other) const {
+ LIBC_ALWAYS_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]);
+ // Fast hi part of the 128-bit product:
+ // 0 <= ((a * b) >> 64) - (Frac64(a) * Frac64(b)).val[0] <= 1.
+ LIBC_ALWAYS_INLINE constexpr Frac64 operator*(Frac64 other) const {
+ return Frac64(UInt<64>::quick_mul_hi(other).val);
}
- LIBC_INLINE constexpr Frac64 &operator+=(Frac64 other) {
+ LIBC_ALWAYS_INLINE constexpr Frac64 &operator+=(Frac64 other) {
*this = *this + other;
return *this;
}
- LIBC_INLINE constexpr Frac64 &operator-=(Frac64 other) {
+ LIBC_ALWAYS_INLINE constexpr Frac64 &operator-=(Frac64 other) {
*this = *this - other;
return *this;
}
- LIBC_INLINE constexpr Frac64 &operator*=(Frac64 other) {
+ LIBC_ALWAYS_INLINE constexpr Frac64 &operator*=(Frac64 other) {
*this = *this * other;
return *this;
}
- LIBC_INLINE constexpr Frac64 operator<<(size_t s) const {
+ LIBC_ALWAYS_INLINE constexpr Frac64 operator<<(size_t s) const {
return Frac64(val[0] << s);
}
- LIBC_INLINE constexpr Frac64 operator>>(size_t s) const {
+ LIBC_ALWAYS_INLINE constexpr Frac64 operator>>(size_t s) const {
return Frac64(val[0] >> s);
}
};
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index baac96e76d5435..c0a3430bac2b9d 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3110,6 +3110,21 @@ add_header_library(
libc.src.__support.macros.optimization
)
+add_header_library(
+ exp2f_integer_utils
+ HDRS
+ exp2f_integer_utils.h
+ DEPENDS
+ libc.hdr.fenv_macros
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.frac64
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
+)
+
add_header_library(
expf_double_eval
HDRS
@@ -3179,12 +3194,14 @@ add_header_library(
HDRS
expf_integer_eval.h
DEPENDS
- .exp_constants
+ .exp2f_integer_utils
+ libc.hdr.fenv_macros
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.frac64
+ libc.src.__support.macros.attributes
libc.src.__support.macros.config
libc.src.__support.macros.optimization
- libc.src.__support.frac64
)
add_header_library(
diff --git a/libc/src/__support/math/exp2f_integer_utils.h b/libc/src/__support/math/exp2f_integer_utils.h
new file mode 100644
index 00000000000000..569674436a0ed8
--- /dev/null
+++ b/libc/src/__support/math/exp2f_integer_utils.h
@@ -0,0 +1,182 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// Shared utilities and constants for integer-only single-precision exponential
+/// functions.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_INTEGER_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_INTEGER_UTILS_H
+
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/multiply_add.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/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace static_rounding {
+
+// Degree-11 minimax polynomial approximating 2^u - 1 on [0, 1], generated by
+// Sollya with:
+// > P = fpminimax(2^x, 11, [|1, 64...|], [0, 1], absolute, fixed);
+// > dirtyinfnorm(2^x - P(x), [0, 1]);
+// 0x1.6238...p-58
+// Store the fractional part of the coefficients below with LSB = 2^-64:
+LIBC_INLINE_VAR constexpr Frac64 EXP2F_COEFFS[11] = {
+ Frac64(0xb172'17f7'd1cf'b7cf), // u
+ Frac64(0x3d7f'7bff'057d'4a5e), // u^2
+ Frac64(0x0e35'846b'8363'9484), // u^3
+ Frac64(0x0276'556d'ec97'dcd4), // u^4
+ Frac64(0x0057'61ff'dc04'c7ff), // u^5
+ Frac64(0x000a'1847'b6e7'92ec), // u^6
+ Frac64(0x0000'ffe8'14e5'7033), // u^7
+ Frac64(0x0000'1628'b6e9'70c8), // u^8
+ Frac64(0x0000'01b8'8ce7'4088), // u^9
+ Frac64(0x0000'001c'18d5'cb29), // u^10
+ Frac64(0x0000'0002'b43f'4490), // u^11
+};
+
+// Compute |x| * log2(b) in unsigned 64-bit fixed-point format with binary
+// point at bit 54 (LSB = 2^-54), for x_e in [max_e - 31, max_e].
+// The binary point of log2_b is at bit (63 - (133 - max_e)).
+LIBC_ALWAYS_INLINE uint64_t exp_b_scale(uint32_t x_u, uint32_t x_e,
+ Frac64 log2_b, uint32_t max_e) {
+#if defined(LIBC_TYPES_HAS_INT128) && !defined(LIBC_MATH_HAS_SMALL_TABLES)
+ // When 128-bit integers are available (64-bit targets) and small code size is
+ // not requested, variable 64-bit shifts and 64x64 -> 128 multiplications are
+ // single native instructions. Shift the 24-bit significand of |x| into a
+ // 64-bit word and take the upper 64 bits of its 128-bit product with log2_b.
+ uint64_t m = static_cast<uint64_t>((x_u & 0x007f'ffffU) | 0x0080'0000U)
+ << (x_e - (max_e - 38U));
+ return (Frac64(m) * log2_b).val[0];
+#else // !defined(LIBC_TYPES_HAS_INT128) || defined(LIBC_MATH_HAS_SMALL_TABLES)
+ // On 32-bit targets without native 128-bit integers (or when optimizing for
+ // size), avoid 64-bit variable shifts (which lower to branches on 32-bit
+ // targets) and 64x64 -> 128 multiplications:
+ // 1. Place the 24-bit significand of |x| at bits [8..31] (leading 1 at
+ // bit 31) and compute the upper 64 bits of the 32x64 -> 96 product with
+ // (log2_b >> 2) using only two 32x32 -> 64 multiplications.
+ // 2. Since x_e is in [max_e - 31, max_e], s = max_e - x_e is in [0, 31], so
+ // shifting the 64-bit product right by s using two 32-bit shifts places
+ // the binary point of |x| * log2(b) at bit 54 without branches or UB at
+ // s = 0.
+ uint64_t m = (x_u << 8) | 0x8000'0000U;
+ uint64_t log2_b_scaled = log2_b.val[0] >> 2;
+ uint64_t prod = m * (log2_b_scaled >> 32) +
+ ((m * static_cast<uint32_t>(log2_b_scaled)) >> 32);
+ uint32_t s = max_e - x_e;
+ uint32_t prod_lo = static_cast<uint32_t>(prod);
+ uint32_t prod_hi = static_cast<uint32_t>(prod >> 32);
+ uint32_t u_lo = (prod_lo >> s) | ((prod_hi << 1) << (31U - s));
+ uint32_t u_hi = prod_hi >> s;
+ return (static_cast<uint64_t>(u_hi) << 32) | u_lo;
+#endif // defined(LIBC_TYPES_HAS_INT128) && !defined(LIBC_MATH_HAS_SMALL_TABLES)
+}
+
+// Evaluate 2^(x * log2(b)) from u_bits = |x| * log2(b) (with binary point at
+// bit 54) and return the IEEE 754 single-precision bit representation.
+LIBC_ALWAYS_INLINE uint32_t exp2f_eval_bits(uint64_t u_bits, bool is_neg,
+ [[maybe_unused]] int rounding) {
+ // Negate in two's complement when x < 0 so that floor(x * log2(b)) is the
+ // arithmetic right shift by 54, and the fractional part u in [0, 1) is the
+ // low 54 bits.
+ // Setting the LSB (| 1U) prevents InstCombine's demanded-bits analysis from
+ // narrowing the free 32-bit low-word extraction mask 0xffff'ffffU to
+ // 0xffff'fc00U inside quick_mul_hi on 32-bit targets. Even when the
+ // fractional part is 0, u = Frac64(1) satisfies (a * 1) >> 64 == 0, so p
+ // remains 0 in all rounding modes.
+ int64_t u_signed =
+ is_neg ? -static_cast<int64_t>(u_bits) : static_cast<int64_t>(u_bits);
+ int k =
+ (static_cast<int32_t>(static_cast<uint64_t>(u_signed) >> 32) >> 22) + 126;
+ Frac64 u((static_cast<uint64_t>(u_signed) << 10) | 1U);
+
+#if defined(LIBC_TYPES_HAS_INT128) && !defined(LIBC_MATH_HAS_SMALL_TABLES)
+ // When 128-bit integers are available (64-bit targets) and small code size is
+ // not requested, each Frac64 multiplication is a single native 64x64 -> 128
+ // multiply-high instruction, and out-of-order 64-bit CPUs have ample
+ // registers and multiple integer multiplier pipelines:
+ // 1. Evaluate p ~ 2^u - 1 on [0, 1) using an unrolled 2-way Estrin/Horner
+ // scheme (evaluating even and odd powers in parallel in u^2) to halve the
+ // critical-path latency of the polynomial evaluation.
+ // 2. Take a fast return path for the common case of normal outputs in
+ // round-to-nearest mode (k >= 0 && rounding == FE_TONEAREST), bypassing
+ // the subnormal and directed-rounding logic below.
+ Frac64 u2 = u * u;
+ Frac64 p_even =
+ fputil::polyeval(u2, EXP2F_COEFFS[0], EXP2F_COEFFS[2], EXP2F_COEFFS[4],
+ EXP2F_COEFFS[6], EXP2F_COEFFS[8], EXP2F_COEFFS[10]);
+ Frac64 p_odd =
+ fputil::polyeval(u2, EXP2F_COEFFS[1], EXP2F_COEFFS[3], EXP2F_COEFFS[5],
+ EXP2F_COEFFS[7], EXP2F_COEFFS[9]);
+ Frac64 p = fputil::multiply_add(p_odd, u, p_even) * u;
+
+ if (LIBC_LIKELY(k >= 0 && rounding == FE_TONEAREST))
+ return (((static_cast<uint32_t>(p.val[0] >> 40) + 1U) >> 1) +
+ (static_cast<uint32_t>(k + 1) << 23));
+#else // !defined(LIBC_TYPES_HAS_INT128) || defined(LIBC_MATH_HAS_SMALL_TABLES)
+ // On 32-bit targets without native 128-bit integers (such as 32-bit ARM
+ // microcontrollers without hardware double-precision FPU), or when
+ // LIBC_MATH_HAS_SMALL_TABLES is enabled, unrolling or using Estrin's scheme
+ // would increase register pressure and bloat code size (by ~470-810 bytes on
+ // 32-bit targets and ~70-360 bytes on 64-bit targets). Instead, evaluate
+ // p ~ 2^u - 1 using a compact, non-unrolled Horner loop over EXP2F_COEFFS and
+ // share the unified normal/subnormal rounding epilogue below.
+ Frac64 p(0);
+ LIBC_LOOP_NOUNROLL
+ for (int i = 10; i >= 0; --i)
+ p = (p + EXP2F_COEFFS[i]) * u;
+#endif // defined(LIBC_TYPES_HAS_INT128) && !defined(LIBC_MATH_HAS_SMALL_TABLES)
+
+ // Combine 2^u = 1 + p with 2^(k - 126):
+ // Place the explicit leading 1 of 2^u at bit 31 and the top 31 bits of p in
+ // bits [0..30].
+ // 1. Normal result (k >= 0, biased exponent k + 1 in [1, 254]):
+ // Shift right by 7 to place the leading 1 at bit 24 (which after >> 1
+ // contributes +1 to the biased exponent field, added to k << 23).
+ // 2. Subnormal result (k < 0, i.e., k in [-24, -1]):
+ // Shift right by 7 - k in [8, 31] to place both the leading 1 and the
+ // fractional bits directly at their subnormal bit positions.
+ uint32_t mant = static_cast<uint32_t>(p.val[0] >> 33) | 0x8000'0000U;
+ uint32_t shift = 7;
+ uint32_t exp_bits = static_cast<uint32_t>(k) << 23;
+
+ if (LIBC_UNLIKELY(k < 0)) {
+ shift = static_cast<uint32_t>(7 - k);
+ exp_bits = 0;
+ }
+
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (LIBC_UNLIKELY(rounding != FE_TONEAREST)) {
+ uint32_t round_up = 0;
+ if (LIBC_UNLIKELY(rounding == FE_UPWARD)) {
+ uint32_t mask = ~0U >> (31U - shift);
+ round_up = static_cast<uint32_t>(
+ ((mant & mask) | static_cast<uint32_t>(p.val[0])) != 0);
+ }
+
+ return ((mant >> shift) >> 1) + round_up + exp_bits;
+ }
+#endif // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+
+ // Round to nearest.
+ return (((mant >> shift) + 1U) >> 1) + exp_bits;
+}
+
+} // namespace static_rounding
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2F_INTEGER_UTILS_H
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 52736323d8d7ef..89b61870c4e98f 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -7,8 +7,7 @@
//===----------------------------------------------------------------------===//
///
/// \file
-/// This file contains the statically rounded, integer-only implementation of
-/// expf(x)
+/// Integer-only implementation of expf.
///
//===----------------------------------------------------------------------===//
@@ -18,278 +17,103 @@
#include "hdr/fenv_macros.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/attributes.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/exp2f_integer_utils.h"
namespace LIBC_NAMESPACE_DECL {
-
-namespace shared {
-
namespace math {
-
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 cons...
[truncated]
``````````
</details>
https://github.com/llvm/llvm-project/pull/233118
More information about the libc-commits
mailing list