[libc-commits] [libc] initial commit (PR #225803)
Hoàng Minh Thiên via libc-commits
libc-commits at lists.llvm.org
Sat Sep 26 08:26:55 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/2] 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/2] 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
More information about the libc-commits
mailing list