[libc-commits] [libc] [libc][math] Integer-only, statically rounded implementation of exp (PR #225803)
Hoàng Minh Thiên via libc-commits
libc-commits at lists.llvm.org
Sun Sep 27 04:31:30 PDT 2026
https://github.com/hmthien050209 updated https://github.com/llvm/llvm-project/pull/225803
>From 651cd266e94b1f3a1d7a88ae8ce4d6eb167bab5d Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Wed, 23 Sep 2026 21:50:17 +0700
Subject: [PATCH 1/6] initial commit
---
libc/shared/math/static_rounding/exp.h | 28 ++
libc/shared/static_rounding_math.h | 1 +
libc/src/__support/math/CMakeLists.txt | 13 +
libc/src/__support/math/exp_integer_eval.h | 367 ++++++++++++++++++
libc/src/math/generic/exp.cpp | 13 +-
libc/test/src/math/CMakeLists.txt | 15 +
.../src/math/exp_static_rounding_test.cpp | 143 +++++++
libc/test/src/math/smoke/CMakeLists.txt | 13 +
.../math/smoke/exp_static_rounding_test.cpp | 45 +++
9 files changed, 637 insertions(+), 1 deletion(-)
create mode 100644 libc/shared/math/static_rounding/exp.h
create mode 100644 libc/src/__support/math/exp_integer_eval.h
create mode 100644 libc/test/src/math/exp_static_rounding_test.cpp
create mode 100644 libc/test/src/math/smoke/exp_static_rounding_test.cpp
diff --git a/libc/shared/math/static_rounding/exp.h b/libc/shared/math/static_rounding/exp.h
new file mode 100644
index 0000000000000..1b9ecd94e9532
--- /dev/null
+++ b/libc/shared/math/static_rounding/exp.h
@@ -0,0 +1,28 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 shared statically-rounded exp(x) function
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_EXP_H
+#define LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_EXP_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/exp_integer_eval.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::static_rounding::exp;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_EXP_H
diff --git a/libc/shared/static_rounding_math.h b/libc/shared/static_rounding_math.h
index 5488354a9b552..5f6bb411f7468 100644
--- a/libc/shared/static_rounding_math.h
+++ b/libc/shared/static_rounding_math.h
@@ -15,6 +15,7 @@
#define LLVM_LIBC_SHARED_STATIC_ROUNDING_MATH_H
#include "shared/libc_common.h"
+#include "shared/math/static_rounding/exp.h"
#include "shared/math/static_rounding/expf.h"
#endif // LLVM_LIBC_SHARED_STATIC_ROUNDING_MATH_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 8d6a1495e1c11..c534dcfed0fbc 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -4199,6 +4199,19 @@ add_header_library(
libc.src.__support.macros.optimization
)
+add_header_library(
+ exp_integer_eval
+ HDRS
+ exp_integer_eval.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.frac128
+)
+
add_header_library(
exp2
HDRS
diff --git a/libc/src/__support/math/exp_integer_eval.h b/libc/src/__support/math/exp_integer_eval.h
new file mode 100644
index 0000000000000..f19d725f68405
--- /dev/null
+++ b/libc/src/__support/math/exp_integer_eval.h
@@ -0,0 +1,367 @@
+//===----------------------------------------------------------------------===//
+//
+// 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, integer-only implementation of
+/// exp(x)
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
+
+#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/frac128.h"
+#include "src/__support/integer_literals.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/check/exp_exceptions.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace shared {
+
+namespace math {
+
+namespace static_rounding {
+
+using LIBC_NAMESPACE::operator""_u128;
+
+// print(2+round(1/log(2), 128, RN));
+// LSB(INV_LN2) = 2^-127
+LIBC_INLINE_VAR constexpr Frac128 INV_LN2_F128 =
+ Frac128(0xb8aa3b29'5c17f0bb'be87fed0'691d3e89_u128);
+
+// exp(x) for x from 0 to (0b1111/2^4) = 15/16
+// > for i from 0 to 15 do {
+// print(1+round(exp(1/(16-i)), 128, RN));
+// };
+// LSB(EXP_MID4[i]) = 2^-127
+LIBC_INLINE_VAR constexpr Frac128 EXP_MID4[] = {
+ Frac128(1u),
+ Frac128(0x08415abb'e9a76bea'd8d00cf1'12e4d4a9_u128),
+ Frac128(0x08d2ff22'4cc9cd75'daf50224'e6405613_u128),
+ Frac128(0x097a3089'bcde8358'06617ab2'2577bd36_u128),
+ Frac128(0x0a3c18b3'c8501827'a038960e'6554890c_u128),
+ Frac128(0x0b1fac01'4ab4cc6e'f83ad1c5'0d2fb016_u128),
+ Frac128(0x0c2e831f'eb8fa72a'da22cfa4'efa8e05c_u128),
+ Frac128(0x0d763d9a'd0069cd7'c11c604c'2b558c0e_u128),
+ Frac128(0x0f0add66'738b06a8'2a34e44f'7af67e71_u128),
+ Frac128(0x110b022d'b7ae67ce'76b441c2'7035c6a1_u128),
+ Frac128(0x13a8048b'71443b31'12f72650'b356e2df_u128),
+ Frac128(0x1736d169'0604545d'44b774b0'16630f77_u128),
+ Frac128(0x1c56ecf2'c5646740'b2bb19c5'bbfe54e1_u128),
+ Frac128(0x245af1e1'f40c333b'3de1db4d'd55f29a7_u128),
+ Frac128(0x32a36d8d'd1689885'3fb63793'fb4f42a0_u128),
+ Frac128(0x53094c70'f034de4b'96ff7d5b'6f99fcd9_u128),
+ Frac128(0xdbf0a8b1'45769535'5fb8ac40'4e7a79e4_u128),
+};
+
+// 128-bit polynomial approximation of e^x coefficients generated with Sollya:
+// > P = fpminimax(exp(x), 12, [|1, 128...|], [0, 1/16], absolute, fixed);
+// Store the fractional part of the coefficients below
+// > dirtyinfnorm(exp(x) - P(x), [0, 1/16]);
+// 0x1.9295...p-110
+// LSB(EXPF_COEFFS[i]) = 2^-128
+LIBC_INLINE_VAR constexpr Frac128 EXP_COEFFS[] = {
+ // degree-0 = 1, add back afterwards to reduce calc ops
+ Frac128(0xffffffff'ffffffff'ffffffff'acc77512_u128),
+ Frac128(0x80000000'00000000'0000015d'68ca804a_u128),
+ Frac128(0x2aaaaaaa'aaaaaaaa'aaa8a818'b0e4f3c0_u128),
+ Frac128(0x0aaaaaaa'aaaaaaaa'ac27ae0b'4b8e5377_u128),
+ Frac128(0x02222222'22222221'7c98b384'6025fb89_u128),
+ Frac128(0x005b05b0'5b05b088'd6cbde88'10bd52a0_u128),
+ Frac128(0x000d00d0'0d00c797'c6f4c09c'21513949_u128),
+ Frac128(0x0001a01a'01a12ab5'01d2062c'93e5cc17_u128),
+ Frac128(0x00002e3b'c7332843'170d9ff8'520123e8_u128),
+ Frac128(0x0000049f'954bbaf1'd7e43929'0c19789f_u128),
+ Frac128(0x0000006b'8c001fdf'5a004af5'0d6f47b4_u128),
+ Frac128(0x00000009'401bb484'e5aa7a08'b5a17885_u128),
+};
+
+// TODO: make 128-bit exp version work first, fast 64-bit and Ziv's test version
+// later on.
+// When rolling out the 64-bit fast path, take a look at expf also
+//
+// TODO: Ziv's test + pushing from Frac64 to full Frac128 pipeline on
+// hard-to-round cases Ziv test by extracting the last 12 bits of the result
+// (result & ((1u << 13) - 1)), cast into uint32_t, +4, and then | (1u << 12).
+// If the result is > bound, we have a hard to round case, then dial to Frac128
+// pipeline
+//
+// TODO: current implementation is mostly ported over from expf. Tricky input
+// tests are failing. Debug.
+//
+// TODO: test against CORE-MATH
+
+LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ bool is_neg = xbits.is_neg();
+ uint64_t x_val = xbits.uintval();
+ uint64_t x_val_abs = xbits.abs().uintval();
+
+ // TODO: dedupe? Ported over from the base LLVM-libc exp(x) implementation
+ // with same checks, just without exceptions/errnos.
+ //
+ // x < log(2^-1075) or x >= 0x1.6232bdd7abcd3p+9 or |x| < 2^-53.
+ if (LIBC_UNLIKELY(
+ x_val >= 0xc0874910d52d3052 ||
+ (x_val < 0xbca0000000000000 && x_val >= 0x40862e42fefa39f0) ||
+ x_val < 0x3ca0000000000000)) {
+ // |x| <= 2^-53
+ if (x_val_abs <= 0x3ca0'0000'0000'0000ULL) {
+ // exp(x) ~ 1 + x
+ return 1 + x;
+ }
+
+ // x <= log(2^-1075) || x >= 0x1.6232bdd7abcd3p+9 or inf/nan.
+
+ // x <= log(2^-1075) or -inf/nan
+ if (x_val >= 0xc087'4910'd52d'3052ULL) {
+ // exp(-Inf) = 0
+ if (xbits.is_inf())
+ return 0.0;
+
+ // exp(nan) = nan
+ if (xbits.is_nan())
+ return x;
+
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (rounding == FE_UPWARD)
+ return FPBits::min_subnormal().get_val();
+#endif
+ return 0.0;
+ }
+
+ // x >= round(log(MAX_NORMAL), D, RU) = 0x1.62e42fefa39fp+9 or +inf/nan
+ // x is finite
+ if (x_val < 0x7ff0'0000'0000'0000ULL) {
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
+ return FPBits::max_normal().get_val();
+#endif
+ }
+ // x is +inf or nan
+ return x + FPBits::inf().get_val();
+ }
+
+ // Main calculations
+
+ uint16_t x_e = xbits.get_biased_exponent();
+ uint64_t x_s = xbits.get_mantissa();
+
+ // The idea of the algorithm below is that:
+ // For x = hi + mid + low, with:
+ // - hi is an integer
+ // - hi = N * ln(2)
+ // - mid * 2^4 is an integer
+ // - lo is the remainder bits
+ // Then:
+ // exp(x) = 2^N * exp(mid) * exp(lo)
+ // With this formula:
+ // - Multiplying by 2^N is exact and cheap, via adding into the exponent
+ // field
+ // - exp(mid) can be calculated via the LUT declared above
+ // - exp(lo) ~ 1 + lo + a0 * lo^2 + ...
+ // Then we can construct exp(x) pretty easily, as hi, mid, lo bits can be
+ // separate and be used independently, then we only need to reconstruct in the
+ // final steps, which makes our life easier.
+ //
+ // Errors in computing in double-precision
+ // TODO
+
+ // Range reduction
+
+ // Recall binary64:
+ // 1 sign bit
+ // 11 exp bits
+ // 52 sig bits
+
+ // add leading bit = 1
+ x_s |= uint64_t(1) << FPBits::FRACTION_LEN;
+
+ // shift to top 64 bits --> decimal point at hidden bit that we've added
+ int x_e_unbiased = static_cast<int>(x_e) - FPBits::EXP_BIAS;
+
+ // Shift left by 128 - 53 = 75 bits for the fractional part of the double
+ // to align to the correct positions inside Frac128
+ // LSB(x_s_frac) = 2^-75
+ UInt128 x_s_shifted = static_cast<UInt128>(x_s) << 75;
+
+ // Apply the exponent: shift so the binary point is at the hidden bit
+ if (x_e_unbiased > 0) {
+ x_s_shifted <<= x_e_unbiased;
+ } else if (x_e_unbiased < 0) {
+ x_s_shifted >>= -x_e_unbiased;
+ }
+
+ // LSB(x_s_frac) = 2^-75
+ Frac128 x_s_frac(x_s_shifted);
+
+ // LSB(x_ln2) = 2^-74 (product of 2^-75 * 2^127 / 2^128 ~ 2^-75 effectively,
+ // but the integer part of x*log2(e) lands in the high word)
+ Frac128 x_ln2 = x_s_frac * INV_LN2_F128;
+
+ // We use the top 128-bit word of x_ln2 to extract N (the integer part of
+ // x * log2(e)). The fractional part drives the polynomial/LUT evaluation.
+ constexpr uint64_t FRAC_MASK = (uint64_t(1) << 52) - 1;
+ uint64_t x_ln2_hi = x_ln2.val[1];
+
+ uint64_t e_y, l2y_r_hi;
+ uint32_t e_y_unbiased;
+
+ // As Frac128 can't store the sign, we need to handle the sign separately:
+ // - Both branches are computing floor(x * log2(e)).
+ // - For negative x, we round up to the next multiple of 2^52, then clear the
+ // last 52 bits.
+ // - For positive x, we round down (just clear) the last 52 bits.
+ //
+ // Then, l2y_r_hi is the remainder of x * log2(e) after removing the integer
+ // part, which is used to look up EXP_MID4 and compute exp(lo) - 1.
+ //
+ // e_y_unbiased is biased exponent field, but already bit-positioned to the
+ // exponent field of the double representation.
+ if (LIBC_UNLIKELY(is_neg)) {
+ e_y = (x_ln2_hi + FRAC_MASK) & ~FRAC_MASK;
+ l2y_r_hi = e_y - x_ln2_hi;
+ e_y_unbiased = (FPBits::EXP_BIAS << 20) - static_cast<uint32_t>(e_y >> 43);
+ } else {
+ e_y = x_ln2_hi & ~FRAC_MASK;
+ l2y_r_hi = x_ln2_hi - e_y;
+ e_y_unbiased = (FPBits::EXP_BIAS << 20) + static_cast<uint32_t>(e_y >> 43);
+ }
+
+ uint32_t k = static_cast<uint32_t>(e_y >> 52);
+ int d = static_cast<int>(k) - FPBits::EXP_BIAS;
+
+ // d >= 53 --> k >= 1076
+ // --> guaranteed to be below 2^-1074
+ //
+ // underflow
+ if (LIBC_UNLIKELY(is_neg && d >= 53)) {
+ return 0.0;
+ }
+
+ // Extract the 4 mid bits (bits [51:48] of l2y_r_hi) for LUT index
+ uint16_t x_mid = static_cast<uint16_t>((l2y_r_hi >> 48) & 0xf);
+
+ // Extract the low 48 bits for the polynomial approximation of exp(lo).
+ // LSB(x_lo) = 2^-48
+ uint64_t x_lo = l2y_r_hi & ((uint64_t(1) << 48) - 1);
+
+ // Shift left to move all bits to the high part
+ // LSB(x_lo_frac) = 2^-128
+ Frac128 x_lo_frac(static_cast<UInt128>(x_lo) << 80);
+
+ Frac128 p =
+ x_lo_frac * fputil::polyeval(x_lo_frac, EXP_COEFFS[0], EXP_COEFFS[1],
+ EXP_COEFFS[2], EXP_COEFFS[3], EXP_COEFFS[4],
+ EXP_COEFFS[5], EXP_COEFFS[6], EXP_COEFFS[7],
+ EXP_COEFFS[8], EXP_COEFFS[9], EXP_COEFFS[10],
+ EXP_COEFFS[11]);
+
+ // With:
+ // - p = exp(lo) - 1 --> exp(lo) = p + 1
+ // - mid_val = exp(mid) = EXP_MID4[x_mid]
+ // We have:
+ // exp(mid) * exp(lo) = mid_val * (p + 1)
+ // (Workaround because we're dealing with fractional representation of things)
+ Frac128 mid_val = EXP_MID4[x_mid];
+ Frac128 result128 = mid_val * p + mid_val;
+
+ // We're computing with errors < worst-case errors, so tie-rounding never
+ // happens. Hence, round-to-nearest, tie-to-even is equivalent to
+ // round-to-nearest, tie-to-away. Which is what we're implementing below
+ // in the following order:
+ //
+ // 1. Shift so that the rounding bit is at bit-0
+ // 2. Add 1 for rounding
+ // 3. Perform another shift by 1
+ // 4. Depending on the rounding modes, adjust accordingly:
+ // a. Add 1 if rounding-up (0 if not)
+ // b. Add e_y_unbiased to the result (0 if the result is subnormal)
+
+ uint32_t shift_length = 72;
+ uint32_t leading_one = 0;
+
+ // subnormal
+ if (LIBC_UNLIKELY(is_neg && d >= 0)) {
+ e_y_unbiased = 0;
+ leading_one = 1 << (52 - d);
+
+ // In the below shifts, we're shifting by (shift_length + 1) at max, while
+ // shift_length is already 72, and if d = 52 --> shift_length + d = 124,
+ // and we'll shift by whole 128 bits, which is undefined behavior in C++.
+ //
+ // So, we'll truncate the last 2 bits.
+ if (d >= 51) {
+ d -= 2;
+ result128.val[1] >>= 2;
+ }
+
+ shift_length += d + 1;
+ }
+
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ uint64_t result = (static_cast<uint64_t>(result128.val[1] >> shift_length) +
+ (static_cast<uint64_t>(leading_one) + 1));
+ result >>= 1;
+ result += static_cast<uint64_t>(e_y_unbiased) << 32;
+
+ return cpp::bit_cast<double>(result);
+#else // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (rounding == FE_TONEAREST) {
+ uint64_t result = (static_cast<uint64_t>(result128.val[1] >> shift_length) +
+ (static_cast<uint64_t>(leading_one) + 1));
+ result >>= 1;
+ result += static_cast<uint64_t>(e_y_unbiased) << 32;
+
+ return cpp::bit_cast<double>(result);
+ }
+
+ uint64_t should_round_up = 0;
+
+ if (LIBC_UNLIKELY(rounding == FE_UPWARD)) {
+ uint64_t round_up_mask = (uint64_t(1) << (shift_length + 1)) - 1;
+ should_round_up =
+ static_cast<uint64_t>((result128.val[1] & round_up_mask) != 0);
+ }
+
+ uint64_t result =
+ (static_cast<uint64_t>(result128.val[1] >> (shift_length + 1)) +
+ should_round_up + (static_cast<uint64_t>(leading_one) >> 1));
+ result += static_cast<uint64_t>(e_y_unbiased) << 32;
+
+ return cpp::bit_cast<double>(result);
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+}
+
+} // namespace static_rounding
+
+} // namespace math
+
+} // namespace shared
+
+namespace math {
+namespace integer_eval {
+
+LIBC_INLINE double exp(double x) {
+ return shared::math::static_rounding::exp(x, FE_TONEAREST);
+}
+
+} // namespace integer_eval
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
diff --git a/libc/src/math/generic/exp.cpp b/libc/src/math/generic/exp.cpp
index dc4d2ca480cb8..45f105aa25a90 100644
--- a/libc/src/math/generic/exp.cpp
+++ b/libc/src/math/generic/exp.cpp
@@ -7,9 +7,20 @@
//===----------------------------------------------------------------------===//
#include "src/math/exp.h"
+#include "shared/math/static_rounding/exp.h"
#include "src/__support/math/exp.h"
namespace LIBC_NAMESPACE_DECL {
-LLVM_LIBC_FUNCTION(double, exp, (double x)) { return math::exp(x); }
+LLVM_LIBC_FUNCTION(double, exp, (double x)) {
+// TODO: update to follow the new structure (following
+// https://github.com/llvm/llvm-project/pull/224735)
+#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::exp(x, FE_TONEAREST);
+#else
+ return math::exp(x);
+#endif
+}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 1fd0772aa5031..0c3b7a2fca73f 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1367,6 +1367,21 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ exp_static_rounding_test
+ SUITE
+ libc-math-unittests
+ SRCS
+ exp_static_rounding_test.cpp
+ DEPENDS
+ libc.hdr.math_macros
+ libc.hdr.stdint_proxy
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.libc_errno
+ libc.src.__support.macros.optimization
+ libc.src.__support.math.exp
+ libc.src.__support.math.exp_integer_eval
+)
add_fp_unittest(
expf16_test
diff --git a/libc/test/src/math/exp_static_rounding_test.cpp b/libc/test/src/math/exp_static_rounding_test.cpp
new file mode 100644
index 0000000000000..f3fc4a5796ee2
--- /dev/null
+++ b/libc/test/src/math/exp_static_rounding_test.cpp
@@ -0,0 +1,143 @@
+//===----------------------------------------------------------------------===//
+//
+// 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::exp(x)
+///
+//===----------------------------------------------------------------------===//
+
+#include "hdr/math_macros.h"
+#include "hdr/stdint_proxy.h"
+#include "shared/static_rounding_math.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/libc_errno.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/exp.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/RoundingModeUtils.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcExpStaticRoundingTest = LIBC_NAMESPACE::testing::FPTest<double>;
+using RoundingMode = LIBC_NAMESPACE::fputil::testing::RoundingMode;
+using ForceRoundingMode = LIBC_NAMESPACE::fputil::testing::ForceRoundingMode;
+using LIBC_NAMESPACE::testing::tlog;
+
+namespace static_rounding = LIBC_NAMESPACE::shared::math::static_rounding;
+namespace math = LIBC_NAMESPACE::math;
+
+TEST_F(LlvmLibcExpStaticRoundingTest, SpecialNumbers) {
+ constexpr double VALUES[] = {aNaN, inf, neg_inf, -0x1.0p20,
+ 0x1.0p20, zero, neg_zero};
+
+ for (auto rounding : ROUNDING_MODES) {
+ const int fenv_rounding = get_fe_rounding(rounding);
+
+ // Statically rounded exp doesn't raise exceptions
+ for (auto x : VALUES) {
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::exp(x), static_rounding::exp(x, fenv_rounding), rounding);
+ }
+ }
+}
+
+TEST_F(LlvmLibcExpStaticRoundingTest, TrickyInputs) {
+ constexpr uint64_t VALUES[] = {
+ 0x3FD79289C6E6A5C0,
+ 0x3FD05DE80A173EA0, // 0x1.05de80a173eap-2
+ 0xbf1eb7a4cb841fcc, // -0x1.eb7a4cb841fccp-14
+ 0xbf19a61fb925970d,
+ 0x3fda7b764e2cf47a, // 0x1.a7b764e2cf47ap-2
+ 0xc04757852a4b93aa, // -0x1.757852a4b93aap+5
+ 0x4044c19e5712e377, // x=0x1.4c19e5712e377p+5
+ 0xbf19a61fb925970d, // x=-0x1.9a61fb925970dp-14
+ 0xc039a74cdab36c28, // x=-0x1.9a74cdab36c28p+4
+ 0xc085b3e4e2e3bba9, // x=-0x1.5b3e4e2e3bba9p+9
+ 0xc086960d591aec34, // x=-0x1.6960d591aec34p+9
+ 0xc086232c09d58d91, // x=-0x1.6232c09d58d91p+9
+ 0xc0874910d52d3051, // x=-0x1.74910d52d3051p9
+ 0xc0867a172ceb0990, // x=-0x1.67a172ceb099p+9
+ };
+
+ for (auto rounding : ROUNDING_MODES) {
+ const int fenv_rounding = get_fe_rounding(rounding);
+
+ // Statically rounded exp doesn't raise exceptions
+ for (auto val : VALUES) {
+ double x = FPBits(val).get_val();
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::exp(x), static_rounding::exp(x, fenv_rounding), rounding);
+ }
+ }
+}
+
+// TODO: double has a very big range. Define test functions that allows for
+// finding max ULPs over a range. EXPECT_* just isn't enough for this.
+TEST_F(LlvmLibcExpStaticRoundingTest, InDoubleRange) {
+ constexpr uint64_t COUNT = 1'231;
+ constexpr uint64_t START = FPBits(0.25).uintval();
+ constexpr uint64_t STOP = FPBits(4.0).uintval();
+ constexpr uint64_t STEP = (STOP - START) / COUNT;
+
+ auto test = [&](RoundingMode rounding) {
+ const int fenv_rounding = get_fe_rounding(rounding);
+
+ uint64_t fails = 0;
+ uint64_t count = 0;
+ uint64_t cc = 0;
+ double mx, mr = 0.0;
+ double tol = 0.5;
+
+ for (uint64_t i = 0, v = START; i <= COUNT; ++i, v += STEP) {
+ double x = FPBits(v).get_val();
+ if (FPBits(v).is_nan() || FPBits(v).is_inf() || x < 0.0)
+ continue;
+ double result = math::exp(x);
+ ++cc;
+ if (FPBits(result).is_nan() || FPBits(result).is_inf())
+ continue;
+
+ ++count;
+ // ASSERT_MPFR_MATCH(mpfr::Operation::Log, x, result, 0.5);
+ // if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Exp, x, result,
+ // TOLERANCE + 0.5, rounding_mode))
+ // {
+ // ++fails;
+ // while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Exp, x,
+ // result, tol,
+ // rounding_mode)) {
+ // mx = x;
+ // mr = result;
+
+ // if (tol > 1000.0)
+ // break;
+
+ // tol *= 2.0;
+ // }
+ // }
+ }
+ if (fails) {
+ tlog << " Statically rounded exp failed: " << fails << "/" << count << "/"
+ << cc << " tests.\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
+ // EXPECT_MPFR_MATCH(mpfr::Operation::Exp, mx, mr, 0.5, rounding_mode);
+ }
+ };
+
+ tlog << " Test Rounding To Nearest...\n";
+ test(RoundingMode::Nearest);
+
+ tlog << " Test Rounding Downward...\n";
+ test(RoundingMode::Downward);
+
+ tlog << " Test Rounding Upward...\n";
+ test(RoundingMode::Upward);
+
+ tlog << " Test Rounding Toward Zero...\n";
+ test(RoundingMode::TowardZero);
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 0bfe73a2e643a..588aead6d78d7 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -1428,6 +1428,19 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ exp_static_rounding_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ exp_static_rounding_test.cpp
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.src.__support.math.exp
+ libc.src.__support.math.exp_integer_eval
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
expf_test
SUITE
diff --git a/libc/test/src/math/smoke/exp_static_rounding_test.cpp b/libc/test/src/math/smoke/exp_static_rounding_test.cpp
new file mode 100644
index 0000000000000..0fd7c4066ed43
--- /dev/null
+++ b/libc/test/src/math/smoke/exp_static_rounding_test.cpp
@@ -0,0 +1,45 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 smoke tests for static_rounding::exp(x)
+///
+//===----------------------------------------------------------------------===//
+
+#include "hdr/errno_macros.h"
+#include "hdr/math_macros.h"
+#include "hdr/stdint_proxy.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/math/exp.h"
+#include "src/__support/math/exp_integer_eval.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcExpStaticRoundingTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+namespace static_rounding = LIBC_NAMESPACE::shared::math::static_rounding;
+namespace math = LIBC_NAMESPACE::math;
+
+TEST_F(LlvmLibcExpStaticRoundingTest, SpecialNumbers) {
+ using LIBC_NAMESPACE::fputil::testing::get_fe_rounding;
+
+ constexpr double VALUES[] = {sNaN, aNaN, inf, neg_inf,
+ -0x1.0p20, 0x1.0p20, zero, neg_zero};
+
+ for (auto rounding : ROUNDING_MODES) {
+ const int fenv_rounding = get_fe_rounding(rounding);
+
+ for (auto x : VALUES) {
+ EXPECT_FP_EQ_ROUNDING_MODE(
+ math::exp(x), static_rounding::exp(x, fenv_rounding), rounding);
+ // Statically rounded exp doesn't raise exceptions, but the baseline
+ // exp may raise overflow exception.
+ // So, we won't check for that here.
+ }
+ }
+}
>From 4498d179b5a0d1075dc20a869ebd558de30b0e84 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Sat, 26 Sep 2026 22:26:31 +0700
Subject: [PATCH 2/6] feat: hybrid pipeline
still having errors, debugging in progress
---
libc/src/__support/math/exp_integer_eval.h | 428 ++++++++++++---------
libc/src/math/generic/exp.cpp | 20 +-
2 files changed, 260 insertions(+), 188 deletions(-)
diff --git a/libc/src/__support/math/exp_integer_eval.h b/libc/src/__support/math/exp_integer_eval.h
index f19d725f68405..ae06e9ed44fec 100644
--- a/libc/src/__support/math/exp_integer_eval.h
+++ b/libc/src/__support/math/exp_integer_eval.h
@@ -17,6 +17,8 @@
#include "hdr/fenv_macros.h"
#include "src/__support/CPP/bit.h"
+#include "src/__support/CPP/type_traits/enable_if.h"
+#include "src/__support/CPP/type_traits/is_same.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/frac128.h"
@@ -25,6 +27,63 @@
#include "src/__support/macros/optimization.h"
#include "src/__support/math/check/exp_exceptions.h"
+// #include <iomanip>
+// #include <iostream>
+// #include <string>
+// #include <type_traits>
+// #include <vector>
+//
+// std::vector<std::string> split_top_level(const std::string &s) {
+// std::vector<std::string> parts;
+// int depth = 0;
+// size_t start = 0;
+// for (size_t i = 0; i < s.size(); ++i) {
+// char c = s[i];
+// if (c == '(' || c == '[' || c == '{')
+// ++depth;
+// else if (c == ')' || c == ']' || c == '}')
+// --depth;
+// else if (c == ',' && depth == 0) {
+// parts.push_back(s.substr(start, i - start));
+// start = i + 1;
+// }
+// }
+// parts.push_back(s.substr(start));
+// for (auto &p : parts) {
+// size_t a = p.find_first_not_of(" \t\n");
+// size_t b = p.find_last_not_of(" \t\n");
+// p = (a == std::string::npos) ? "" : p.substr(a, b - a + 1);
+// }
+// return parts;
+// }
+//
+// template <typename T>
+// void debug_value(const std::string &name, const T &value) {
+// using U = std::remove_cv_t<std::remove_reference_t<T>>;
+// std::cout << std::left << std::setw(24) << name << " = ";
+// if constexpr (std::is_floating_point_v<U>) {
+// std::cout << std::defaultfloat << std::setw(24) << value << " ("
+// << std::hexfloat << std::right << std::setw(24) << value
+// << std::defaultfloat << std::left << ")";
+// } else if constexpr (std::is_integral_v<U>) {
+// std::cout << std::dec << std::setw(24) << value << " (" << std::right
+// << std::setw(24) << std::showbase << std::hex << value
+// << std::noshowbase << std::dec << ")";
+// } else {
+// std::cout << value;
+// }
+// std::cout << '\n';
+// }
+//
+// template <typename... Ts> void debug_all(const char *names, Ts &&...values) {
+// auto parts = split_top_level(names);
+// std::cout << "===== DEBUG =====\n";
+// size_t i = 0;
+// (debug_value(parts[i++], values), ...);
+// }
+//
+// #define DEBUG(...) debug_all(#__VA_ARGS__, __VA_ARGS__)
+
namespace LIBC_NAMESPACE_DECL {
namespace shared {
@@ -33,74 +92,130 @@ namespace math {
namespace static_rounding {
-using LIBC_NAMESPACE::operator""_u128;
+// TODO: tricky input tests are failing. Debug. Getting the last bits
+// truncating to all-0 in all rounding modes.
+// TODO: test against CORE-MATH
+// TODO: refactor to follow the new structure; dedupe codes
// print(2+round(1/log(2), 128, RN));
// LSB(INV_LN2) = 2^-127
LIBC_INLINE_VAR constexpr Frac128 INV_LN2_F128 =
- Frac128(0xb8aa3b29'5c17f0bb'be87fed0'691d3e89_u128);
+ Frac128({0xbe87'fed0'691d'3e89ULL, 0xb8aa'3b29'5c17'f0bbULL});
-// exp(x) for x from 0 to (0b1111/2^4) = 15/16
+// 2^x for x from 0 to (0b1111/2^4) = 15/16
// > for i from 0 to 15 do {
-// print(1+round(exp(1/(16-i)), 128, RN));
+// print(1+round(2^(i/16), 128, RN));
// };
// LSB(EXP_MID4[i]) = 2^-127
-LIBC_INLINE_VAR constexpr Frac128 EXP_MID4[] = {
- Frac128(1u),
- Frac128(0x08415abb'e9a76bea'd8d00cf1'12e4d4a9_u128),
- Frac128(0x08d2ff22'4cc9cd75'daf50224'e6405613_u128),
- Frac128(0x097a3089'bcde8358'06617ab2'2577bd36_u128),
- Frac128(0x0a3c18b3'c8501827'a038960e'6554890c_u128),
- Frac128(0x0b1fac01'4ab4cc6e'f83ad1c5'0d2fb016_u128),
- Frac128(0x0c2e831f'eb8fa72a'da22cfa4'efa8e05c_u128),
- Frac128(0x0d763d9a'd0069cd7'c11c604c'2b558c0e_u128),
- Frac128(0x0f0add66'738b06a8'2a34e44f'7af67e71_u128),
- Frac128(0x110b022d'b7ae67ce'76b441c2'7035c6a1_u128),
- Frac128(0x13a8048b'71443b31'12f72650'b356e2df_u128),
- Frac128(0x1736d169'0604545d'44b774b0'16630f77_u128),
- Frac128(0x1c56ecf2'c5646740'b2bb19c5'bbfe54e1_u128),
- Frac128(0x245af1e1'f40c333b'3de1db4d'd55f29a7_u128),
- Frac128(0x32a36d8d'd1689885'3fb63793'fb4f42a0_u128),
- Frac128(0x53094c70'f034de4b'96ff7d5b'6f99fcd9_u128),
- Frac128(0xdbf0a8b1'45769535'5fb8ac40'4e7a79e4_u128),
+LIBC_INLINE_VAR constexpr Frac128 EXP_MID[] = {
+ Frac128({0x0000'0000'0000'0000ULL, 0x8000'0000'0000'0000ULL}),
+ Frac128({0xc5c9'5b8c'2154'c1b2ULL, 0x85aa'c367'cc48'7b14ULL}),
+ Frac128({0xfbe4'6287'58a5'3c90ULL, 0x8b95'c1e3'ea8b'd6e6ULL}),
+ Frac128({0x0fd6'd8e0'ae5a'c9d8ULL, 0x91c3'd373'ab11'c336ULL}),
+ Frac128({0x46ad'2318'2e42'f6f6ULL, 0x9837'f051'8db8'a96fULL}),
+ Frac128({0xa091'1f09'ebb9'fdd1ULL, 0x9ef5'3260'91a1'11adULL}),
+ Frac128({0x1cbd'7f62'1710'701bULL, 0xa5fe'd6a9'b151'38eaULL}),
+ Frac128({0x4980'a8c8'f59a'2ec4ULL, 0xad58'3eea'42a1'4ac6ULL}),
+ Frac128({0x597d'89b3'754a'be9fULL, 0xb504'f333'f9de'6484ULL}),
+ Frac128({0xa881'1fb6'6d0f'af7aULL, 0xbd08'a39f'580c'36beULL}),
+ Frac128({0x3e2a'd0c9'64dd'9f37ULL, 0xc567'2a11'5506'daddULL}),
+ Frac128({0xe235'838f'95f2'c6edULL, 0xce24'8c15'1f84'80e3ULL}),
+ Frac128({0x39a6'8bb9'902d'3fdeULL, 0xd744'fcca'd69d'6af4ULL}),
+ Frac128({0x0658'9504'8dd3'33caULL, 0xe0cc'deec'2a94'e111ULL}),
+ Frac128({0xd02d'75b3'706e'54fbULL, 0xeac0'c6e7'dd24'392eULL}),
+ Frac128({0x7b9d'0c7a'ed98'0fc3ULL, 0xf525'7d15'2486'cc2cULL}),
};
-// 128-bit polynomial approximation of e^x coefficients generated with Sollya:
-// > P = fpminimax(exp(x), 12, [|1, 128...|], [0, 1/16], absolute, fixed);
+// 128-bit polynomial approximation of 2^x coefficients generated with Sollya:
+// > P = fpminimax(2^x, 12, [|1, 128...|], [0, 1/16], absolute, fixed);
// Store the fractional part of the coefficients below
-// > dirtyinfnorm(exp(x) - P(x), [0, 1/16]);
-// 0x1.9295...p-110
+// > dirtyinfnorm(2^x - P(x), [0, 1/16]);
+// 0x1.b328...p-117
// LSB(EXPF_COEFFS[i]) = 2^-128
LIBC_INLINE_VAR constexpr Frac128 EXP_COEFFS[] = {
// degree-0 = 1, add back afterwards to reduce calc ops
- Frac128(0xffffffff'ffffffff'ffffffff'acc77512_u128),
- Frac128(0x80000000'00000000'0000015d'68ca804a_u128),
- Frac128(0x2aaaaaaa'aaaaaaaa'aaa8a818'b0e4f3c0_u128),
- Frac128(0x0aaaaaaa'aaaaaaaa'ac27ae0b'4b8e5377_u128),
- Frac128(0x02222222'22222221'7c98b384'6025fb89_u128),
- Frac128(0x005b05b0'5b05b088'd6cbde88'10bd52a0_u128),
- Frac128(0x000d00d0'0d00c797'c6f4c09c'21513949_u128),
- Frac128(0x0001a01a'01a12ab5'01d2062c'93e5cc17_u128),
- Frac128(0x00002e3b'c7332843'170d9ff8'520123e8_u128),
- Frac128(0x0000049f'954bbaf1'd7e43929'0c19789f_u128),
- Frac128(0x0000006b'8c001fdf'5a004af5'0d6f47b4_u128),
- Frac128(0x00000009'401bb484'e5aa7a08'b5a17885_u128),
+ Frac128({0xc9e3'b398'033f'0902ULL, 0xb172'17f7'd1cf'79abULL}),
+ Frac128({0xde2d'60e0'866c'6365ULL, 0x3d7f'7bff'058b'1d50ULL}),
+ Frac128({0x99d3'ad04'd47d'0efeULL, 0x0e35'846b'8250'5fc5ULL}),
+ Frac128({0x399a'b423'1644'1b55ULL, 0x0276'556d'f749'cee5ULL}),
+ Frac128({0x405f'ebef'0f76'011aULL, 0x0057'61ff'9e29'9cc4ULL}),
+ Frac128({0x1ac4'9999'c0ef'b9abULL, 0x000a'1848'97c3'63c4ULL}),
+ Frac128({0xe6f8'52f9'914f'2d9aULL, 0x0000'ffe5'fe2c'4573ULL}),
+ Frac128({0x6056'28a6'5061'b43fULL, 0x0000'162c'0223'a816ULL}),
+ Frac128({0x6b99'3efe'2193'e54cULL, 0x0000'01b5'253d'0671ULL}),
+ Frac128({0x9c76'f49c'ed65'743aULL, 0x0000'001e'4cf8'0b70ULL}),
+ Frac128({0x1ac1'3330'6c4e'3d33ULL, 0x0000'0001'e8ae'6938ULL}),
+ Frac128({0x219f'9904'1f29'5f13ULL, 0x0000'0000'1cd9'af73ULL}),
};
-// TODO: make 128-bit exp version work first, fast 64-bit and Ziv's test version
-// later on.
-// When rolling out the 64-bit fast path, take a look at expf also
-//
-// TODO: Ziv's test + pushing from Frac64 to full Frac128 pipeline on
-// hard-to-round cases Ziv test by extracting the last 12 bits of the result
-// (result & ((1u << 13) - 1)), cast into uint32_t, +4, and then | (1u << 12).
-// If the result is > bound, we have a hard to round case, then dial to Frac128
-// pipeline
-//
-// TODO: current implementation is mostly ported over from expf. Tricky input
-// tests are failing. Debug.
-//
-// TODO: test against CORE-MATH
+// Handle the final roundings of the result. Depends on whether Frac64 or
+// Frac128 being used.
+template <typename TFrac, typename TUInt,
+ cpp::enable_if_t<cpp::is_same<TFrac, Frac64>::value ||
+ cpp::is_same<TFrac, Frac128>::value,
+ int> = 0>
+LIBC_INLINE double exp_handle_rounding(TFrac result_frac, bool is_neg, int d,
+ TUInt e_y_unbiased,
+ [[maybe_unused]] int rounding) {
+ constexpr bool IS_FAST_PATH = cpp::is_same<TFrac, Frac64>::value;
+
+ uint32_t shift_length = IS_FAST_PATH ? 11 : 72;
+ uint64_t leading_one = 0;
+
+ // subnormal
+ if (LIBC_UNLIKELY(is_neg && d >= 0)) {
+ e_y_unbiased = 0;
+ leading_one = uint64_t(1) << (52 - d);
+
+ if (d >= 51) {
+ d -= 2;
+ if constexpr (IS_FAST_PATH)
+ result_frac.val[0] >>= 2;
+ else
+ result_frac.val[1] >>= 2;
+ }
+
+ shift_length += d + 1;
+ }
+
+ auto frac_bits = [&]() -> uint64_t {
+ if constexpr (IS_FAST_PATH)
+ return result_frac.val[0];
+ else
+ return result_frac.val[1];
+ };
+
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ TUInt result =
+ (static_cast<TUInt>(frac_bits() >> shift_length) + (leading_one + 1));
+ result >>= 1;
+ result += static_cast<TUInt>(e_y_unbiased) << 32;
+
+ return cpp::bit_cast<double>(result);
+#else // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (rounding == FE_TONEAREST) {
+ TUInt result =
+ (static_cast<TUInt>(frac_bits() >> shift_length) + (leading_one + 1));
+ result >>= 1;
+ result += static_cast<TUInt>(e_y_unbiased) << 32;
+
+ return cpp::bit_cast<double>(result);
+ }
+
+ TUInt should_round_up = 0;
+
+ if (LIBC_UNLIKELY(rounding == FE_UPWARD)) {
+ uint64_t round_up_mask = (uint64_t(1) << (shift_length + 1)) - 1;
+ should_round_up = static_cast<TUInt>((frac_bits() & round_up_mask) != 0);
+ }
+
+ TUInt result = (static_cast<TUInt>(frac_bits() >> (shift_length + 1)) +
+ should_round_up + (leading_one >> 1));
+ result += static_cast<TUInt>(e_y_unbiased) << 32;
+
+ return cpp::bit_cast<double>(result);
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+}
LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
using FPBits = typename fputil::FPBits<double>;
@@ -161,87 +276,78 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
uint64_t x_s = xbits.get_mantissa();
// The idea of the algorithm below is that:
- // For x = hi + mid + low, with:
+ // For x = 2^(hi + mid + low), with:
// - hi is an integer
- // - hi = N * ln(2)
// - mid * 2^4 is an integer
// - lo is the remainder bits
// Then:
- // exp(x) = 2^N * exp(mid) * exp(lo)
+ // exp(x) = 2^hi * 2^mid * 2^lo
// With this formula:
- // - Multiplying by 2^N is exact and cheap, via adding into the exponent
+ // - Multiplying by 2^hi is exact and cheap, via adding into the exponent
// field
- // - exp(mid) can be calculated via the LUT declared above
- // - exp(lo) ~ 1 + lo + a0 * lo^2 + ...
+ // - 2^mid can be calculated via the LUT declared above
+ // - 2^lo ~ 1 + lo + a0 * lo^2 + ...
// Then we can construct exp(x) pretty easily, as hi, mid, lo bits can be
// separate and be used independently, then we only need to reconstruct in the
// final steps, which makes our life easier.
- //
- // Errors in computing in double-precision
- // TODO
// Range reduction
- // Recall binary64:
+ // binary64 recall:
// 1 sign bit
// 11 exp bits
// 52 sig bits
- // add leading bit = 1
- x_s |= uint64_t(1) << FPBits::FRACTION_LEN;
+ uint64_t x_s_shifted = (x_s << 11) | (uint64_t(1) << 63);
+ // LSB(x_s_frac) = 2^-52
+ Frac64 x_s_frac(x_s_shifted);
- // shift to top 64 bits --> decimal point at hidden bit that we've added
- int x_e_unbiased = static_cast<int>(x_e) - FPBits::EXP_BIAS;
+ constexpr Frac64 INV_LN2_HI = INV_LN2_F128.to_frac64();
+ Frac64 x_ln2 = x_s_frac * INV_LN2_HI;
+ uint64_t x_ln2_bits = x_ln2.val[0];
- // Shift left by 128 - 53 = 75 bits for the fractional part of the double
- // to align to the correct positions inside Frac128
- // LSB(x_s_frac) = 2^-75
- UInt128 x_s_shifted = static_cast<UInt128>(x_s) << 75;
+ uint64_t k = 0;
+ uint64_t frac_bits = 0;
+ int shift = 0;
- // Apply the exponent: shift so the binary point is at the hidden bit
- if (x_e_unbiased > 0) {
- x_s_shifted <<= x_e_unbiased;
- } else if (x_e_unbiased < 0) {
- x_s_shifted >>= -x_e_unbiased;
+ int x_e_unbiased = static_cast<int>(x_e) - FPBits::EXP_BIAS;
+ if (x_e_unbiased >= 0) {
+ shift = 62 - x_e_unbiased;
+ k = x_ln2_bits >> shift;
+ frac_bits = x_ln2_bits << (64 - shift);
+ } else {
+ k = 0;
+ shift = -x_e_unbiased;
+ frac_bits = (shift < 64) ? (x_ln2_bits >> shift) : 0;
}
- // LSB(x_s_frac) = 2^-75
- Frac128 x_s_frac(x_s_shifted);
-
- // LSB(x_ln2) = 2^-74 (product of 2^-75 * 2^127 / 2^128 ~ 2^-75 effectively,
- // but the integer part of x*log2(e) lands in the high word)
- Frac128 x_ln2 = x_s_frac * INV_LN2_F128;
-
- // We use the top 128-bit word of x_ln2 to extract N (the integer part of
- // x * log2(e)). The fractional part drives the polynomial/LUT evaluation.
- constexpr uint64_t FRAC_MASK = (uint64_t(1) << 52) - 1;
- uint64_t x_ln2_hi = x_ln2.val[1];
-
- uint64_t e_y, l2y_r_hi;
- uint32_t e_y_unbiased;
-
// As Frac128 can't store the sign, we need to handle the sign separately:
// - Both branches are computing floor(x * log2(e)).
- // - For negative x, we round up to the next multiple of 2^52, then clear the
- // last 52 bits.
+ // - For negative x, we round up to the next multiple of 2^52, then clear
+ // the last 52 bits.
// - For positive x, we round down (just clear) the last 52 bits.
//
- // Then, l2y_r_hi is the remainder of x * log2(e) after removing the integer
- // part, which is used to look up EXP_MID4 and compute exp(lo) - 1.
+ // Then, l2y_r_hi is the remainder of x * log2(e) after removing the
+ // integer part, which is used to look up EXP_MID4 and compute exp(lo) - 1.
//
// e_y_unbiased is biased exponent field, but already bit-positioned to the
// exponent field of the double representation.
+ uint64_t l2y_r_hi;
+ uint64_t e_y_unbiased;
+
if (LIBC_UNLIKELY(is_neg)) {
- e_y = (x_ln2_hi + FRAC_MASK) & ~FRAC_MASK;
- l2y_r_hi = e_y - x_ln2_hi;
- e_y_unbiased = (FPBits::EXP_BIAS << 20) - static_cast<uint32_t>(e_y >> 43);
+ if (frac_bits != 0) {
+ k = k + 1;
+ l2y_r_hi = ~frac_bits + 1; // 1 - r
+ } else {
+ l2y_r_hi = 0;
+ }
+ e_y_unbiased = (FPBits::EXP_BIAS << 20) - static_cast<uint32_t>(k << 20);
} else {
- e_y = x_ln2_hi & ~FRAC_MASK;
- l2y_r_hi = x_ln2_hi - e_y;
- e_y_unbiased = (FPBits::EXP_BIAS << 20) + static_cast<uint32_t>(e_y >> 43);
+ l2y_r_hi = frac_bits;
+ e_y_unbiased = (FPBits::EXP_BIAS << 20) + static_cast<uint32_t>(k << 20);
}
- uint32_t k = static_cast<uint32_t>(e_y >> 52);
int d = static_cast<int>(k) - FPBits::EXP_BIAS;
// d >= 53 --> k >= 1076
@@ -253,11 +359,53 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
}
// Extract the 4 mid bits (bits [51:48] of l2y_r_hi) for LUT index
- uint16_t x_mid = static_cast<uint16_t>((l2y_r_hi >> 48) & 0xf);
+ uint16_t x_mid = static_cast<uint16_t>((l2y_r_hi >> 60) & 0xf);
// Extract the low 48 bits for the polynomial approximation of exp(lo).
// LSB(x_lo) = 2^-48
- uint64_t x_lo = l2y_r_hi & ((uint64_t(1) << 48) - 1);
+ uint64_t x_lo = l2y_r_hi & ((uint64_t(1) << 60) - 1);
+
+ // Fast path: 64-bit calculations first
+
+ Frac64 x_lo_frac64(x_lo << 4); // aligning LSB to Frac64
+
+ Frac64 p64 =
+ x_lo_frac64 *
+ fputil::polyeval(x_lo_frac64, EXP_COEFFS[0].to_frac64(),
+ EXP_COEFFS[1].to_frac64(), EXP_COEFFS[2].to_frac64(),
+ EXP_COEFFS[3].to_frac64(), EXP_COEFFS[4].to_frac64(),
+ EXP_COEFFS[5].to_frac64(), EXP_COEFFS[6].to_frac64(),
+ EXP_COEFFS[7].to_frac64(), EXP_COEFFS[8].to_frac64(),
+ EXP_COEFFS[9].to_frac64(), EXP_COEFFS[10].to_frac64());
+
+ // With:
+ // - p = 2^lo - 1 --> exp(lo) = p + 1
+ // - mid_val = 2^mid = EXP_MID[x_mid]
+ // We have:
+ // 2^mid * 2^lo = mid_val * (p + 1)
+ // The same applies for both of the 64-bit and 128-bit paths
+ // (Workaround because we're dealing with fractional representation of things)
+ Frac64 mid_val64 = EXP_MID[x_mid].to_frac64();
+ Frac64 result64 = mid_val64 * p64 + mid_val64;
+
+ // Testing
+ uint64_t result64_bits = result64.val[0];
+ constexpr uint64_t LAST_BITS_MASK = ((1u << 13) - 1);
+ uint32_t result64_last_bits =
+ static_cast<uint32_t>(result64_bits & LAST_BITS_MASK);
+ bool is_hard = (result64_last_bits + 4) &
+ LAST_BITS_MASK; // will be != 0 if hard-to-round
+
+ // DEBUG(x, is_neg, rounding, x_e_unbiased, x_s, x_ln2_bits, k, l2y_r_hi,
+ // e_y_unbiased, d, x_mid, x_lo, result64_bits, result64_last_bits,
+ // is_hard);
+
+ // Execute the fast path if possible
+ if (LIBC_LIKELY(!is_hard)) {
+ return exp_handle_rounding(result64, is_neg, d, e_y_unbiased, rounding);
+ }
+
+ // Else, dial to the 128-bit path
// Shift left to move all bits to the high part
// LSB(x_lo_frac) = 2^-128
@@ -270,80 +418,10 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
EXP_COEFFS[8], EXP_COEFFS[9], EXP_COEFFS[10],
EXP_COEFFS[11]);
- // With:
- // - p = exp(lo) - 1 --> exp(lo) = p + 1
- // - mid_val = exp(mid) = EXP_MID4[x_mid]
- // We have:
- // exp(mid) * exp(lo) = mid_val * (p + 1)
- // (Workaround because we're dealing with fractional representation of things)
- Frac128 mid_val = EXP_MID4[x_mid];
+ Frac128 mid_val = EXP_MID[x_mid];
Frac128 result128 = mid_val * p + mid_val;
- // We're computing with errors < worst-case errors, so tie-rounding never
- // happens. Hence, round-to-nearest, tie-to-even is equivalent to
- // round-to-nearest, tie-to-away. Which is what we're implementing below
- // in the following order:
- //
- // 1. Shift so that the rounding bit is at bit-0
- // 2. Add 1 for rounding
- // 3. Perform another shift by 1
- // 4. Depending on the rounding modes, adjust accordingly:
- // a. Add 1 if rounding-up (0 if not)
- // b. Add e_y_unbiased to the result (0 if the result is subnormal)
-
- uint32_t shift_length = 72;
- uint32_t leading_one = 0;
-
- // subnormal
- if (LIBC_UNLIKELY(is_neg && d >= 0)) {
- e_y_unbiased = 0;
- leading_one = 1 << (52 - d);
-
- // In the below shifts, we're shifting by (shift_length + 1) at max, while
- // shift_length is already 72, and if d = 52 --> shift_length + d = 124,
- // and we'll shift by whole 128 bits, which is undefined behavior in C++.
- //
- // So, we'll truncate the last 2 bits.
- if (d >= 51) {
- d -= 2;
- result128.val[1] >>= 2;
- }
-
- shift_length += d + 1;
- }
-
-#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- uint64_t result = (static_cast<uint64_t>(result128.val[1] >> shift_length) +
- (static_cast<uint64_t>(leading_one) + 1));
- result >>= 1;
- result += static_cast<uint64_t>(e_y_unbiased) << 32;
-
- return cpp::bit_cast<double>(result);
-#else // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- if (rounding == FE_TONEAREST) {
- uint64_t result = (static_cast<uint64_t>(result128.val[1] >> shift_length) +
- (static_cast<uint64_t>(leading_one) + 1));
- result >>= 1;
- result += static_cast<uint64_t>(e_y_unbiased) << 32;
-
- return cpp::bit_cast<double>(result);
- }
-
- uint64_t should_round_up = 0;
-
- if (LIBC_UNLIKELY(rounding == FE_UPWARD)) {
- uint64_t round_up_mask = (uint64_t(1) << (shift_length + 1)) - 1;
- should_round_up =
- static_cast<uint64_t>((result128.val[1] & round_up_mask) != 0);
- }
-
- uint64_t result =
- (static_cast<uint64_t>(result128.val[1] >> (shift_length + 1)) +
- should_round_up + (static_cast<uint64_t>(leading_one) >> 1));
- result += static_cast<uint64_t>(e_y_unbiased) << 32;
-
- return cpp::bit_cast<double>(result);
-#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ return exp_handle_rounding(result128, is_neg, d, e_y_unbiased, rounding);
}
} // namespace static_rounding
diff --git a/libc/src/math/generic/exp.cpp b/libc/src/math/generic/exp.cpp
index 45f105aa25a90..4443c3bd23c86 100644
--- a/libc/src/math/generic/exp.cpp
+++ b/libc/src/math/generic/exp.cpp
@@ -1,26 +1,20 @@
-//===-- Double-precision e^x function -------------------------------------===//
+//===----------------------------------------------------------------------===//
//
// 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 double-precision e^x function
+///
+//===----------------------------------------------------------------------===//
#include "src/math/exp.h"
-#include "shared/math/static_rounding/exp.h"
#include "src/__support/math/exp.h"
namespace LIBC_NAMESPACE_DECL {
-LLVM_LIBC_FUNCTION(double, exp, (double x)) {
-// TODO: update to follow the new structure (following
-// https://github.com/llvm/llvm-project/pull/224735)
-#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::exp(x, FE_TONEAREST);
-#else
- return math::exp(x);
-#endif
-}
+LLVM_LIBC_FUNCTION(double, exp, (double x)) { return math::exp(x); }
} // namespace LIBC_NAMESPACE_DECL
>From 0730ae77b92157609b77b9043247b8b36d653af8 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Sun, 27 Sep 2026 15:31:53 +0700
Subject: [PATCH 3/6] test: updated InDoubleRange test & refactor: refactored
fast/accurate paths
still having big inaccuracies, still debugging
---
libc/src/__support/math/exp_integer_eval.h | 175 ++++++++++++------
libc/test/UnitTest/FPMatcher.h | 10 +
.../src/math/exp_static_rounding_test.cpp | 45 ++---
3 files changed, 149 insertions(+), 81 deletions(-)
diff --git a/libc/src/__support/math/exp_integer_eval.h b/libc/src/__support/math/exp_integer_eval.h
index ae06e9ed44fec..ebf29118006f6 100644
--- a/libc/src/__support/math/exp_integer_eval.h
+++ b/libc/src/__support/math/exp_integer_eval.h
@@ -80,6 +80,7 @@
// std::cout << "===== DEBUG =====\n";
// size_t i = 0;
// (debug_value(parts[i++], values), ...);
+// std::cout << "=================\n\n";
// }
//
// #define DEBUG(...) debug_all(#__VA_ARGS__, __VA_ARGS__)
@@ -96,6 +97,7 @@ namespace static_rounding {
// truncating to all-0 in all rounding modes.
// TODO: test against CORE-MATH
// TODO: refactor to follow the new structure; dedupe codes
+// TODO: LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// print(2+round(1/log(2), 128, RN));
// LSB(INV_LN2) = 2^-127
@@ -148,25 +150,26 @@ LIBC_INLINE_VAR constexpr Frac128 EXP_COEFFS[] = {
Frac128({0x219f'9904'1f29'5f13ULL, 0x0000'0000'1cd9'af73ULL}),
};
-// Handle the final roundings of the result. Depends on whether Frac64 or
-// Frac128 being used.
+// Round the fractional result and combine it with its exponent.
template <typename TFrac, typename TUInt,
cpp::enable_if_t<cpp::is_same<TFrac, Frac64>::value ||
cpp::is_same<TFrac, Frac128>::value,
int> = 0>
LIBC_INLINE double exp_handle_rounding(TFrac result_frac, bool is_neg, int d,
- TUInt e_y_unbiased,
+ TUInt e_y,
[[maybe_unused]] int rounding) {
constexpr bool IS_FAST_PATH = cpp::is_same<TFrac, Frac64>::value;
- uint32_t shift_length = IS_FAST_PATH ? 11 : 72;
+ uint32_t shift_length = 11;
uint64_t leading_one = 0;
// subnormal
if (LIBC_UNLIKELY(is_neg && d >= 0)) {
- e_y_unbiased = 0;
+ e_y = 0;
leading_one = uint64_t(1) << (52 - d);
+ // Truncate the last 2 bits to avoid undefined behavior when shifting by 64
+ // bits
if (d >= 51) {
d -= 2;
if constexpr (IS_FAST_PATH)
@@ -189,7 +192,7 @@ LIBC_INLINE double exp_handle_rounding(TFrac result_frac, bool is_neg, int d,
TUInt result =
(static_cast<TUInt>(frac_bits() >> shift_length) + (leading_one + 1));
result >>= 1;
- result += static_cast<TUInt>(e_y_unbiased) << 32;
+ result += static_cast<TUInt>(e_y) << 32;
return cpp::bit_cast<double>(result);
#else // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
@@ -197,7 +200,7 @@ LIBC_INLINE double exp_handle_rounding(TFrac result_frac, bool is_neg, int d,
TUInt result =
(static_cast<TUInt>(frac_bits() >> shift_length) + (leading_one + 1));
result >>= 1;
- result += static_cast<TUInt>(e_y_unbiased) << 32;
+ result += static_cast<TUInt>(e_y) << 32;
return cpp::bit_cast<double>(result);
}
@@ -206,17 +209,78 @@ LIBC_INLINE double exp_handle_rounding(TFrac result_frac, bool is_neg, int d,
if (LIBC_UNLIKELY(rounding == FE_UPWARD)) {
uint64_t round_up_mask = (uint64_t(1) << (shift_length + 1)) - 1;
- should_round_up = static_cast<TUInt>((frac_bits() & round_up_mask) != 0);
+ bool has_remainder = (frac_bits() & round_up_mask) != 0;
+ if constexpr (!IS_FAST_PATH) {
+ has_remainder = has_remainder || (result_frac.val[0] != 0);
+ }
+ should_round_up = static_cast<TUInt>(has_remainder);
}
TUInt result = (static_cast<TUInt>(frac_bits() >> (shift_length + 1)) +
should_round_up + (leading_one >> 1));
- result += static_cast<TUInt>(e_y_unbiased) << 32;
+ result += static_cast<TUInt>(e_y) << 32;
return cpp::bit_cast<double>(result);
#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
}
+LIBC_INLINE double exp_accurate_path(uint64_t x_s_shifted, int x_e_unbiased,
+ bool is_neg, int rounding) {
+ using FPBits = typename fputil::FPBits<double>;
+
+ // Recalculate everything in 128-bit precision, with the same idea as the
+ // 64-bit path.
+
+ Frac128 x_s_frac({0, x_s_shifted});
+ Frac128 x_ln2 = x_s_frac * INV_LN2_F128;
+
+ uint64_t k = 0;
+ Frac128 l2y_r;
+ if (x_e_unbiased >= -1) {
+ int shift = 62 - x_e_unbiased;
+ k = (x_ln2 >> (64 + shift)).val[0];
+ l2y_r = x_ln2 << (64 - shift);
+ } else {
+ int shift = -x_e_unbiased - 2;
+ l2y_r = (shift < 128) ? (x_ln2 >> shift) : Frac128{};
+ }
+
+ if (LIBC_UNLIKELY(is_neg)) {
+ if (l2y_r.val[0] != 0 || l2y_r.val[1] != 0) {
+ ++k;
+ l2y_r = ~l2y_r + Frac128(1);
+ }
+ }
+
+ uint64_t e_y;
+ if (is_neg)
+ e_y = (FPBits::EXP_BIAS << 20) - static_cast<uint32_t>(k << 20);
+ else
+ e_y = (FPBits::EXP_BIAS << 20) + static_cast<uint32_t>(k << 20);
+
+ int d = static_cast<int>(k) - FPBits::EXP_BIAS;
+
+ if (LIBC_UNLIKELY(is_neg && d >= 53))
+ return 0.0;
+
+ uint16_t x_mid = static_cast<uint16_t>((l2y_r.val[1] >> 60) & 0xf);
+
+ Frac128 x_lo_frac = l2y_r;
+ x_lo_frac.val[1] &= 0x0fff'ffff'ffff'ffffULL;
+
+ Frac128 p =
+ x_lo_frac * fputil::polyeval(x_lo_frac, EXP_COEFFS[0], EXP_COEFFS[1],
+ EXP_COEFFS[2], EXP_COEFFS[3], EXP_COEFFS[4],
+ EXP_COEFFS[5], EXP_COEFFS[6], EXP_COEFFS[7],
+ EXP_COEFFS[8], EXP_COEFFS[9], EXP_COEFFS[10],
+ EXP_COEFFS[11]);
+
+ Frac128 mid_val = EXP_MID[x_mid];
+ Frac128 result = mid_val * p + mid_val;
+
+ return exp_handle_rounding(result, is_neg, d, e_y, rounding);
+}
+
LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
using FPBits = typename fputil::FPBits<double>;
FPBits xbits(x);
@@ -311,29 +375,29 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
int shift = 0;
int x_e_unbiased = static_cast<int>(x_e) - FPBits::EXP_BIAS;
- if (x_e_unbiased >= 0) {
+ if (x_e_unbiased >= -1) {
shift = 62 - x_e_unbiased;
k = x_ln2_bits >> shift;
frac_bits = x_ln2_bits << (64 - shift);
} else {
k = 0;
- shift = -x_e_unbiased;
+ shift = -x_e_unbiased - 2;
frac_bits = (shift < 64) ? (x_ln2_bits >> shift) : 0;
}
- // As Frac128 can't store the sign, we need to handle the sign separately:
+ // As Frac64/Frac128 can't store the sign, we need to handle the sign
+ // separately:
// - Both branches are computing floor(x * log2(e)).
// - For negative x, we round up to the next multiple of 2^52, then clear
// the last 52 bits.
// - For positive x, we round down (just clear) the last 52 bits.
//
// Then, l2y_r_hi is the remainder of x * log2(e) after removing the
- // integer part, which is used to look up EXP_MID4 and compute exp(lo) - 1.
+ // integer part, which is used to look up EXP_MID and compute 2^lo.
//
- // e_y_unbiased is biased exponent field, but already bit-positioned to the
- // exponent field of the double representation.
+ // e_y is the biased exponent field, positioned for the final bit assembly.
uint64_t l2y_r_hi;
- uint64_t e_y_unbiased;
+ uint64_t e_y;
if (LIBC_UNLIKELY(is_neg)) {
if (frac_bits != 0) {
@@ -342,10 +406,10 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
} else {
l2y_r_hi = 0;
}
- e_y_unbiased = (FPBits::EXP_BIAS << 20) - static_cast<uint32_t>(k << 20);
+ e_y = (FPBits::EXP_BIAS << 20) - static_cast<uint32_t>(k << 20);
} else {
l2y_r_hi = frac_bits;
- e_y_unbiased = (FPBits::EXP_BIAS << 20) + static_cast<uint32_t>(k << 20);
+ e_y = (FPBits::EXP_BIAS << 20) + static_cast<uint32_t>(k << 20);
}
int d = static_cast<int>(k) - FPBits::EXP_BIAS;
@@ -358,20 +422,22 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
return 0.0;
}
- // Extract the 4 mid bits (bits [51:48] of l2y_r_hi) for LUT index
+ // Extract the top 4 fractional bits for the LUT index.
uint16_t x_mid = static_cast<uint16_t>((l2y_r_hi >> 60) & 0xf);
- // Extract the low 48 bits for the polynomial approximation of exp(lo).
- // LSB(x_lo) = 2^-48
+ // The remaining 60 bits are the polynomial input.
uint64_t x_lo = l2y_r_hi & ((uint64_t(1) << 60) - 1);
// Fast path: 64-bit calculations first
- Frac64 x_lo_frac64(x_lo << 4); // aligning LSB to Frac64
+ // Don't << 4, as the polynomial approximation is correct in range [0, 1/16],
+ // and x_lo is already in that range.
+ // LSB(x_lo_frac) = 2^-64
+ Frac64 x_lo_frac(x_lo);
- Frac64 p64 =
- x_lo_frac64 *
- fputil::polyeval(x_lo_frac64, EXP_COEFFS[0].to_frac64(),
+ Frac64 p =
+ x_lo_frac *
+ fputil::polyeval(x_lo_frac, EXP_COEFFS[0].to_frac64(),
EXP_COEFFS[1].to_frac64(), EXP_COEFFS[2].to_frac64(),
EXP_COEFFS[3].to_frac64(), EXP_COEFFS[4].to_frac64(),
EXP_COEFFS[5].to_frac64(), EXP_COEFFS[6].to_frac64(),
@@ -379,49 +445,46 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
EXP_COEFFS[9].to_frac64(), EXP_COEFFS[10].to_frac64());
// With:
- // - p = 2^lo - 1 --> exp(lo) = p + 1
+ // - p = 2^lo - 1 --> 2^lo = p + 1
// - mid_val = 2^mid = EXP_MID[x_mid]
// We have:
// 2^mid * 2^lo = mid_val * (p + 1)
// The same applies for both of the 64-bit and 128-bit paths
// (Workaround because we're dealing with fractional representation of things)
- Frac64 mid_val64 = EXP_MID[x_mid].to_frac64();
- Frac64 result64 = mid_val64 * p64 + mid_val64;
+ Frac64 mid_val = EXP_MID[x_mid].to_frac64();
+ Frac64 result = mid_val * p + mid_val;
+
+ uint64_t result_bits = result.val[0];
- // Testing
- uint64_t result64_bits = result64.val[0];
- constexpr uint64_t LAST_BITS_MASK = ((1u << 13) - 1);
- uint32_t result64_last_bits =
- static_cast<uint32_t>(result64_bits & LAST_BITS_MASK);
- bool is_hard = (result64_last_bits + 4) &
- LAST_BITS_MASK; // will be != 0 if hard-to-round
+ // Rounding test
+ constexpr uint32_t LAST_BITS = 12;
+ constexpr uint32_t ROUNDING_ERROR = 4;
+ uint32_t result_last_bits =
+ static_cast<uint32_t>(result_bits & ((1u << LAST_BITS) - 1));
+ bool is_hard;
+ if (rounding == FE_TONEAREST) {
+ uint32_t rounded_lo =
+ (result_last_bits + (1u << (LAST_BITS - 1)) - ROUNDING_ERROR) >>
+ LAST_BITS;
+ uint32_t rounded_hi =
+ (result_last_bits + (1u << (LAST_BITS - 1)) + ROUNDING_ERROR) >>
+ LAST_BITS;
+ is_hard = rounded_lo != rounded_hi;
+ } else {
+ is_hard = result_last_bits <= ROUNDING_ERROR ||
+ result_last_bits >= (1u << LAST_BITS) - ROUNDING_ERROR;
+ }
// DEBUG(x, is_neg, rounding, x_e_unbiased, x_s, x_ln2_bits, k, l2y_r_hi,
- // e_y_unbiased, d, x_mid, x_lo, result64_bits, result64_last_bits,
+ // e_y, d, x_mid, x_lo, result64_bits, result_last_bits,
// is_hard);
- // Execute the fast path if possible
- if (LIBC_LIKELY(!is_hard)) {
- return exp_handle_rounding(result64, is_neg, d, e_y_unbiased, rounding);
+ if (false && LIBC_LIKELY(!is_hard)) {
+ return exp_handle_rounding(result, is_neg, d, e_y, rounding);
}
- // Else, dial to the 128-bit path
-
- // Shift left to move all bits to the high part
- // LSB(x_lo_frac) = 2^-128
- Frac128 x_lo_frac(static_cast<UInt128>(x_lo) << 80);
-
- Frac128 p =
- x_lo_frac * fputil::polyeval(x_lo_frac, EXP_COEFFS[0], EXP_COEFFS[1],
- EXP_COEFFS[2], EXP_COEFFS[3], EXP_COEFFS[4],
- EXP_COEFFS[5], EXP_COEFFS[6], EXP_COEFFS[7],
- EXP_COEFFS[8], EXP_COEFFS[9], EXP_COEFFS[10],
- EXP_COEFFS[11]);
-
- Frac128 mid_val = EXP_MID[x_mid];
- Frac128 result128 = mid_val * p + mid_val;
-
- return exp_handle_rounding(result128, is_neg, d, e_y_unbiased, rounding);
+ // Dial back to 128-bit path for hard-to-round cases
+ return exp_accurate_path(x_s_shifted, x_e_unbiased, is_neg, rounding);
}
} // namespace static_rounding
diff --git a/libc/test/UnitTest/FPMatcher.h b/libc/test/UnitTest/FPMatcher.h
index 314d8f4f1ef1e..cf09e8d75c66e 100644
--- a/libc/test/UnitTest/FPMatcher.h
+++ b/libc/test/UnitTest/FPMatcher.h
@@ -170,6 +170,16 @@ template <TestCond C, typename T> FPMatcher<T, C> getMatcher(T expectedValue) {
return FPMatcher<T, C>(expectedValue);
}
+template <typename T>
+typename fputil::FPBits<T>::StorageType ulp_distance(T x, T y) {
+ using FPBits = fputil::FPBits<T>;
+ using StorageType = typename FPBits::StorageType;
+
+ StorageType x_bits = FPBits(x).uintval();
+ StorageType y_bits = FPBits(y).uintval();
+ return x_bits > y_bits ? x_bits - y_bits : y_bits - x_bits;
+}
+
template <TestCond C, typename T>
CFPMatcher<T, C> getMatcherComplex(T expectedValue) {
return CFPMatcher<T, C>(expectedValue);
diff --git a/libc/test/src/math/exp_static_rounding_test.cpp b/libc/test/src/math/exp_static_rounding_test.cpp
index f3fc4a5796ee2..22438dd7f57db 100644
--- a/libc/test/src/math/exp_static_rounding_test.cpp
+++ b/libc/test/src/math/exp_static_rounding_test.cpp
@@ -76,8 +76,6 @@ TEST_F(LlvmLibcExpStaticRoundingTest, TrickyInputs) {
}
}
-// TODO: double has a very big range. Define test functions that allows for
-// finding max ULPs over a range. EXPECT_* just isn't enough for this.
TEST_F(LlvmLibcExpStaticRoundingTest, InDoubleRange) {
constexpr uint64_t COUNT = 1'231;
constexpr uint64_t START = FPBits(0.25).uintval();
@@ -85,47 +83,44 @@ TEST_F(LlvmLibcExpStaticRoundingTest, InDoubleRange) {
constexpr uint64_t STEP = (STOP - START) / COUNT;
auto test = [&](RoundingMode rounding) {
+ ForceRoundingMode __r(rounding);
+ if (!__r.success)
+ return;
const int fenv_rounding = get_fe_rounding(rounding);
uint64_t fails = 0;
uint64_t count = 0;
uint64_t cc = 0;
- double mx, mr = 0.0;
- double tol = 0.5;
+ uint64_t max_ulp = 0;
+ double me = 0.0, mr = 0.0;
for (uint64_t i = 0, v = START; i <= COUNT; ++i, v += STEP) {
double x = FPBits(v).get_val();
if (FPBits(v).is_nan() || FPBits(v).is_inf() || x < 0.0)
continue;
- double result = math::exp(x);
+ double expected = math::exp(x);
+ double result = static_rounding::exp(x, fenv_rounding);
++cc;
- if (FPBits(result).is_nan() || FPBits(result).is_inf())
+ if (FPBits(expected).is_nan() || FPBits(expected).is_inf() ||
+ FPBits(result).is_nan() || FPBits(result).is_inf())
continue;
++count;
- // ASSERT_MPFR_MATCH(mpfr::Operation::Log, x, result, 0.5);
- // if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Exp, x, result,
- // TOLERANCE + 0.5, rounding_mode))
- // {
- // ++fails;
- // while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Exp, x,
- // result, tol,
- // rounding_mode)) {
- // mx = x;
- // mr = result;
-
- // if (tol > 1000.0)
- // break;
-
- // tol *= 2.0;
- // }
- // }
+ uint64_t ulp = LIBC_NAMESPACE::testing::ulp_distance(expected, result);
+ if (ulp != 0) {
+ ++fails;
+ if (ulp > max_ulp) {
+ max_ulp = ulp;
+ me = expected;
+ mr = result;
+ }
+ }
}
if (fails) {
tlog << " Statically rounded exp failed: " << fails << "/" << count << "/"
<< cc << " tests.\n";
- tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
- // EXPECT_MPFR_MATCH(mpfr::Operation::Exp, mx, mr, 0.5, rounding_mode);
+ tlog << " Max ULPs is: " << max_ulp << ".\n";
+ EXPECT_FP_EQ(me, mr);
}
};
>From 4a2dd08024261f2e838b86d0be4d5de8d8e63e71 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Sun, 27 Sep 2026 17:03:46 +0700
Subject: [PATCH 4/6] fix: all failing tests now pass
---
libc/src/__support/math/exp_integer_eval.h | 11 +++++++----
1 file changed, 7 insertions(+), 4 deletions(-)
diff --git a/libc/src/__support/math/exp_integer_eval.h b/libc/src/__support/math/exp_integer_eval.h
index ebf29118006f6..b328bb3ce2c3b 100644
--- a/libc/src/__support/math/exp_integer_eval.h
+++ b/libc/src/__support/math/exp_integer_eval.h
@@ -21,6 +21,7 @@
#include "src/__support/CPP/type_traits/is_same.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/frac128.h"
#include "src/__support/integer_literals.h"
#include "src/__support/macros/config.h"
@@ -184,8 +185,10 @@ LIBC_INLINE double exp_handle_rounding(TFrac result_frac, bool is_neg, int d,
auto frac_bits = [&]() -> uint64_t {
if constexpr (IS_FAST_PATH)
return result_frac.val[0];
- else
- return result_frac.val[1];
+ else {
+ // Discarding the leading bit
+ return (result_frac.val[1] << 1) | (result_frac.val[0] >> 63);
+ }
};
#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
@@ -276,7 +279,7 @@ LIBC_INLINE double exp_accurate_path(uint64_t x_s_shifted, int x_e_unbiased,
EXP_COEFFS[11]);
Frac128 mid_val = EXP_MID[x_mid];
- Frac128 result = mid_val * p + mid_val;
+ Frac128 result = fputil::multiply_add(mid_val, p, mid_val);
return exp_handle_rounding(result, is_neg, d, e_y, rounding);
}
@@ -452,7 +455,7 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
// The same applies for both of the 64-bit and 128-bit paths
// (Workaround because we're dealing with fractional representation of things)
Frac64 mid_val = EXP_MID[x_mid].to_frac64();
- Frac64 result = mid_val * p + mid_val;
+ Frac64 result = fputil::multiply_add(mid_val, p, mid_val);
uint64_t result_bits = result.val[0];
>From e22266afc95317e7bdb392afca5d6dfb3dcea238 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Sun, 27 Sep 2026 18:07:47 +0700
Subject: [PATCH 5/6] fix: accurately enable fast path
was disabled before
---
libc/src/__support/math/exp_integer_eval.h | 112 ++++++---------------
1 file changed, 30 insertions(+), 82 deletions(-)
diff --git a/libc/src/__support/math/exp_integer_eval.h b/libc/src/__support/math/exp_integer_eval.h
index b328bb3ce2c3b..d39eaeff381ce 100644
--- a/libc/src/__support/math/exp_integer_eval.h
+++ b/libc/src/__support/math/exp_integer_eval.h
@@ -23,68 +23,9 @@
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/frac128.h"
-#include "src/__support/integer_literals.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 <iomanip>
-// #include <iostream>
-// #include <string>
-// #include <type_traits>
-// #include <vector>
-//
-// std::vector<std::string> split_top_level(const std::string &s) {
-// std::vector<std::string> parts;
-// int depth = 0;
-// size_t start = 0;
-// for (size_t i = 0; i < s.size(); ++i) {
-// char c = s[i];
-// if (c == '(' || c == '[' || c == '{')
-// ++depth;
-// else if (c == ')' || c == ']' || c == '}')
-// --depth;
-// else if (c == ',' && depth == 0) {
-// parts.push_back(s.substr(start, i - start));
-// start = i + 1;
-// }
-// }
-// parts.push_back(s.substr(start));
-// for (auto &p : parts) {
-// size_t a = p.find_first_not_of(" \t\n");
-// size_t b = p.find_last_not_of(" \t\n");
-// p = (a == std::string::npos) ? "" : p.substr(a, b - a + 1);
-// }
-// return parts;
-// }
-//
-// template <typename T>
-// void debug_value(const std::string &name, const T &value) {
-// using U = std::remove_cv_t<std::remove_reference_t<T>>;
-// std::cout << std::left << std::setw(24) << name << " = ";
-// if constexpr (std::is_floating_point_v<U>) {
-// std::cout << std::defaultfloat << std::setw(24) << value << " ("
-// << std::hexfloat << std::right << std::setw(24) << value
-// << std::defaultfloat << std::left << ")";
-// } else if constexpr (std::is_integral_v<U>) {
-// std::cout << std::dec << std::setw(24) << value << " (" << std::right
-// << std::setw(24) << std::showbase << std::hex << value
-// << std::noshowbase << std::dec << ")";
-// } else {
-// std::cout << value;
-// }
-// std::cout << '\n';
-// }
-//
-// template <typename... Ts> void debug_all(const char *names, Ts &&...values) {
-// auto parts = split_top_level(names);
-// std::cout << "===== DEBUG =====\n";
-// size_t i = 0;
-// (debug_value(parts[i++], values), ...);
-// std::cout << "=================\n\n";
-// }
-//
-// #define DEBUG(...) debug_all(#__VA_ARGS__, __VA_ARGS__)
namespace LIBC_NAMESPACE_DECL {
@@ -94,11 +35,8 @@ namespace math {
namespace static_rounding {
-// TODO: tricky input tests are failing. Debug. Getting the last bits
-// truncating to all-0 in all rounding modes.
// TODO: test against CORE-MATH
// TODO: refactor to follow the new structure; dedupe codes
-// TODO: LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// print(2+round(1/log(2), 128, RN));
// LSB(INV_LN2) = 2^-127
@@ -151,6 +89,21 @@ LIBC_INLINE_VAR constexpr Frac128 EXP_COEFFS[] = {
Frac128({0x219f'9904'1f29'5f13ULL, 0x0000'0000'1cd9'af73ULL}),
};
+// 64-bit polynomial approximation of 2^x coefficients generated with Sollya:
+// > P = fpminimax(2^x, 6, [|1, 64...|], [0, 1/16], absolute, fixed);
+// Store the fractional part of the coefficients below
+// > dirtyinfnorm(2^x - P(x), [0, 1/16]);
+// 0x1.2ac9...p-57
+// This is different from EXPF_COEFFS: EXPF_COEFFS is for approximating for x in
+// range of [0, 1]
+// LSB(EXP_COEFFS[i]) = 2^-64
+LIBC_INLINE_VAR constexpr Frac64 EXP_64_COEFFS[] = {
+ // degree-0 = 1, add back afterwards to reduce calc ops
+ Frac64(0xb172'17f7'd1cd'3e3cULL), Frac64(0x3d7f'7bff'0838'4d4aULL),
+ Frac64(0x0e35'846a'6f46'd462ULL), Frac64(0x0276'55a0'd536'1dd9ULL),
+ Frac64(0x0057'5d3d'4b45'0056ULL), Frac64(0x000a'504b'13fe'e008ULL),
+};
+
// Round the fractional result and combine it with its exponent.
template <typename TFrac, typename TUInt,
cpp::enable_if_t<cpp::is_same<TFrac, Frac64>::value ||
@@ -182,13 +135,12 @@ LIBC_INLINE double exp_handle_rounding(TFrac result_frac, bool is_neg, int d,
shift_length += d + 1;
}
+ // Get the bits, discarding the leading bit
auto frac_bits = [&]() -> uint64_t {
if constexpr (IS_FAST_PATH)
- return result_frac.val[0];
- else {
- // Discarding the leading bit
+ return result_frac.val[0] << 1;
+ else
return (result_frac.val[1] << 1) | (result_frac.val[0] >> 63);
- }
};
#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
@@ -438,14 +390,10 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
// LSB(x_lo_frac) = 2^-64
Frac64 x_lo_frac(x_lo);
- Frac64 p =
- x_lo_frac *
- fputil::polyeval(x_lo_frac, EXP_COEFFS[0].to_frac64(),
- EXP_COEFFS[1].to_frac64(), EXP_COEFFS[2].to_frac64(),
- EXP_COEFFS[3].to_frac64(), EXP_COEFFS[4].to_frac64(),
- EXP_COEFFS[5].to_frac64(), EXP_COEFFS[6].to_frac64(),
- EXP_COEFFS[7].to_frac64(), EXP_COEFFS[8].to_frac64(),
- EXP_COEFFS[9].to_frac64(), EXP_COEFFS[10].to_frac64());
+ Frac64 p = x_lo_frac * fputil::polyeval(x_lo_frac, EXP_64_COEFFS[0],
+ EXP_64_COEFFS[1], EXP_64_COEFFS[2],
+ EXP_64_COEFFS[3], EXP_64_COEFFS[4],
+ EXP_64_COEFFS[5]);
// With:
// - p = 2^lo - 1 --> 2^lo = p + 1
@@ -457,11 +405,15 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
Frac64 mid_val = EXP_MID[x_mid].to_frac64();
Frac64 result = fputil::multiply_add(mid_val, p, mid_val);
- uint64_t result_bits = result.val[0];
+ uint64_t result_bits = result.val[0] << 1;
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return exp_handle_rounding(result, is_neg, d, e_y, rounding);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
// Rounding test
constexpr uint32_t LAST_BITS = 12;
- constexpr uint32_t ROUNDING_ERROR = 4;
+ constexpr uint32_t ROUNDING_ERROR = 0x800;
uint32_t result_last_bits =
static_cast<uint32_t>(result_bits & ((1u << LAST_BITS) - 1));
bool is_hard;
@@ -478,11 +430,7 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
result_last_bits >= (1u << LAST_BITS) - ROUNDING_ERROR;
}
- // DEBUG(x, is_neg, rounding, x_e_unbiased, x_s, x_ln2_bits, k, l2y_r_hi,
- // e_y, d, x_mid, x_lo, result64_bits, result_last_bits,
- // is_hard);
-
- if (false && LIBC_LIKELY(!is_hard)) {
+ if (LIBC_LIKELY(!is_hard)) {
return exp_handle_rounding(result, is_neg, d, e_y, rounding);
}
>From ce0a3e33ed5d29221069030546cde0f4b2c61a20 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Sun, 27 Sep 2026 18:20:53 +0700
Subject: [PATCH 6/6] chore: refactor
---
libc/src/__support/math/CMakeLists.txt | 12 +-
.../__support/math/exp_integer_constants.h | 128 ++++++++++++++++++
libc/src/__support/math/exp_integer_eval.h | 84 +-----------
libc/src/__support/math/expf_integer_eval.h | 33 +----
4 files changed, 149 insertions(+), 108 deletions(-)
create mode 100644 libc/src/__support/math/exp_integer_constants.h
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index c534dcfed0fbc..5d3ada1e243b3 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -3071,6 +3071,14 @@ add_header_library(
libc.src.__support.macros.config
)
+add_header_library(
+ exp_integer_constants
+ HDRS
+ exp_integer_constants.h
+ DEPENDS
+ libc.src.__support.macros.config
+)
+
add_header_library(
exp2f_float_utils
HDRS
@@ -3137,7 +3145,7 @@ add_header_library(
HDRS
expf_integer_eval.h
DEPENDS
- .exp_constants
+ .exp_integer_constants
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.fp_bits
libc.src.__support.macros.config
@@ -4204,7 +4212,7 @@ add_header_library(
HDRS
exp_integer_eval.h
DEPENDS
- .exp_constants
+ .exp_integer_constants
libc.src.__support.CPP.bit
libc.src.__support.FPUtil.fp_bits
libc.src.__support.macros.config
diff --git a/libc/src/__support/math/exp_integer_constants.h b/libc/src/__support/math/exp_integer_constants.h
new file mode 100644
index 0000000000000..467595026fb4d
--- /dev/null
+++ b/libc/src/__support/math/exp_integer_constants.h
@@ -0,0 +1,128 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 look-up tables for integer-only, statically rounded
+/// exp*(x) functions
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_CONSTANTS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_CONSTANTS_H
+
+#include "src/__support/frac128.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace shared {
+
+namespace math {
+
+namespace static_rounding {
+
+// print(2+round(1/log(2), 128, RN));
+// LSB(INV_LN2) = 2^-127
+LIBC_INLINE_VAR constexpr Frac128 INV_LN2_F128 =
+ Frac128({0xbe87'fed0'691d'3e89ULL, 0xb8aa'3b29'5c17'f0bbULL});
+
+// 2^x for x from 0 to (0b1111/2^4) = 15/16
+// > for i from 0 to 15 do {
+// print(1+round(2^(i/16), 128, RN));
+// };
+// LSB(EXP_MID4[i]) = 2^-127
+LIBC_INLINE_VAR constexpr Frac128 EXP_MID[] = {
+ Frac128({0x0000'0000'0000'0000ULL, 0x8000'0000'0000'0000ULL}),
+ Frac128({0xc5c9'5b8c'2154'c1b2ULL, 0x85aa'c367'cc48'7b14ULL}),
+ Frac128({0xfbe4'6287'58a5'3c90ULL, 0x8b95'c1e3'ea8b'd6e6ULL}),
+ Frac128({0x0fd6'd8e0'ae5a'c9d8ULL, 0x91c3'd373'ab11'c336ULL}),
+ Frac128({0x46ad'2318'2e42'f6f6ULL, 0x9837'f051'8db8'a96fULL}),
+ Frac128({0xa091'1f09'ebb9'fdd1ULL, 0x9ef5'3260'91a1'11adULL}),
+ Frac128({0x1cbd'7f62'1710'701bULL, 0xa5fe'd6a9'b151'38eaULL}),
+ Frac128({0x4980'a8c8'f59a'2ec4ULL, 0xad58'3eea'42a1'4ac6ULL}),
+ Frac128({0x597d'89b3'754a'be9fULL, 0xb504'f333'f9de'6484ULL}),
+ Frac128({0xa881'1fb6'6d0f'af7aULL, 0xbd08'a39f'580c'36beULL}),
+ Frac128({0x3e2a'd0c9'64dd'9f37ULL, 0xc567'2a11'5506'daddULL}),
+ Frac128({0xe235'838f'95f2'c6edULL, 0xce24'8c15'1f84'80e3ULL}),
+ Frac128({0x39a6'8bb9'902d'3fdeULL, 0xd744'fcca'd69d'6af4ULL}),
+ Frac128({0x0658'9504'8dd3'33caULL, 0xe0cc'deec'2a94'e111ULL}),
+ Frac128({0xd02d'75b3'706e'54fbULL, 0xeac0'c6e7'dd24'392eULL}),
+ Frac128({0x7b9d'0c7a'ed98'0fc3ULL, 0xf525'7d15'2486'cc2cULL}),
+};
+
+// 128-bit polynomial approximation of 2^x coefficients generated with Sollya:
+// > P = fpminimax(2^x, 12, [|1, 128...|], [0, 1/16], absolute, fixed);
+// Store the fractional part of the coefficients below
+// > dirtyinfnorm(2^x - P(x), [0, 1/16]);
+// 0x1.b328...p-117
+// LSB(EXPF_COEFFS[i]) = 2^-128
+LIBC_INLINE_VAR constexpr Frac128 EXP_COEFFS[] = {
+ // degree-0 = 1, add back afterwards to reduce calc ops
+ Frac128({0xc9e3'b398'033f'0902ULL, 0xb172'17f7'd1cf'79abULL}),
+ Frac128({0xde2d'60e0'866c'6365ULL, 0x3d7f'7bff'058b'1d50ULL}),
+ Frac128({0x99d3'ad04'd47d'0efeULL, 0x0e35'846b'8250'5fc5ULL}),
+ Frac128({0x399a'b423'1644'1b55ULL, 0x0276'556d'f749'cee5ULL}),
+ Frac128({0x405f'ebef'0f76'011aULL, 0x0057'61ff'9e29'9cc4ULL}),
+ Frac128({0x1ac4'9999'c0ef'b9abULL, 0x000a'1848'97c3'63c4ULL}),
+ Frac128({0xe6f8'52f9'914f'2d9aULL, 0x0000'ffe5'fe2c'4573ULL}),
+ Frac128({0x6056'28a6'5061'b43fULL, 0x0000'162c'0223'a816ULL}),
+ Frac128({0x6b99'3efe'2193'e54cULL, 0x0000'01b5'253d'0671ULL}),
+ Frac128({0x9c76'f49c'ed65'743aULL, 0x0000'001e'4cf8'0b70ULL}),
+ Frac128({0x1ac1'3330'6c4e'3d33ULL, 0x0000'0001'e8ae'6938ULL}),
+ Frac128({0x219f'9904'1f29'5f13ULL, 0x0000'0000'1cd9'af73ULL}),
+};
+
+// 64-bit polynomial approximation of 2^x coefficients generated with Sollya:
+// > P = fpminimax(2^x, 6, [|1, 64...|], [0, 1/16], absolute, fixed);
+// Store the fractional part of the coefficients below
+// > dirtyinfnorm(2^x - P(x), [0, 1/16]);
+// 0x1.2ac9...p-57
+// This is different from EXPF_COEFFS: EXPF_COEFFS is for approximating for x in
+// range of [0, 1]
+// LSB(EXP_COEFFS[i]) = 2^-64
+LIBC_INLINE_VAR constexpr Frac64 EXP_64_COEFFS[] = {
+ // degree-0 = 1, add back afterwards to reduce calc ops
+ Frac64(0xb172'17f7'd1cd'3e3cULL), Frac64(0x3d7f'7bff'0838'4d4aULL),
+ Frac64(0x0e35'846a'6f46'd462ULL), Frac64(0x0276'55a0'd536'1dd9ULL),
+ Frac64(0x0057'5d3d'4b45'0056ULL), Frac64(0x000a'504b'13fe'e008ULL),
+};
+
+// print(2+round(1/log(2), 64, RN));
+// LSB(INV_LN2) = 2^-63
+LIBC_INLINE_VAR constexpr Frac64 INV_LN2_F64 = 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 static_rounding
+
+} // namespace math
+
+} // namespace shared
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_CONSTANTS_H
diff --git a/libc/src/__support/math/exp_integer_eval.h b/libc/src/__support/math/exp_integer_eval.h
index d39eaeff381ce..002e4886cd18f 100644
--- a/libc/src/__support/math/exp_integer_eval.h
+++ b/libc/src/__support/math/exp_integer_eval.h
@@ -15,6 +15,7 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_EVAL_H
+#include "exp_integer_constants.h" // LUTs
#include "hdr/fenv_macros.h"
#include "src/__support/CPP/bit.h"
#include "src/__support/CPP/type_traits/enable_if.h"
@@ -35,75 +36,6 @@ namespace math {
namespace static_rounding {
-// TODO: test against CORE-MATH
-// TODO: refactor to follow the new structure; dedupe codes
-
-// print(2+round(1/log(2), 128, RN));
-// LSB(INV_LN2) = 2^-127
-LIBC_INLINE_VAR constexpr Frac128 INV_LN2_F128 =
- Frac128({0xbe87'fed0'691d'3e89ULL, 0xb8aa'3b29'5c17'f0bbULL});
-
-// 2^x for x from 0 to (0b1111/2^4) = 15/16
-// > for i from 0 to 15 do {
-// print(1+round(2^(i/16), 128, RN));
-// };
-// LSB(EXP_MID4[i]) = 2^-127
-LIBC_INLINE_VAR constexpr Frac128 EXP_MID[] = {
- Frac128({0x0000'0000'0000'0000ULL, 0x8000'0000'0000'0000ULL}),
- Frac128({0xc5c9'5b8c'2154'c1b2ULL, 0x85aa'c367'cc48'7b14ULL}),
- Frac128({0xfbe4'6287'58a5'3c90ULL, 0x8b95'c1e3'ea8b'd6e6ULL}),
- Frac128({0x0fd6'd8e0'ae5a'c9d8ULL, 0x91c3'd373'ab11'c336ULL}),
- Frac128({0x46ad'2318'2e42'f6f6ULL, 0x9837'f051'8db8'a96fULL}),
- Frac128({0xa091'1f09'ebb9'fdd1ULL, 0x9ef5'3260'91a1'11adULL}),
- Frac128({0x1cbd'7f62'1710'701bULL, 0xa5fe'd6a9'b151'38eaULL}),
- Frac128({0x4980'a8c8'f59a'2ec4ULL, 0xad58'3eea'42a1'4ac6ULL}),
- Frac128({0x597d'89b3'754a'be9fULL, 0xb504'f333'f9de'6484ULL}),
- Frac128({0xa881'1fb6'6d0f'af7aULL, 0xbd08'a39f'580c'36beULL}),
- Frac128({0x3e2a'd0c9'64dd'9f37ULL, 0xc567'2a11'5506'daddULL}),
- Frac128({0xe235'838f'95f2'c6edULL, 0xce24'8c15'1f84'80e3ULL}),
- Frac128({0x39a6'8bb9'902d'3fdeULL, 0xd744'fcca'd69d'6af4ULL}),
- Frac128({0x0658'9504'8dd3'33caULL, 0xe0cc'deec'2a94'e111ULL}),
- Frac128({0xd02d'75b3'706e'54fbULL, 0xeac0'c6e7'dd24'392eULL}),
- Frac128({0x7b9d'0c7a'ed98'0fc3ULL, 0xf525'7d15'2486'cc2cULL}),
-};
-
-// 128-bit polynomial approximation of 2^x coefficients generated with Sollya:
-// > P = fpminimax(2^x, 12, [|1, 128...|], [0, 1/16], absolute, fixed);
-// Store the fractional part of the coefficients below
-// > dirtyinfnorm(2^x - P(x), [0, 1/16]);
-// 0x1.b328...p-117
-// LSB(EXPF_COEFFS[i]) = 2^-128
-LIBC_INLINE_VAR constexpr Frac128 EXP_COEFFS[] = {
- // degree-0 = 1, add back afterwards to reduce calc ops
- Frac128({0xc9e3'b398'033f'0902ULL, 0xb172'17f7'd1cf'79abULL}),
- Frac128({0xde2d'60e0'866c'6365ULL, 0x3d7f'7bff'058b'1d50ULL}),
- Frac128({0x99d3'ad04'd47d'0efeULL, 0x0e35'846b'8250'5fc5ULL}),
- Frac128({0x399a'b423'1644'1b55ULL, 0x0276'556d'f749'cee5ULL}),
- Frac128({0x405f'ebef'0f76'011aULL, 0x0057'61ff'9e29'9cc4ULL}),
- Frac128({0x1ac4'9999'c0ef'b9abULL, 0x000a'1848'97c3'63c4ULL}),
- Frac128({0xe6f8'52f9'914f'2d9aULL, 0x0000'ffe5'fe2c'4573ULL}),
- Frac128({0x6056'28a6'5061'b43fULL, 0x0000'162c'0223'a816ULL}),
- Frac128({0x6b99'3efe'2193'e54cULL, 0x0000'01b5'253d'0671ULL}),
- Frac128({0x9c76'f49c'ed65'743aULL, 0x0000'001e'4cf8'0b70ULL}),
- Frac128({0x1ac1'3330'6c4e'3d33ULL, 0x0000'0001'e8ae'6938ULL}),
- Frac128({0x219f'9904'1f29'5f13ULL, 0x0000'0000'1cd9'af73ULL}),
-};
-
-// 64-bit polynomial approximation of 2^x coefficients generated with Sollya:
-// > P = fpminimax(2^x, 6, [|1, 64...|], [0, 1/16], absolute, fixed);
-// Store the fractional part of the coefficients below
-// > dirtyinfnorm(2^x - P(x), [0, 1/16]);
-// 0x1.2ac9...p-57
-// This is different from EXPF_COEFFS: EXPF_COEFFS is for approximating for x in
-// range of [0, 1]
-// LSB(EXP_COEFFS[i]) = 2^-64
-LIBC_INLINE_VAR constexpr Frac64 EXP_64_COEFFS[] = {
- // degree-0 = 1, add back afterwards to reduce calc ops
- Frac64(0xb172'17f7'd1cd'3e3cULL), Frac64(0x3d7f'7bff'0838'4d4aULL),
- Frac64(0x0e35'846a'6f46'd462ULL), Frac64(0x0276'55a0'd536'1dd9ULL),
- Frac64(0x0057'5d3d'4b45'0056ULL), Frac64(0x000a'504b'13fe'e008ULL),
-};
-
// Round the fractional result and combine it with its exponent.
template <typename TFrac, typename TUInt,
cpp::enable_if_t<cpp::is_same<TFrac, Frac64>::value ||
@@ -244,14 +176,11 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
uint64_t x_val = xbits.uintval();
uint64_t x_val_abs = xbits.abs().uintval();
- // TODO: dedupe? Ported over from the base LLVM-libc exp(x) implementation
- // with same checks, just without exceptions/errnos.
- //
// x < log(2^-1075) or x >= 0x1.6232bdd7abcd3p+9 or |x| < 2^-53.
- if (LIBC_UNLIKELY(
- x_val >= 0xc0874910d52d3052 ||
- (x_val < 0xbca0000000000000 && x_val >= 0x40862e42fefa39f0) ||
- x_val < 0x3ca0000000000000)) {
+ if (LIBC_UNLIKELY(x_val >= 0xc087'4910'd52d'3052ULL ||
+ (x_val < 0xbca0'0000'0000'0000ULL &&
+ x_val >= 0x4086'2e42'fefa'39f0ULL) ||
+ x_val < 0x3ca0'0000'0000'0000ULL)) {
// |x| <= 2^-53
if (x_val_abs <= 0x3ca0'0000'0000'0000ULL) {
// exp(x) ~ 1 + x
@@ -321,8 +250,7 @@ LIBC_INLINE double exp(double x, [[maybe_unused]] int rounding) {
// LSB(x_s_frac) = 2^-52
Frac64 x_s_frac(x_s_shifted);
- constexpr Frac64 INV_LN2_HI = INV_LN2_F128.to_frac64();
- Frac64 x_ln2 = x_s_frac * INV_LN2_HI;
+ Frac64 x_ln2 = x_s_frac * INV_LN2_F64;
uint64_t x_ln2_bits = x_ln2.val[0];
uint64_t k = 0;
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 52736323d8d7e..096bb1facf7e4 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -15,6 +15,7 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_INTEGER_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPF_INTEGER_EVAL_H
+#include "exp_integer_constants.h" // LUTs
#include "hdr/fenv_macros.h"
#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FPBits.h"
@@ -32,30 +33,6 @@ 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 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, [[maybe_unused]] int rounding) {
@@ -159,7 +136,7 @@ LIBC_INLINE float expf(float x, [[maybe_unused]] int rounding) {
Frac64 x_u_frac(x_u);
// LSB(x_ln2) = 2^-54
- Frac64 x_ln2 = x_u_frac * INV_LN2;
+ Frac64 x_ln2 = x_u_frac * INV_LN2_F64;
constexpr uint64_t FRAC_MASK = (uint64_t(1) << 54) - 1;
uint64_t x_ln2_bit = x_ln2.val[0];
@@ -174,7 +151,7 @@ LIBC_INLINE float expf(float x, [[maybe_unused]] int rounding) {
// - For positive x, we round down (just clear) the last 54 bits.
//
// Then, l2y_r is the remainder of x * log2(e) after removing the integer
- // part, which is used to compute 2^l2y_r_frac - 1.
+ // part, which is used to compute 2^l2y_r_frac.
//
// e_y_unbiased is biased exponent field, but already bit-positioned to the
// exponent field of the float representation.
@@ -202,7 +179,7 @@ LIBC_INLINE float expf(float x, [[maybe_unused]] int rounding) {
// LSB(l2y_r_frac) = LSB(l2y_r) * 2^-10 = 2^-64
Frac64 l2y_r_frac(l2y_r << 10);
- // p = 2^l2y_r_frac - 1
+ // p = 2^l2y_r_frac
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],
@@ -212,7 +189,7 @@ LIBC_INLINE float expf(float x, [[maybe_unused]] int rounding) {
uint32_t shift_length = 40;
uint32_t leading_one = 0;
- // We're computing with errors < worst-cast errors, so tie-rounding never
+ // We're computing with errors < worst-case errors, so tie-rounding never
// happens. Hence, round-to-nearest, tie-to-even is equivalent to
// round-to-nearest, tie-to-away. Which is what we're implementing below
// in the following order:
More information about the libc-commits
mailing list