[libc-commits] [libc] [llvm] [libc][math] Optimize exp2 and refactor shared double-precision exp utilities. (PR #231041)
via libc-commits
libc-commits at lists.llvm.org
Sun Oct 11 15:14:39 PDT 2026
https://github.com/lntue created https://github.com/llvm/llvm-project/pull/231041
Refactor double-precision `exp2` to improve latency, throughput, and code
sharing across `exp`, `exp2`, `exp10`, and `pow`, while adding support for
`LIBC_MATH_HAS_SMALL_TABLES` and `LIBC_MATH_HAS_SKIP_ACCURATE_PASS`.
Algorithm & implementation overview:
- Fast path (`exp2_double_eval.h`, `exp_utils.h`):
- Single-step range reduction `x = k * 2^-N + dx` (`N = 7` by default, or
`N = 4` when `LIBC_MATH_HAS_SMALL_TABLES` is enabled) with
`hi = k >> N`, `idx = k & (2^N - 1)`, and `u = dx * log(2)`.
- Without FMA, use fixed-point splitting (`SPLIT_18` on `dx`, `LOG_2_HI_12`
on `log(2)`, and `SPLIT_22` on `EXP2_MID[idx]`) so `m_hi * (1 + u_hi)` is
exact in `double` without needing Veltkamp splitting or `DoubleDouble`
multiplications.
- Approximate `(exp(u) - 1 - u) / u^2` with a degree-3 minimax polynomial
evaluated with Estrin's scheme (`exp_poly`), and reconstruct
`2^(idx / 2^N) * exp(u)` followed by Ziv's rounding test.
- Scale `2^hi * r` directly in vector registers (`scale_exp2`) using SSE2
and NEON integer additions on the exponent field to avoid
float-to-GPR-to-float round-trips.
- Single-precision / skip-accurate-pass variant (`exp_double_utils.h`):
- When `LIBC_MATH_HAS_SKIP_ACCURATE_PASS` is enabled, evaluate `exp2_eval`
purely in `double` precision (~0.51 ULP error, <= 1 ULP in RN and <= 2 ULPs
in directed rounding modes).
- 128-bit accurate path (`exp2_accurate_128.h`):
- Replace the previous 2-step `DyadicFloat<128>` lookup tables (`2^-6` and
`2^-12`) and degree-6 `DyadicFloat<128>` Taylor polynomial with a single
16-entry `Frac128` lookup table (`EXP2_MID_FRAC128`, `2^(i/16)`) and a
degree-14 minimax polynomial (`EXP2_POLY_128`) in 128-bit fixed-point
arithmetic.
- Provide `exp_accurate_128` and `exp2_df128` to be shared across `exp2`,
`exp`, `exp10`, and `pow`.
- File organization & tests:
- Split `double_eval::exp2` into `libc/src/__support/math/exp2_double_eval.h`
so `double_eval::exp2` can be tested directly alongside `libc::exp2`.
- Add a `LIBC_NOINLINE` wrapper around `double_eval::exp2` in unit and smoke
tests to prevent excessive inlining into test functions.
Assisted-by: Gemini was used for performance analysis, error bound verification, and refactoring.
>From bb0b0f1a728295ac6a00db6012eb232f6875ea5b Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Sun, 11 Oct 2026 22:05:07 +0000
Subject: [PATCH] [libc][math] Optimize exp2 and refactor shared
double-precision exp utilities.
---
libc/src/__support/FPUtil/PolyEval.h | 13 +-
libc/src/__support/FPUtil/multiply_add.h | 13 +-
libc/src/__support/math/CMakeLists.txt | 75 ++-
libc/src/__support/math/exp2.h | 437 +---------------
libc/src/__support/math/exp2_accurate_128.h | 228 +++++++++
libc/src/__support/math/exp2_double_eval.h | 133 +++++
libc/src/__support/math/exp_double_utils.h | 175 +++++++
libc/src/__support/math/exp_utils.h | 475 +++++++++++++++++-
libc/test/src/math/CMakeLists.txt | 2 +
libc/test/src/math/exp2_test.cpp | 255 ++++++----
libc/test/src/math/smoke/CMakeLists.txt | 2 +
libc/test/src/math/smoke/exp2_test.cpp | 122 +++--
.../llvm-project-overlay/libc/BUILD.bazel | 72 ++-
.../libc/test/src/math/BUILD.bazel | 7 +-
.../libc/test/src/math/smoke/BUILD.bazel | 7 +-
15 files changed, 1425 insertions(+), 591 deletions(-)
create mode 100644 libc/src/__support/math/exp2_accurate_128.h
create mode 100644 libc/src/__support/math/exp2_double_eval.h
create mode 100644 libc/src/__support/math/exp_double_utils.h
diff --git a/libc/src/__support/FPUtil/PolyEval.h b/libc/src/__support/FPUtil/PolyEval.h
index 58a3f4e24fdeb..55583cca70c3b 100644
--- a/libc/src/__support/FPUtil/PolyEval.h
+++ b/libc/src/__support/FPUtil/PolyEval.h
@@ -1,10 +1,15 @@
-//===-- Common header for PolyEval implementations --------------*- C++ -*-===//
+//===----------------------------------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Common header for PolyEval implementations.
+///
+//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_POLYEVAL_H
#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_POLYEVAL_H
@@ -67,13 +72,15 @@ altpolyeval(T, T a0) {
template <typename T, typename... Ts>
LIBC_INLINE constexpr cpp::enable_if_t<(sizeof(T) > sizeof(void *)), T>
altpolyeval(const T &x, const T &a0, const Ts &...a) {
- return a0 - x * altpolyeval(x, a...);
+ T prod = x * altpolyeval(x, a...);
+ return a0 - prod;
}
template <typename T, typename... Ts>
LIBC_INLINE constexpr cpp::enable_if_t<(sizeof(T) <= sizeof(void *)), T>
altpolyeval(T x, T a0, Ts... a) {
- return a0 - x * altpolyeval(x, a...);
+ T prod = x * altpolyeval(x, a...);
+ return a0 - prod;
}
} // namespace fputil
diff --git a/libc/src/__support/FPUtil/multiply_add.h b/libc/src/__support/FPUtil/multiply_add.h
index dce09622cd668..fe87877d2443d 100644
--- a/libc/src/__support/FPUtil/multiply_add.h
+++ b/libc/src/__support/FPUtil/multiply_add.h
@@ -1,10 +1,15 @@
-//===-- Common header for multiply-add implementations ----------*- C++ -*-===//
+//===----------------------------------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Common header for multiply-add implementations.
+///
+//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_FPUTIL_MULTIPLY_ADD_H
#define LLVM_LIBC_SRC___SUPPORT_FPUTIL_MULTIPLY_ADD_H
@@ -25,13 +30,15 @@ namespace fputil {
template <typename T>
LIBC_INLINE constexpr cpp::enable_if_t<(sizeof(T) > sizeof(void *)), T>
multiply_add(const T &x, const T &y, const T &z) {
- return x * y + z;
+ T prod = x * y;
+ return prod + z;
}
template <typename T>
LIBC_INLINE constexpr cpp::enable_if_t<(sizeof(T) <= sizeof(void *)), T>
multiply_add(T x, T y, T z) {
- return x * y + z;
+ T prod = x * y;
+ return prod + z;
}
} // namespace fputil
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index ef33e3f8ada7b..3a19ee880fe06 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -4195,14 +4195,43 @@ add_header_library(
libc.src.__support.FPUtil.triple_double
)
+add_header_library(
+ exp2_accurate_128
+ HDRS
+ exp2_accurate_128.h
+ DEPENDS
+ libc.src.__support.FPUtil.dyadic_float
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.common
+ libc.src.__support.frac128
+ libc.src.__support.frac64
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.uint128
+)
+
add_header_library(
exp_utils
HDRS
exp_utils.h
DEPENDS
- libc.src.__support.CPP.optional
+ .exp2_accurate_128
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
libc.src.__support.CPP.bit
+ libc.src.__support.CPP.optional
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.rounding_mode
+ libc.src.__support.common
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.architectures
+ libc.src.__support.macros.properties.cpu_features
)
add_header_library(
@@ -4227,25 +4256,49 @@ add_header_library(
)
add_header_library(
- exp2
+ exp_double_utils
HDRS
- exp2.h
+ exp_double_utils.h
DEPENDS
- .common_constants
.exp_utils
+ libc.hdr.errno_macros
libc.src.__support.CPP.bit
- libc.src.__support.CPP.optional
- libc.src.__support.FPUtil.dyadic_float
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.common
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+)
+
+add_header_library(
+ exp2_double_eval
+ HDRS
+ exp2_double_eval.h
+ DEPENDS
+ .exp_double_utils
+ .exp_utils
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.nearest_integer
- libc.src.__support.FPUtil.polyeval
- libc.src.__support.FPUtil.rounding_mode
- libc.src.__support.FPUtil.triple_double
- libc.src.__support.integer_literals
+ libc.src.__support.common
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
+)
+
+add_header_library(
+ exp2
+ HDRS
+ exp2.h
+ DEPENDS
+ .exp2_double_eval
+ libc.src.__support.common
+ libc.src.__support.macros.config
libc.src.__support.macros.optimization
- libc.src.errno.errno
)
add_header_library(
diff --git a/libc/src/__support/math/exp2.h b/libc/src/__support/math/exp2.h
index 084578796e422..12bad5db2a9b2 100644
--- a/libc/src/__support/math/exp2.h
+++ b/libc/src/__support/math/exp2.h
@@ -1,445 +1,34 @@
-//===-- Implementation header for exp2 --------------------------*- C++ -*-===//
+//===----------------------------------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Implementation header for exp2.
+///
+//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_H
-#include "common_constants.h" // Lookup tables EXP2_MID1 and EXP_M2.
-#include "exp_constants.h"
-#include "exp_utils.h" // ziv_test_denorm.
-#include "src/__support/CPP/bit.h"
-#include "src/__support/CPP/optional.h"
-#include "src/__support/FPUtil/FEnvImpl.h"
-#include "src/__support/FPUtil/FPBits.h"
-#include "src/__support/FPUtil/PolyEval.h"
-#include "src/__support/FPUtil/double_double.h"
-#include "src/__support/FPUtil/dyadic_float.h"
-#include "src/__support/FPUtil/multiply_add.h"
-#include "src/__support/FPUtil/nearest_integer.h"
-#include "src/__support/FPUtil/rounding_mode.h"
-#include "src/__support/FPUtil/triple_double.h"
#include "src/__support/common.h"
-#include "src/__support/integer_literals.h"
#include "src/__support/macros/config.h"
-#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/macros/optimization.h"
-namespace LIBC_NAMESPACE_DECL {
+#include "src/__support/math/exp2_double_eval.h"
+#define LIBC_MATH_EXP2_IMPL double_eval
+namespace LIBC_NAMESPACE_DECL {
namespace math {
-namespace exp2_internal {
-
-using namespace common_constants_internal;
-
-using fputil::DoubleDouble;
-using fputil::TripleDouble;
-using DFloat128 = typename fputil::DyadicFloat<128>;
-
-using LIBC_NAMESPACE::operator""_u128;
-
-// Error bounds:
-// Errors when using double precision.
-#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
-LIBC_INLINE_VAR constexpr double ERR_D = 0x1.0p-63;
-#else
-LIBC_INLINE_VAR constexpr double ERR_D = 0x1.8p-63;
-#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-// Errors when using double-double precision.
-LIBC_INLINE_VAR constexpr double ERR_DD = 0x1.0p-100;
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
-// Polynomial approximations with double precision. Generated by Sollya with:
-// > P = fpminimax((2^x - 1)/x, 3, [|D...|], [-2^-13 - 2^-30, 2^-13 + 2^-30]);
-// > P;
-// Error bounds:
-// | output - (2^dx - 1) / dx | < 1.5 * 2^-52.
-LIBC_INLINE double poly_approx_d(double dx) {
- // dx^2
- double dx2 = dx * dx;
- double c0 =
- fputil::multiply_add(dx, 0x1.ebfbdff82c58ep-3, 0x1.62e42fefa39efp-1);
- double c1 =
- fputil::multiply_add(dx, 0x1.3b2aba7a95a89p-7, 0x1.c6b08e8fc0c0ep-5);
- double p = fputil::multiply_add(dx2, c1, c0);
- return p;
-}
-
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-// Polynomial approximation with double-double precision. Generated by Solya
-// with:
-// > P = fpminimax((2^x - 1)/x, 5, [|DD...|], [-2^-13 - 2^-30, 2^-13 + 2^-30]);
-// Error bounds:
-// | output - 2^(dx) | < 2^-101
-LIBC_INLINE DoubleDouble poly_approx_dd(const DoubleDouble &dx) {
- // Taylor polynomial.
- constexpr DoubleDouble COEFFS[] = {
- {0, 0x1p0},
- {0x1.abc9e3b39824p-56, 0x1.62e42fefa39efp-1},
- {-0x1.5e43a53e4527bp-57, 0x1.ebfbdff82c58fp-3},
- {-0x1.d37963a9444eep-59, 0x1.c6b08d704a0cp-5},
- {0x1.4eda1a81133dap-62, 0x1.3b2ab6fba4e77p-7},
- {-0x1.c53fd1ba85d14p-64, 0x1.5d87fe7a265a5p-10},
- {0x1.d89250b013eb8p-70, 0x1.430912f86cb8ep-13},
- };
-
- DoubleDouble p = fputil::polyeval(dx, COEFFS[0], COEFFS[1], COEFFS[2],
- COEFFS[3], COEFFS[4], COEFFS[5], COEFFS[6]);
- return p;
-}
-
-// Polynomial approximation with 128-bit precision:
-// Return exp(dx) ~ 1 + a0 * dx + a1 * dx^2 + ... + a6 * dx^7
-// For |dx| < 2^-13 + 2^-30:
-// | output - exp(dx) | < 2^-126.
-LIBC_INLINE DFloat128 poly_approx_f128(const DFloat128 &dx) {
- constexpr DFloat128 COEFFS_128[]{
- {Sign::POS, -127, 0x80000000'00000000'00000000'00000000_u128}, // 1.0
- {Sign::POS, -128, 0xb17217f7'd1cf79ab'c9e3b398'03f2f6af_u128},
- {Sign::POS, -128, 0x3d7f7bff'058b1d50'de2d60dd'9c9a1d9f_u128},
- {Sign::POS, -132, 0xe35846b8'2505fc59'9d3b15d9'e7fb6897_u128},
- {Sign::POS, -134, 0x9d955b7d'd273b94e'184462f6'bcd2b9e7_u128},
- {Sign::POS, -137, 0xaec3ff3c'53398883'39ea1bb9'64c51a89_u128},
- {Sign::POS, -138, 0x2861225f'345c396a'842c5341'8fa8ae61_u128},
- {Sign::POS, -144, 0xffe5fe2d'109a319d'7abeb5ab'd5ad2079_u128},
- };
-
- DFloat128 p = fputil::polyeval(dx, COEFFS_128[0], COEFFS_128[1],
- COEFFS_128[2], COEFFS_128[3], COEFFS_128[4],
- COEFFS_128[5], COEFFS_128[6], COEFFS_128[7]);
- return p;
-}
-
-// Compute 2^(x) using 128-bit precision.
-// TODO(lntue): investigate triple-double precision implementation for this
-// step.
-LIBC_INLINE DFloat128 exp2_f128(double x, int hi, int idx1, int idx2) {
- DFloat128 dx = DFloat128(x);
-
- // TODO: Skip recalculating exp_mid1 and exp_mid2.
- DFloat128 exp_mid1 =
- fputil::quick_add(DFloat128(EXP2_MID1[idx1].hi),
- fputil::quick_add(DFloat128(EXP2_MID1[idx1].mid),
- DFloat128(EXP2_MID1[idx1].lo)));
-
- DFloat128 exp_mid2 =
- fputil::quick_add(DFloat128(EXP2_MID2[idx2].hi),
- fputil::quick_add(DFloat128(EXP2_MID2[idx2].mid),
- DFloat128(EXP2_MID2[idx2].lo)));
-
- DFloat128 exp_mid = fputil::quick_mul(exp_mid1, exp_mid2);
-
- DFloat128 p = poly_approx_f128(dx);
-
- DFloat128 r = fputil::quick_mul(exp_mid, p);
-
- r.exponent += hi;
-
- return r;
-}
-
-// Compute 2^x with double-double precision.
-LIBC_INLINE DoubleDouble exp2_double_double(double x,
- const DoubleDouble &exp_mid) {
- DoubleDouble dx({0, x});
-
- // Degree-6 polynomial approximation in double-double precision.
- // | p - 2^x | < 2^-103.
- DoubleDouble p = poly_approx_dd(dx);
-
- // Error bounds: 2^-102.
- DoubleDouble r = fputil::quick_mult(exp_mid, p);
-
- return r;
-}
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-
-// When output is denormal.
-LIBC_INLINE double exp2_denorm(double x) {
- // Range reduction.
- int k =
- static_cast<int>(cpp::bit_cast<uint64_t>(x + 0x1.8000'0000'4p21) >> 19);
- double kd = static_cast<double>(k);
-
- uint32_t idx1 = (k >> 6) & 0x3f;
- uint32_t idx2 = k & 0x3f;
-
- int hi = k >> 12;
-
- DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi};
- DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi};
- DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2);
-
- // |dx| < 2^-13 + 2^-30.
- double dx = fputil::multiply_add(kd, -0x1.0p-12, x); // exact
-
- double mid_lo = dx * exp_mid.hi;
-
- // Approximate (2^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4.
- double p = poly_approx_d(dx);
-
- double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
-
-#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- return ziv_test_denorm</*SKIP_ZIV_TEST=*/true>(hi, exp_mid.hi, lo, ERR_D)
- .value();
-#else
- using FPBits = fputil::FPBits<double>;
-
- if (auto r = ziv_test_denorm(hi, exp_mid.hi, lo, ERR_D);
- LIBC_LIKELY(r.has_value())) {
- double res = r.value();
- if (LIBC_UNLIKELY(FPBits(res).is_normal()))
- return res;
- if (x != fputil::nearest_integer(x)) {
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_underflow_except_if_required<double>();
- }
- return res;
- }
-
- // Use double-double
- DoubleDouble r_dd = exp2_double_double(dx, exp_mid);
-
- if (auto r = ziv_test_denorm(hi, r_dd.hi, r_dd.lo, ERR_DD);
- LIBC_LIKELY(r.has_value())) {
- double res = r.value();
- if (LIBC_UNLIKELY(FPBits(res).is_normal()))
- return res;
- if (x != fputil::nearest_integer(x)) {
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_underflow_except_if_required<double>();
- }
- return res;
- }
-
- // Use 128-bit precision
- DFloat128 r_f128 = exp2_f128(dx, hi, idx1, idx2);
-
- return static_cast<double>(r_f128);
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-}
-
-// Check for exceptional cases when:
-// * log2(1 - 2^-54) < x < log2(1 + 2^-53)
-// * x >= 1024
-// * x <= -1022
-// * x is inf or nan
-LIBC_INLINE double set_exceptional(double x) {
- using FPBits = typename fputil::FPBits<double>;
- FPBits xbits(x);
-
- uint64_t x_u = xbits.uintval();
- uint64_t x_abs = xbits.abs().uintval();
-
- // |x| < log2(1 + 2^-53)
- if (x_abs <= 0x3ca71547652b82fd) {
- // 2^(x) ~ 1 + x/2
- return fputil::multiply_add(x, 0.5, 1.0);
- }
-
- // x <= -1022 || x >= 1024 or inf/nan.
- if (x_u > 0xc08ff00000000000) {
- // x <= -1075 or -inf/nan
- if (x_u >= 0xc090cc0000000000) {
- // exp(-Inf) = 0
- if (xbits.is_inf())
- return 0.0;
-
- // exp(nan) = nan
- if (xbits.is_nan())
- return x + FPBits::inf().get_val();
-
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_underflow_except_if_required<double>();
-#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- if (fputil::quick_get_round() == FE_UPWARD)
- return FPBits::min_subnormal().get_val();
-#endif
- return 0.0;
- }
-
- return exp2_denorm(x);
- }
-
- // x >= 1024 or +inf/nan
- // x is finite
- if (x_u < 0x7ff0'0000'0000'0000ULL) {
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_overflow_except_if_required<double>();
-#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
- int rounding = fputil::quick_get_round();
- if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
- return FPBits::max_normal().get_val();
-#endif
- return FPBits::inf().get_val();
- }
- // x is +inf or nan
- return x + FPBits::inf().get_val();
-}
-
-} // namespace exp2_internal
-
-LIBC_INLINE double exp2(double x) {
- using namespace exp2_internal;
- using FPBits = typename fputil::FPBits<double>;
- FPBits xbits(x);
-
- uint64_t x_u = xbits.uintval();
-
- // x < -1022 or x >= 1024 or log2(1 - 2^-54) < x < log2(1 + 2^-53).
- if (LIBC_UNLIKELY(x_u > 0xc08ff00000000000 ||
- (x_u <= 0xbc971547652b82fe && x_u >= 0x4090000000000000) ||
- x_u <= 0x3ca71547652b82fd)) {
- return set_exceptional(x);
- }
-
- // Now -1075 < x <= log2(1 - 2^-54) or log2(1 + 2^-53) < x < 1024
-
- // Range reduction:
- // Let x = (hi + mid1 + mid2) + lo
- // in which:
- // hi is an integer
- // mid1 * 2^6 is an integer
- // mid2 * 2^12 is an integer
- // then:
- // 2^(x) = 2^hi * 2^(mid1) * 2^(mid2) * 2^(lo).
- // With this formula:
- // - multiplying by 2^hi is exact and cheap, simply by adding the exponent
- // field.
- // - 2^(mid1) and 2^(mid2) are stored in 2 x 64-element tables.
- // - 2^(lo) ~ 1 + a0*lo + a1 * lo^2 + ...
- //
- // We compute (hi + mid1 + mid2) together by perform the rounding on x * 2^12.
- // Since |x| < |-1075)| < 2^11,
- // |x * 2^12| < 2^11 * 2^12 < 2^23,
- // So we can fit the rounded result round(x * 2^12) in int32_t.
- // Thus, the goal is to be able to use an additional addition and fixed width
- // shift to get an int32_t representing round(x * 2^12).
- //
- // Assuming int32_t using 2-complement representation, since the mantissa part
- // of a double precision is unsigned with the leading bit hidden, if we add an
- // extra constant C = 2^e1 + 2^e2 with e1 > e2 >= 2^25 to the product, the
- // part that are < 2^e2 in resulted mantissa of (x*2^12*L2E + C) can be
- // considered as a proper 2-complement representations of x*2^12.
- //
- // One small problem with this approach is that the sum (x*2^12 + C) in
- // double precision is rounded to the least significant bit of the dorminant
- // factor C. In order to minimize the rounding errors from this addition, we
- // want to minimize e1. Another constraint that we want is that after
- // shifting the mantissa so that the least significant bit of int32_t
- // corresponds to the unit bit of (x*2^12*L2E), the sign is correct without
- // any adjustment. So combining these 2 requirements, we can choose
- // C = 2^33 + 2^32, so that the sign bit corresponds to 2^31 bit, and hence
- // after right shifting the mantissa, the resulting int32_t has correct sign.
- // With this choice of C, the number of mantissa bits we need to shift to the
- // right is: 52 - 33 = 19.
- //
- // Moreover, since the integer right shifts are equivalent to rounding down,
- // we can add an extra 0.5 so that it will become round-to-nearest, tie-to-
- // +infinity. So in particular, we can compute:
- // hmm = x * 2^12 + C,
- // where C = 2^33 + 2^32 + 2^-1, then if
- // k = int32_t(lower 51 bits of double(x * 2^12 + C) >> 19),
- // the reduced argument:
- // lo = x - 2^-12 * k is bounded by:
- // |lo| <= 2^-13 + 2^-12*2^-19
- // = 2^-13 + 2^-31.
- //
- // Finally, notice that k only uses the mantissa of x * 2^12, so the
- // exponent 2^12 is not needed. So we can simply define
- // C = 2^(33 - 12) + 2^(32 - 12) + 2^(-13 - 12), and
- // k = int32_t(lower 51 bits of double(x + C) >> 19).
-
- // Rounding errors <= 2^-31.
- int k =
- static_cast<int>(cpp::bit_cast<uint64_t>(x + 0x1.8000'0000'4p21) >> 19);
- double kd = static_cast<double>(k);
-
- uint32_t idx1 = (k >> 6) & 0x3f;
- uint32_t idx2 = k & 0x3f;
-
- int hi = k >> 12;
-
- DoubleDouble exp_mid1{EXP2_MID1[idx1].mid, EXP2_MID1[idx1].hi};
- DoubleDouble exp_mid2{EXP2_MID2[idx2].mid, EXP2_MID2[idx2].hi};
- DoubleDouble exp_mid = fputil::quick_mult(exp_mid1, exp_mid2);
-
- // |dx| < 2^-13 + 2^-30.
- double dx = fputil::multiply_add(kd, -0x1.0p-12, x); // exact
-
- // We use the degree-4 polynomial to approximate 2^(lo):
- // 2^(lo) ~ 1 + a0 * lo + a1 * lo^2 + a2 * lo^3 + a3 * lo^4 = 1 + lo * P(lo)
- // So that the errors are bounded by:
- // |P(lo) - (2^lo - 1)/lo| < |lo|^4 / 64 < 2^(-13 * 4) / 64 = 2^-58
- // Let P_ be an evaluation of P where all intermediate computations are in
- // double precision. Using either Horner's or Estrin's schemes, the evaluated
- // errors can be bounded by:
- // |P_(lo) - P(lo)| < 2^-51
- // => |lo * P_(lo) - (2^lo - 1) | < 2^-64
- // => 2^(mid1 + mid2) * |lo * P_(lo) - expm1(lo)| < 2^-63.
- // Since we approximate
- // 2^(mid1 + mid2) ~ exp_mid.hi + exp_mid.lo,
- // We use the expression:
- // (exp_mid.hi + exp_mid.lo) * (1 + dx * P_(dx)) ~
- // ~ exp_mid.hi + (exp_mid.hi * dx * P_(dx) + exp_mid.lo)
- // with errors bounded by 2^-63.
-
- double mid_lo = dx * exp_mid.hi;
-
- // Approximate (2^dx - 1)/dx ~ 1 + a0*dx + a1*dx^2 + a2*dx^3 + a3*dx^4.
- double p = poly_approx_d(dx);
-
- double lo = fputil::multiply_add(p, mid_lo, exp_mid.lo);
-
-#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- // To multiply by 2^hi, a fast way is to simply add hi to the exponent
- // field.
- int64_t exp_hi = static_cast<int64_t>(hi) * (1LL << FPBits::FRACTION_LEN);
- double r =
- cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(exp_mid.hi + lo));
- return r;
-#else
- double upper = exp_mid.hi + (lo + ERR_D);
- double lower = exp_mid.hi + (lo - ERR_D);
-
- if (LIBC_LIKELY(upper == lower)) {
- // To multiply by 2^hi, a fast way is to simply add hi to the exponent
- // field.
- int64_t exp_hi = static_cast<int64_t>(hi) * (1LL << FPBits::FRACTION_LEN);
- double r = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(upper));
- return r;
- }
-
- // Use double-double
- DoubleDouble r_dd = exp2_double_double(dx, exp_mid);
-
- double upper_dd = r_dd.hi + (r_dd.lo + ERR_DD);
- double lower_dd = r_dd.hi + (r_dd.lo - ERR_DD);
-
- if (LIBC_LIKELY(upper_dd == lower_dd)) {
- // To multiply by 2^hi, a fast way is to simply add hi to the exponent
- // field.
- int64_t exp_hi = static_cast<int64_t>(hi) * (1LL << FPBits::FRACTION_LEN);
- double r = cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(upper_dd));
- return r;
- }
-
- // Use 128-bit precision
- DFloat128 r_f128 = exp2_f128(dx, hi, idx1, idx2);
-
- return static_cast<double>(r_f128);
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-}
+using LIBC_MATH_EXP2_IMPL::exp2;
} // namespace math
-
} // namespace LIBC_NAMESPACE_DECL
+#undef LIBC_MATH_EXP2_IMPL
+
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_H
diff --git a/libc/src/__support/math/exp2_accurate_128.h b/libc/src/__support/math/exp2_accurate_128.h
new file mode 100644
index 0000000000000..07fc870a5f3e7
--- /dev/null
+++ b/libc/src/__support/math/exp2_accurate_128.h
@@ -0,0 +1,228 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// 128-bit accurate path for double-precision exp2(x).
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_ACCURATE_128_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_ACCURATE_128_H
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/dyadic_float.h"
+#include "src/__support/common.h"
+#include "src/__support/frac128.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/uint128.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace double_eval {
+
+using DFloat128 = typename fputil::DyadicFloat<128>;
+
+// Table of 2^(i / 16) with i = 0..15 in 128-bit fixed-point format:
+// Generated by Sollya with:
+// > display = hexadecimal!;
+// > prec = 256;
+// > for i from 0 to 15 do {
+// a = 2^(i / 16);
+// c_int = round(a * 2^127, 256, RN);
+// hi = floor(c_int * 2^(-64));
+// lo = c_int - hi * 2^64;
+// print(" Frac128({", lo, "ULL,", hi, "ULL}), // 2^(", i, "/16)");
+// };
+LIBC_INLINE_VAR constexpr Frac128 EXP2_MID_FRAC128[16] = {
+ // lsb = 2^-127
+ Frac128({0x0000000000000000ULL, 0x8000000000000000ULL}), // 2^(0/16)
+ Frac128({0xc5c95b8c2154c1b2ULL, 0x85aac367cc487b14ULL}), // 2^(1/16)
+ Frac128({0xfbe4628758a53c90ULL, 0x8b95c1e3ea8bd6e6ULL}), // 2^(2/16)
+ Frac128({0x0fd6d8e0ae5ac9d8ULL, 0x91c3d373ab11c336ULL}), // 2^(3/16)
+ Frac128({0x46ad23182e42f6f6ULL, 0x9837f0518db8a96fULL}), // 2^(4/16)
+ Frac128({0xa0911f09ebb9fdd1ULL, 0x9ef5326091a111adULL}), // 2^(5/16)
+ Frac128({0x1cbd7f621710701bULL, 0xa5fed6a9b15138eaULL}), // 2^(6/16)
+ Frac128({0x4980a8c8f59a2ec4ULL, 0xad583eea42a14ac6ULL}), // 2^(7/16)
+ Frac128({0x597d89b3754abe9fULL, 0xb504f333f9de6484ULL}), // 2^(8/16)
+ Frac128({0xa8811fb66d0faf7aULL, 0xbd08a39f580c36beULL}), // 2^(9/16)
+ Frac128({0x3e2ad0c964dd9f37ULL, 0xc5672a115506daddULL}), // 2^(10/16)
+ Frac128({0xe235838f95f2c6edULL, 0xce248c151f8480e3ULL}), // 2^(11/16)
+ Frac128({0x39a68bb9902d3fdeULL, 0xd744fccad69d6af4ULL}), // 2^(12/16)
+ Frac128({0x065895048dd333caULL, 0xe0ccdeec2a94e111ULL}), // 2^(13/16)
+ Frac128({0xd02d75b3706e54fbULL, 0xeac0c6e7dd24392eULL}), // 2^(14/16)
+ Frac128({0x7b9d0c7aed980fc3ULL, 0xf5257d152486cc2cULL}), // 2^(15/16)
+};
+
+// Polynomial approximation for (2^x - 1) / x:
+// For x in [0, 2^-4], we approximate:
+// P(x) ~ (2^x - 1) / x
+// then:
+// 2^x = 1 + x * P(x).
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > prec = 500;
+// > P = fpminimax((2^x - 1)/x, 14, [|128...|], [2^(-150), 2^(-4)], fixed);
+// > for i from 0 to 14 do {
+// c = coeff(P, i);
+// c_int = round(c * 2^128, 256, RN);
+// hi = floor(c_int * 2^(-64));
+// lo = c_int - hi * 2^64;
+// print("Frac128({", lo, "ULL,", hi, "ULL}), // c_", i);
+// };
+// > dirtyinfnorm(2^x - (1 + x * P), [0, 2^(-4)]);
+// 0x1.278a...p-139 < 2^-138.
+//
+// Expanding P(x) in Horner form:
+// P(x) = c_0 + x * (c_1 + x * (c_2 + ...))
+// where c_k ~ (log(2))^(k+1) / (k+1)! are all positive and in (0, 1).
+LIBC_INLINE_VAR constexpr Frac128 EXP2_POLY_128[15] = {
+ // lsb = 2^-128
+ Frac128({0xc9e3b39803f2f6afULL, 0xb17217f7d1cf79abULL}), // c_0
+ Frac128({0xde2d60dd92e6c1ceULL, 0x3d7f7bff058b1d50ULL}), // c_1
+ Frac128({0x99d3b15d9958cfa0ULL, 0x0e35846b82505fc5ULL}), // c_2
+ Frac128({0x39977c16aee20457ULL, 0x0276556df749cee5ULL}), // c_3
+ Frac128({0x41c5fda1c13b8fb1ULL, 0x005761ff9e299cc4ULL}), // c_4
+ Frac128({0xb7a58759e18c05e7ULL, 0x000a184897c363c3ULL}), // c_5
+ Frac128({0x3434f0e418799c8bULL, 0x0000ffe5fe2c4586ULL}), // c_6
+ Frac128({0x241c41f55c65febbULL, 0x0000162c0223a5c8ULL}), // c_7
+ Frac128({0x77e77364ac415ce2ULL, 0x000001b5253d395eULL}), // c_8
+ Frac128({0xfe79c174054850e2ULL, 0x0000001e4cf5158bULL}), // c_9
+ Frac128({0x29d719d4d61d89c9ULL, 0x00000001e8cac72dULL}), // c_10
+ Frac128({0xd15314bfa6b56886ULL, 0x000000001c3bd6b4ULL}), // c_11
+ Frac128({0xd8de7d5b7789dd3aULL, 0x0000000001815e5dULL}), // c_12
+ Frac128({0x491912c310d5531bULL, 0x0000000000132403ULL}), // c_13
+ Frac128({0x8d1d877caeba1cfcULL, 0x000000000000c220ULL}), // c_14
+};
+
+struct ExpBase192 {
+ int base_exp;
+ UInt<192> log2_b; // lsb = 2^(base_exp - 191)
+};
+
+// Evaluate 2^(idx * 2^-4 + u) in [1, 2) as Frac128 (lsb = 2^-127) for
+// idx in [0, 15] and u in [0, 2^-4) (lsb = 2^-128).
+LIBC_INLINE Frac128 exp2_eval_frac128(unsigned idx, const Frac128 &u) {
+ // p ~ (2^u - 1) / u
+ Frac128 p = EXP2_POLY_128[14]; // lsb = 2^-128
+ LIBC_LOOP_NOUNROLL
+ for (int i = 13; i >= 0; --i)
+ p = p * u + EXP2_POLY_128[i]; // lsb = 2^-128
+
+ // Reconstruction:
+ // u * p ~ 2^u - 1 (lsb = 2^-128)
+ // m = 2^(idx * 2^-4) (lsb = 2^-127)
+ // m + m * (u * p) ~ m * 2^u (lsb = 2^-127)
+ Frac128 m = EXP2_MID_FRAC128[idx]; // lsb = 2^-127
+ return m + m * (u * p); // lsb = 2^-127
+}
+
+// Compute 2^z in 128-bit precision for a DyadicFloat<128> z with
+// 2^-55 < |z| < 2^11 (z.exponent in [-181, -117]) and optional extra low limb
+// lo (lsb = 2^(z.exponent - 64)):
+// z = k * 2^-4 + u, 0 <= u < 2^-4
+// hi = k >> 4
+// idx = k & 0x0f
+// 2^z = 2^hi * EXP2_MID_FRAC128[idx] * (1 + u * P(u)).
+LIBC_INLINE DFloat128 exp2_df128(const DFloat128 &z, uint64_t lo = 0) {
+ // z.mantissa has lsb = 2^z.exponent, so the 2^-4 bit of |z| is at bit
+ // (-4 - z.exponent) of z.mantissa, leaving
+ // shift = 128 - (-4 - z.exponent) = z.exponent + 132 bits
+ // at or above 2^-4 in z.mantissa.val[1].
+ // After range reduction, 2^-55 < |z| < 2^11, so z.exponent is in
+ // [-181, -117] and shift is in [-49, 15].
+ int shift = z.exponent + 132;
+ int k = 0;
+ UInt<128> frac;
+ if (LIBC_LIKELY(shift > 0)) {
+ // |z| >= 2^-4: split |z| = k * 2^-4 + frac * 2^-132.
+ // Here shift is in [1, 15] and shr = 64 - shift is in [49, 63], so all
+ // 64-bit shifts below are strictly in [1, 63] (no undefined behavior).
+ int shr = 64 - shift; // in [49, 63]
+ // k = floor(|z| * 2^4)
+ k = static_cast<int>(z.mantissa.val[1] >> shr); // lsb = 2^-4
+ // Shift z.mantissa left by shift and fill the low shift bits from lo:
+ frac = UInt<128>({(z.mantissa.val[0] << shift) | (lo >> shr),
+ (z.mantissa.val[1] << shift) |
+ (z.mantissa.val[0] >> shr)}); // lsb = 2^-132
+ } else {
+ // |z| < 2^-4: k = 0 and frac = |z| * 2^132, with -shift in [0, 49].
+ frac = (-shift < 128) ? (z.mantissa >> (-shift))
+ : UInt<128>(0); // lsb = 2^-132
+ }
+
+ if (z.sign.is_neg()) {
+ if (LIBC_UNLIKELY(frac.is_zero())) {
+ k = -k;
+ } else {
+ // One's-complement negation on frac (differs from two's-complement by
+ // 2^-132, which is 1/16 ULP of u):
+ // -|z| ~ (-k - 1) * 2^-4 + (2^-4 - (frac + 1) * 2^-132).
+ frac = ~frac;
+ k = ~k;
+ }
+ }
+
+ // Decompose k * 2^-4 = hi + idx * 2^-4 and scale frac from lsb = 2^-132 to
+ // u in [0, 2^-4) with lsb = 2^-128:
+ int hi = k >> 4; // lsb = 1
+ unsigned idx = static_cast<unsigned>(k & 0x0f); // lsb = 2^-4
+ Frac128 u((frac >> 4).val); // lsb = 2^-128
+ Frac128 m_final = exp2_eval_frac128(idx, u); // lsb = 2^-127
+
+ // Since m_final is in [1, 2) with lsb = 2^-127, its leading bit is already
+ // at bit 127, so r is normalized with exponent = hi - 127:
+ DFloat128 r;
+ r.sign = Sign::POS;
+ r.exponent = hi - 127;
+ r.mantissa = UInt<128>(m_final.val);
+ return r;
+}
+
+// Compute b^x (for b = 2, e, or 10) in 128-bit precision using direct base-2
+// 192-bit fixed-point range reduction. Here log2_base->log2_b holds
+// round(log2(b) * 2^(191 - base_exp)), or log2_base is nullptr when b = 2.
+LIBC_NOINLINE LIBC_INLINE double
+exp_accurate_128(double x, const ExpBase192 *log2_base = nullptr) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ int exp_x = xbits.get_exponent();
+
+ // Normalize |x|'s 53-bit mantissa to [2^63, 2^64):
+ uint64_t mant_x = (xbits.get_mantissa() | (1ULL << 52)) << 11;
+
+ DFloat128 z;
+ z.sign = xbits.sign();
+ uint64_t w0 = 0;
+
+ if (log2_base != nullptr) {
+ // 1-limb x 3-limb product:
+ // mant_x has lsb = 2^(exp_x - 63)
+ // log2_base->log2_b has lsb = 2^(base_exp - 191)
+ // p = mant_x * log2_base->log2_b has lsb = 2^(exp_x + base_exp - 254)
+ UInt<256> p = log2_base->log2_b.ful_mul(UInt<64>(mant_x));
+ // Extract the top 192 bits (p.val[3], p.val[2], p.val[1]) as
+ // (z.mantissa, w0):
+ z.exponent = exp_x + log2_base->base_exp - 126;
+ z.mantissa = UInt<128>({p.val[2], p.val[3]}); // lsb = 2^z.exponent
+ w0 = p.val[1]; // lsb = 2^(z.exponent - 64)
+ } else {
+ // For b = 2, log2(b) = 1, so z = x directly:
+ z.exponent = exp_x - 127;
+ z.mantissa = UInt<128>({0, mant_x}); // lsb = 2^(exp_x - 127)
+ }
+
+ return static_cast<double>(exp2_df128(z, w0));
+}
+
+} // namespace double_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_ACCURATE_128_H
diff --git a/libc/src/__support/math/exp2_double_eval.h b/libc/src/__support/math/exp2_double_eval.h
new file mode 100644
index 0000000000000..5263575680cb4
--- /dev/null
+++ b/libc/src/__support/math/exp2_double_eval.h
@@ -0,0 +1,133 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// Double-precision implementation of exp2.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_DOUBLE_EVAL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_DOUBLE_EVAL_H
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#include "src/__support/math/exp_double_utils.h"
+#else
+#include "src/__support/math/exp_utils.h"
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace double_eval {
+
+LIBC_ALWAYS_INLINE double exp2(double x) {
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ uint64_t x_u = xbits.uintval();
+ uint64_t x_abs = xbits.abs().uintval();
+
+ // Bound constants (with biased exponent and mantissa bits):
+ // - Tiny bound: nextafter(2^-54, +inf) = 0x1.0000000000001p-54
+ // (since |x| <= 2^-54 is handled by 1.0 + x in set_exceptional)
+ constexpr uint64_t TINY_BOUND = 0x3c90'0000'0000'0001ULL;
+ // - Overflow bound: 1024.0
+ constexpr uint64_t OVERFLOW_BOUND = 0x4090'0000'0000'0000ULL;
+ // -1024.0
+ constexpr uint64_t NEG_OVERFLOW_BOUND = 0xc090'0000'0000'0000ULL;
+ // - Underflow bound: -1075.0
+ constexpr uint64_t UNDERFLOW_BOUND = 0xc090'cc00'0000'0000ULL;
+
+ // Normal range: 2^-54 <= |x| < 1024.
+ if (LIBC_UNLIKELY(x_abs - TINY_BOUND >= OVERFLOW_BOUND - TINY_BOUND)) {
+ // Exclude non-underflowing negative inputs in (-1075, -1024].
+ if (LIBC_LIKELY(x_u < NEG_OVERFLOW_BOUND || x_u >= UNDERFLOW_BOUND))
+ return set_exceptional(x);
+ }
+
+ // Now -1075 < x < -2^-54 or 2^-54 < x < 1024.
+
+ // Range reduction:
+ // Let x = k * 2^-N + dx
+ // in which:
+ // k = round(x * 2^N)
+ // hi = k >> N
+ // idx = k & (2^N - 1)
+ // then:
+ // 2^x = 2^hi * 2^(idx / 2^N) * 2^dx.
+
+ double kd = fputil::nearest_integer(x * EXP_SCALE);
+ int k = static_cast<int>(kd);
+
+ int hi = k >> EXP_N_BITS;
+ unsigned idx = static_cast<unsigned>(k & EXP_MASK);
+
+ double dx = fputil::multiply_add(kd, -EXP_SCALE_INV, x); // exact
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ double u = dx * LOG_2.hi;
+ return exp2_eval(x, u, hi, idx, /*is_exp2=*/true);
+#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ DoubleDouble u = fputil::exact_mult(dx, LOG_2.hi);
+ u.lo = fputil::multiply_add(dx, LOG_2.lo, u.lo);
+#else // !LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ // Without FMA, we use fixed-point splitting to avoid Veltkamp's splitting:
+ // Let:
+ // m = 2^(idx / 2^7) = m_hi + m_lo
+ // u = dx * log(2) = u_hi + u_lo
+ //
+ // In exp2_eval, we approximate 2^(idx / 2^7) * 2^dx = m * exp(u) by:
+ // m * (1 + u + u^2 * P(u)) ~ r_hi + r_lo
+ // where:
+ // r_hi = m_hi + m_hi * u_hi
+ // r_lo = m * (u_lo + u^2 * P(u)) + m_lo * (1 + u_hi)
+ //
+ // Since r_hi is in [0.997, 2.0), for r_hi to be exact in double precision,
+ // we need:
+ // lsb(m_hi * u_hi) >= 2^-52.
+ // We choose:
+ // lsb(m_hi) >= 2^-22 => |m_lo| <= 2^-23
+ // lsb(u_hi) >= 2^-30
+ // so that lsb(m_hi * u_hi) >= 2^-52, and the rounding errors from
+ // m_lo * (1 + u_hi) are bounded by 2^(-23 - 52) = 2^-75.
+ //
+ // To get lsb(u_hi) >= 2^-30 from u = dx * log(2), since |dx| <= 2^-8 and
+ // log(2) < 1, we split dx at 2^-18 and log(2) at 2^-12:
+ // dx = dx_hi + dx_lo (|dx_lo| <= 2^-19)
+ // log(2) = LOG_2_HI_12 + LOG_2_LO_12 (|LOG_2_LO_12| < 2^-13)
+ // Then:
+ // u_hi = dx_hi * LOG_2_HI_12
+ // is an exact multiple of 2^-30, and the two remainder products in:
+ // u_lo = dx_lo * LOG_2_HI_12 + dx * LOG_2_LO_12
+ // are both bounded by 2^-19, keeping the rounding errors of u_lo below 2^-71.
+ double dx_hi = (dx + SPLIT_18) - SPLIT_18; // lsb = 2^-18
+ double dx_lo = dx - dx_hi;
+ double u_hi = dx_hi * LOG_2_HI_12; // exact, lsb = 2^-30
+ double t0 = dx_lo * LOG_2_HI_12;
+ double t1 = dx * LOG_2_LO_12;
+ double u_lo = t0 + t1;
+ DoubleDouble u = {u_lo, u_hi};
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ return exp2_eval(x, u, hi, idx);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+}
+
+} // namespace double_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_DOUBLE_EVAL_H
diff --git a/libc/src/__support/math/exp_double_utils.h b/libc/src/__support/math/exp_double_utils.h
new file mode 100644
index 0000000000000..9a14c247165fd
--- /dev/null
+++ b/libc/src/__support/math/exp_double_utils.h
@@ -0,0 +1,175 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Shared double-precision utilities and tables for exponential functions when
+/// LIBC_MATH_HAS_SKIP_ACCURATE_PASS is enabled.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_UTILS_H
+
+#include "hdr/errno_macros.h"
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/exp_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace double_eval {
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+#ifdef LIBC_MATH_HAS_SMALL_TABLES
+
+LIBC_INLINE_VAR constexpr int EXP_N_BITS = 3;
+LIBC_INLINE_VAR constexpr double EXP_SCALE = 8.0;
+LIBC_INLINE_VAR constexpr double EXP_SCALE_INV = 0x1.0p-3;
+LIBC_INLINE_VAR constexpr int EXP_MASK = 0x07;
+
+// Degree-5 polynomial approximation for (exp(x) - 1 - x) / x^2 on
+// [-2^-4 * log(2), 2^-4 * log(2)]:
+// Generated by Sollya with:
+// > display = hexadecimal!;
+// > prec = 500;
+// > P = fpminimax((exp(x) - 1 - x)/x^2, 5, [|D...|],
+// [-2^(-4) * log(2), 2^(-4) * log(2)]);
+// > dirtyinfnorm(exp(x) - (1 + x + x^2 * P),
+// [-2^(-4) * log(2), 2^(-4) * log(2)]);
+// 0x1.6828...p-57 < 1.5 * 2^-57.
+LIBC_INLINE_VAR constexpr double EXP_P_COEFFS[6] = {
+ 0x1.000000000002ep-1, // p_0
+ 0x1.555555555549cp-3, // p_1
+ 0x1.5555554e95a5p-5, // p_2
+ 0x1.11111114c1e11p-7, // p_3
+ 0x1.6c1b70af254c4p-10, // p_4
+ 0x1.a01b62c901438p-13, // p_5
+};
+
+// Evaluate P(u) ~ (e^u - 1 - u) / u^2 using Estrin's scheme:
+LIBC_ALWAYS_INLINE double exp_poly(double u, double u2) {
+ double c0 = fputil::multiply_add(u, EXP_P_COEFFS[1], EXP_P_COEFFS[0]);
+ double c1 = fputil::multiply_add(u, EXP_P_COEFFS[3], EXP_P_COEFFS[2]);
+ double u4 = u2 * u2;
+ double c2 = fputil::multiply_add(u, EXP_P_COEFFS[5], EXP_P_COEFFS[4]);
+ double c01 = fputil::multiply_add(u2, c1, c0);
+ return fputil::multiply_add(u4, c2, c01);
+}
+
+// Lookup table for 2^(k * 2^-3) with k = 0..7.
+// Generated by Sollya with:
+// > display = hexadecimal!;
+// > for i from 0 to 7 do print(round(2^(i * 2^-3), D, RN), ",");
+LIBC_INLINE_VAR constexpr double EXP2_MID[8] = {
+ 0x1.0000000000000p0, 0x1.172b83c7d517bp0, 0x1.306fe0a31b715p0,
+ 0x1.4bfdad5362a27p0, 0x1.6a09e667f3bcdp0, 0x1.8ace5422aa0dbp0,
+ 0x1.ae89f995ad3adp0, 0x1.d5818dcfba487p0,
+};
+
+#else // !LIBC_MATH_HAS_SMALL_TABLES
+
+// Lookup table for 2^(k * 2^-7) with k = 0..127.
+// Generated by Sollya with:
+// > display = hexadecimal!;
+// > for i from 0 to 127 do print(round(2^(i * 2^-7), D, RN), ",");
+LIBC_INLINE_VAR constexpr double EXP2_MID[128] = {
+ 0x1.0000000000000p0, 0x1.0163da9fb3335p0, 0x1.02c9a3e778061p0,
+ 0x1.04315e86e7f85p0, 0x1.059b0d3158574p0, 0x1.0706b29ddf6dep0,
+ 0x1.0874518759bc8p0, 0x1.09e3ecac6f383p0, 0x1.0b5586cf9890fp0,
+ 0x1.0cc922b7247f7p0, 0x1.0e3ec32d3d1a2p0, 0x1.0fb66affed31bp0,
+ 0x1.11301d0125b51p0, 0x1.12abdc06c31ccp0, 0x1.1429aaea92dep0,
+ 0x1.15a98c8a58e51p0, 0x1.172b83c7d517bp0, 0x1.18af9388c8deap0,
+ 0x1.1a35beb6fcb75p0, 0x1.1bbe084045cd4p0, 0x1.1d4873168b9aap0,
+ 0x1.1ed5022fcd91dp0, 0x1.2063b88628cd6p0, 0x1.21f49917ddc96p0,
+ 0x1.2387a6e756238p0, 0x1.251ce4fb2a63fp0, 0x1.26b4565e27cddp0,
+ 0x1.284dfe1f56381p0, 0x1.29e9df51fdee1p0, 0x1.2b87fd0dad99p0,
+ 0x1.2d285a6e4030bp0, 0x1.2ecafa93e2f56p0, 0x1.306fe0a31b715p0,
+ 0x1.32170fc4cd831p0, 0x1.33c08b26416ffp0, 0x1.356c55f929ff1p0,
+ 0x1.371a7373aa9cbp0, 0x1.38cae6d05d866p0, 0x1.3a7db34e59ff7p0,
+ 0x1.3c32dc313a8e5p0, 0x1.3dea64c123422p0, 0x1.3fa4504ac801cp0,
+ 0x1.4160a21f72e2ap0, 0x1.431f5d950a897p0, 0x1.44e086061892dp0,
+ 0x1.46a41ed1d0057p0, 0x1.486a2b5c13cdp0, 0x1.4a32af0d7d3dep0,
+ 0x1.4bfdad5362a27p0, 0x1.4dcb299fddd0dp0, 0x1.4f9b2769d2ca7p0,
+ 0x1.516daa2cf6642p0, 0x1.5342b569d4f82p0, 0x1.551a4ca5d920fp0,
+ 0x1.56f4736b527dap0, 0x1.58d12d497c7fdp0, 0x1.5ab07dd485429p0,
+ 0x1.5c9268a5946b7p0, 0x1.5e76f15ad2148p0, 0x1.605e1b976dc09p0,
+ 0x1.6247eb03a5585p0, 0x1.6434634ccc32p0, 0x1.6623882552225p0,
+ 0x1.68155d44ca973p0, 0x1.6a09e667f3bcdp0, 0x1.6c012750bdabfp0,
+ 0x1.6dfb23c651a2fp0, 0x1.6ff7df9519484p0, 0x1.71f75e8ec5f74p0,
+ 0x1.73f9a48a58174p0, 0x1.75feb564267c9p0, 0x1.780694fde5d3fp0,
+ 0x1.7a11473eb0187p0, 0x1.7c1ed0130c132p0, 0x1.7e2f336cf4e62p0,
+ 0x1.80427543e1a12p0, 0x1.82589994cce13p0, 0x1.8471a4623c7adp0,
+ 0x1.868d99b4492edp0, 0x1.88ac7d98a6699p0, 0x1.8ace5422aa0dbp0,
+ 0x1.8cf3216b5448cp0, 0x1.8f1ae99157736p0, 0x1.9145b0b91ffc6p0,
+ 0x1.93737b0cdc5e5p0, 0x1.95a44cbc8520fp0, 0x1.97d829fde4e5p0,
+ 0x1.9a0f170ca07bap0, 0x1.9c49182a3f09p0, 0x1.9e86319e32323p0,
+ 0x1.a0c667b5de565p0, 0x1.a309bec4a2d33p0, 0x1.a5503b23e255dp0,
+ 0x1.a799e1330b358p0, 0x1.a9e6b5579fdbfp0, 0x1.ac36bbfd3f37ap0,
+ 0x1.ae89f995ad3adp0, 0x1.b0e07298db666p0, 0x1.b33a2b84f15fbp0,
+ 0x1.b59728de5593ap0, 0x1.b7f76f2fb5e47p0, 0x1.ba5b030a1064ap0,
+ 0x1.bcc1e904bc1d2p0, 0x1.bf2c25bd71e09p0, 0x1.c199bdd85529cp0,
+ 0x1.c40ab5fffd07ap0, 0x1.c67f12e57d14bp0, 0x1.c8f6d9406e7b5p0,
+ 0x1.cb720dcef9069p0, 0x1.cdf0b555dc3fap0, 0x1.d072d4a07897cp0,
+ 0x1.d2f87080d89f2p0, 0x1.d5818dcfba487p0, 0x1.d80e316c98398p0,
+ 0x1.da9e603db3285p0, 0x1.dd321f301b46p0, 0x1.dfc97337b9b5fp0,
+ 0x1.e264614f5a129p0, 0x1.e502ee78b3ff6p0, 0x1.e7a51fbc74c83p0,
+ 0x1.ea4afa2a490dap0, 0x1.ecf482d8e67f1p0, 0x1.efa1bee615a27p0,
+ 0x1.f252b376bba97p0, 0x1.f50765b6e454p0, 0x1.f7bfdad9cbe14p0,
+ 0x1.fa7c1819e90d8p0, 0x1.fd3c22b8f71f1p0,
+};
+
+#endif // LIBC_MATH_HAS_SMALL_TABLES
+
+LIBC_NOINLINE LIBC_INLINE double exp_denorm(double x, int hi, double r,
+ bool is_exp2) {
+ using FPBits = typename fputil::FPBits<double>;
+ if (is_exp2 && LIBC_UNLIKELY(static_cast<double>(hi) == x))
+ return FPBits(uint64_t(1) << (hi + 1074)).get_val();
+
+ double res = ziv_test_denorm</*SKIP_ZIV_TEST=*/true>(hi, r, 0.0, 0.0).value();
+ if (LIBC_UNLIKELY(FPBits(res).is_normal()))
+ return res;
+
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_underflow_except_if_required<double>();
+ return res;
+}
+
+#ifdef LIBC_MATH_HAS_SMALL_TABLES
+LIBC_NOINLINE LIBC_INLINE double
+#else
+LIBC_ALWAYS_INLINE double
+#endif // LIBC_MATH_HAS_SMALL_TABLES
+exp2_eval(double x, double u, int hi, unsigned idx, bool is_exp2 = false) {
+ double exp2_mid = EXP2_MID[idx];
+ double u2 = u * u;
+ double p = exp_poly(u, u2);
+ // u + u^2 * p
+ double poly_lo = fputil::multiply_add(u2, p, u);
+ double r = fputil::multiply_add(exp2_mid, poly_lo, exp2_mid);
+
+ if (LIBC_UNLIKELY(hi <= -1022))
+ return exp_denorm(x, hi, r, is_exp2);
+
+ return scale_exp2(r, hi);
+}
+
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+} // namespace double_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_UTILS_H
diff --git a/libc/src/__support/math/exp_utils.h b/libc/src/__support/math/exp_utils.h
index e46ede2ed8e26..a7c96ba680661 100644
--- a/libc/src/__support/math/exp_utils.h
+++ b/libc/src/__support/math/exp_utils.h
@@ -1,17 +1,44 @@
-//===-- Common utils for exp function ---------------------------*- C++ -*-===//
+//===----------------------------------------------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
//
//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Shared utilities and tables for double-precision exponential functions.
+///
+//===----------------------------------------------------------------------===//
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
#include "src/__support/CPP/bit.h"
#include "src/__support/CPP/optional.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/rounding_mode.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/architectures.h"
+#include "src/__support/macros/properties/cpu_features.h"
+
+#ifdef LIBC_TARGET_CPU_HAS_SSE2
+#include <immintrin.h>
+#elif defined(LIBC_TARGET_ARCH_IS_AARCH64) && defined(__ARM_NEON)
+#include <arm_neon.h>
+#endif // LIBC_TARGET_CPU_HAS_SSE2
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#include "src/__support/math/exp2_accurate_128.h"
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
namespace LIBC_NAMESPACE_DECL {
@@ -68,6 +95,452 @@ ziv_test_denorm(int hi, double mid, double lo, double err) {
}
}
+namespace math {
+namespace double_eval {
+
+using fputil::DoubleDouble;
+
+// log(2) in double-double precision, generated by Sollya with:
+// > display = hexadecimal!;
+// > prec = 500;
+// > hi = round(log(2), D, RN);
+// > lo = round(log(2) - hi, D, RN);
+// > print("{", lo, ",", hi, "};");
+LIBC_INLINE_VAR constexpr DoubleDouble LOG_2 = {0x1.abc9e3b39803fp-56,
+ 0x1.62e42fefa39efp-1};
+
+// Directly add an integer `hi` to the exponent field of a double precision `r`.
+// It is equivalent to scaling `r` by 2^hi when both `r` and the resulting bits
+// are normal.
+// The default implementation is simply bit cast `r` to int64`, shift `hi` to
+// the exponent field, then add `r` bits and shifted `hi`, before bit cast the
+// result back to double.
+// The addition is done in integer, and most likely generated as addition in the
+// general purpose registers (GPR), causing a full round trip from
+// floating-point/SIMD register -> GPR -> FP/SIMD.
+// When SSE2 or NEON are available, we can use the compiler intrinsic to keep
+// `r` in the FP/SIMD register, and perform the integer addition in the FP/SIMD
+// registers instead. This in general results in a reduction in latency since
+// `r` is normally on the critical path.
+LIBC_ALWAYS_INLINE double scale_exp2(double r, int hi) {
+ using FPBits = typename fputil::FPBits<double>;
+#ifdef LIBC_TARGET_CPU_HAS_SSE2
+ // Move `hi` to SIMD register and shift to exponent field position.
+ __m128i v_exp = _mm_slli_epi64(_mm_cvtsi32_si128(hi), FPBits::FRACTION_LEN);
+ // Load `r` and "bit_cast" it to integer SIMD type.
+ __m128i v_r = _mm_castpd_si128(_mm_set_sd(r));
+ // Perform integer addition in SIMD registers and bit cast the result back to
+ // double.
+ return _mm_cvtsd_f64(_mm_castsi128_pd(_mm_add_epi64(v_exp, v_r)));
+#elif defined(LIBC_TARGET_ARCH_IS_AARCH64) && defined(__ARM_NEON)
+ // Shift `hi` to exponent field position and move it to SIMD register.
+ int64x1_t v_exp =
+ vdup_n_s64(static_cast<int64_t>(hi) << FPBits::FRACTION_LEN);
+ // Bit cast `r` to integer SIMD type.
+ int64x1_t v_r = vreinterpret_s64_f64(vdup_n_f64(r));
+ // Perform integer addition in SIMD registers and bit cast the result back to
+ // double.
+ return vget_lane_f64(vreinterpret_f64_s64(vadd_s64(v_exp, v_r)), 0);
+#else
+ int64_t exp_hi = static_cast<int64_t>(hi) * (1LL << FPBits::FRACTION_LEN);
+ return cpp::bit_cast<double>(exp_hi + cpp::bit_cast<int64_t>(r));
+#endif // LIBC_TARGET_CPU_HAS_SSE2
+}
+
+// Check for exceptional cases (tiny, overflow, underflow, inf, or nan) for
+// exp, exp2, and exp10. Callers have already verified that x is in the exact
+// tiny, underflow, or overflow/inf/nan range for the specific function before
+// calling this helper.
+LIBC_NOINLINE LIBC_INLINE double set_exceptional(double x) {
+ using FPBits = typename fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ // |x| <= 0x1.0p-53 (caller has already checked the exact tiny bound):
+ constexpr uint64_t TINY_BOUND = 0x3ca0'0000'0000'0000ULL;
+ if (xbits.abs().uintval() <= TINY_BOUND) {
+ // b^x ~ 1 + x
+ return 1.0 + x;
+ }
+
+ // Negative underflow or -inf/nan (caller has already excluded
+ // non-underflowing negative inputs):
+ if (xbits.is_neg()) {
+ // b^(-Inf) = 0
+ if (xbits.is_inf())
+ return 0.0;
+
+ // b^(nan) = nan
+ // Adding 0.0 to silence signaling nan.
+ if (xbits.is_nan())
+ return x + 0.0;
+
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_underflow_except_if_required<double>();
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (fputil::quick_get_round() == FE_UPWARD)
+ return FPBits::min_subnormal().get_val();
+#endif
+ return 0.0;
+ }
+
+ // Positive overflow or +inf/nan (caller has already checked the exact
+ // overflow bound)
+
+ // x is finite and positive
+ if (xbits.uintval() < FPBits::inf().uintval()) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_overflow_except_if_required<double>();
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ int rounding = fputil::quick_get_round();
+ if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
+ return FPBits::max_normal().get_val();
+#endif
+ return FPBits::inf().get_val();
+ }
+ // x is +inf or nan
+ // Adding `inf` to silence signaling nan.
+ return x + FPBits::inf().get_val();
+}
+
+#if !defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) || \
+ !defined(LIBC_MATH_HAS_SMALL_TABLES)
+
+LIBC_INLINE_VAR constexpr int EXP_N_BITS = 7;
+LIBC_INLINE_VAR constexpr double EXP_SCALE = 128.0; // 2^EXP_N_BITS = 2^7
+LIBC_INLINE_VAR constexpr double EXP_SCALE_INV = 0x1.0p-7;
+LIBC_INLINE_VAR constexpr int EXP_MASK = 0x7f;
+
+// Polynomial approximation for (exp(x) - 1 - x) / x^2 on
+// [-B, B] with B = 2^-8 * log(2):
+//
+// 1) A degree-3 minimax polynomial P0 for (exp(x) - 1 - x) / x^2 on [-B, B]
+// (which minimizes the error of exp(x) - (1 + x + x^2 * P0(x)) weighted by
+// 1 / x^2, i.e., relative error with respect to x^2):
+// > display = hexadecimal!;
+// > prec = 500;
+// > B = 2^(-8) * log(2);
+// > P0 = fpminimax((exp(x) - 1 - x)/x^2, 3, [|D...|], [-B, B]);
+// > dirtyinfnorm(exp(x) - (1 + x + x^2 * P0), [-B, B]);
+// 0x1.4362...p-64 < 2^-63.66.
+//
+// Some explanation for the errors (from Chebyshev polynomials):
+// Because P0 minimizes |(exp(x) - 1 - x) / x^2 - P0(x)|, its error is governed
+// by the degree-4 Chebyshev polynomial:
+// T_4(x / B) = 8*(x / B)^4 - 8*(x / B)^2 + 1
+// approximating the leading omitted Taylor term x^4 / 720, with magnitude:
+// B^4 / (720 * 8).
+// Multiplying by x^2 gives an absolute error of:
+// x^2 * B^4 / (720 * 8) * T_4(x / B),
+// peaking at x = +-B with:
+// B^6 / (720 * 8) ~ 0x1.4309...p-64 < 2^-63.66.
+//
+// 2) If we instead minimize the absolute error of exp(x) - (1 + x + P(x))
+// directly on [-B, B] with monomials [|2, 3, 4, 5|]:
+// > P = fpminimax(exp(x) - 1 - x, [|2, 3, 4, 5|], [|D...|], [-B, B],
+// absolute);
+// > dirtyinfnorm(exp(x) - (1 + x + P), [-B, B]);
+// 0x1.8e05...p-66 < 2^-65.36.
+//
+// Some explanation for the better errors from Chebyshev polynomials:
+// The free x^2 and x^4 coefficients in P(x) approximate the omitted x^6 / 720
+// term on [-B, B] subject to vanishing at x = 0. Setting t = (x / B)^2 in
+// [0, 1], this is equivalent to approximating t^3 on [0, 1] by
+// a_2 * t^2 + a_1 * t vanishing at t = 0.
+// The degree-3 Chebyshev polynomial T_3(z) = 4*z^3 - 3*z has a root at
+// z = -sqrt(3)/2 and three alternating extrema +-1 on [-sqrt(3)/2, 1] at
+// z = -1/2, 1/2, 1.
+// Mapping t in [0, 1] affinely onto z in [-sqrt(3)/2, 1] via:
+// z = ((2 + sqrt(3)) / 2) * t - sqrt(3) / 2
+// gives an error curve T_3(z) that vanishes at t = 0 and equioscillates three
+// times on (0, 1], with leading coefficient:
+// 4 * ((2 + sqrt(3)) / 2)^3 = (2 + sqrt(3))^3 / 2.
+// Dividing by this leading coefficient gives a maximum error on [-B, B] of:
+// (2 / (2 + sqrt(3))^3) * B^6 / 720 = 2 * (2 - sqrt(3))^3 * B^6 / 720
+// ~ 0x1.8dba...p-66 < 2^-65.36,
+// reducing the maximum error compared to P0 by (2 + sqrt(3))^3 / 16 ~ 3.25x.
+LIBC_INLINE_VAR constexpr double EXP_P_COEFFS[4] = {
+ 0x1.ffffffffffdbdp-2, // p_0
+ 0x1.555555555543dp-3, // p_1
+ 0x1.55555cf16d6eep-5, // p_2
+ 0x1.111116782220ap-7, // p_3
+};
+
+// Evaluate P(u) ~ (e^u - 1 - u) / u^2 using Estrin's scheme:
+LIBC_ALWAYS_INLINE double exp_poly(double u, double u2) {
+ double c0 = fputil::multiply_add(u, EXP_P_COEFFS[1], EXP_P_COEFFS[0]);
+ double c1 = fputil::multiply_add(u, EXP_P_COEFFS[3], EXP_P_COEFFS[2]);
+ return fputil::multiply_add(u2, c1, c0);
+}
+
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS || !LIBC_MATH_HAS_SMALL_TABLES
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// Fixed-point splitting constants (ulp = 2^-22, 2^-18) and pre-split
+// 12-bit + 53-bit logarithm constants for exact non-FMA multiplication without
+// Veltkamp splitting:
+LIBC_INLINE_VAR constexpr double SPLIT_22 = 0x1.8p30; // 1.5 * 2^(-22 + 52)
+LIBC_INLINE_VAR constexpr double SPLIT_18 = 0x1.8p34; // 1.5 * 2^(-18 + 52)
+
+// Splitting log(2) ~ hi + lo with 12-bit-precision hi part.
+// > hi = round(log(2), 12, RN);
+// > lo = round(log(2) - hi, D, RN);
+LIBC_INLINE_VAR constexpr double LOG_2_HI_12 = 0x1.62ep-1;
+LIBC_INLINE_VAR constexpr double LOG_2_LO_12 = 0x1.0bfbe8e7bcd5ep-15;
+
+// Error bound for double precision fast path (N = 7, Ziv failure rate <
+// 0.05%):
+LIBC_INLINE_VAR constexpr double EXP2_ERR_D = 0x1.0p-64;
+
+// Lookup table for 2^(i * 2^-7) with i = 0..127.
+// Generated by Sollya with:
+// > for i from 0 to 127 do {
+// a = 2^(i * 2^-7);
+// b = round(a, D, RN);
+// c = round(a - b, D, RN);
+// print("{", c, ",", b, "},");
+// };
+LIBC_INLINE_VAR constexpr DoubleDouble EXP2_MID[128] = {
+ {0, 0x1p0},
+ {0x1.b61299ab8cdb7p-54, 0x1.0163da9fb3335p0},
+ {-0x1.19083535b085dp-56, 0x1.02c9a3e778061p0},
+ {-0x1.0a31c1977c96ep-54, 0x1.04315e86e7f85p0},
+ {0x1.d73e2a475b465p-55, 0x1.059b0d3158574p0},
+ {-0x1.c91dfe2b13c27p-55, 0x1.0706b29ddf6dep0},
+ {0x1.186be4bb284ffp-57, 0x1.0874518759bc8p0},
+ {0x1.1487818316136p-54, 0x1.09e3ecac6f383p0},
+ {0x1.8a62e4adc610bp-54, 0x1.0b5586cf9890fp0},
+ {0x1.01edc16e24f71p-54, 0x1.0cc922b7247f7p0},
+ {0x1.03a1727c57b53p-59, 0x1.0e3ec32d3d1a2p0},
+ {-0x1.b9bedc44ebd7bp-57, 0x1.0fb66affed31bp0},
+ {-0x1.6c51039449b3ap-54, 0x1.11301d0125b51p0},
+ {-0x1.1b514b36ca5c7p-58, 0x1.12abdc06c31ccp0},
+ {-0x1.32fbf9af1369ep-54, 0x1.1429aaea92dep0},
+ {0x1.2406ab9eeab0ap-55, 0x1.15a98c8a58e51p0},
+ {-0x1.19041b9d78a76p-55, 0x1.172b83c7d517bp0},
+ {-0x1.11023d1970f6cp-54, 0x1.18af9388c8deap0},
+ {0x1.e5b4c7b4968e4p-55, 0x1.1a35beb6fcb75p0},
+ {-0x1.95386352ef607p-54, 0x1.1bbe084045cd4p0},
+ {0x1.e016e00a2643cp-54, 0x1.1d4873168b9aap0},
+ {-0x1.1df98027bb78cp-54, 0x1.1ed5022fcd91dp0},
+ {0x1.dc775814a8495p-55, 0x1.2063b88628cd6p0},
+ {0x1.2a97e9494a5eep-55, 0x1.21f49917ddc96p0},
+ {0x1.9b07eb6c70573p-54, 0x1.2387a6e756238p0},
+ {0x1.ac155bef4f4a4p-55, 0x1.251ce4fb2a63fp0},
+ {0x1.2bd339940e9d9p-55, 0x1.26b4565e27cddp0},
+ {-0x1.a4c3a8c3f0d7ep-54, 0x1.284dfe1f56381p0},
+ {0x1.612e8afad1255p-55, 0x1.29e9df51fdee1p0},
+ {-0x1.10adcd6381aa4p-59, 0x1.2b87fd0dad99p0},
+ {0x1.0024754db41d5p-54, 0x1.2d285a6e4030bp0},
+ {0x1.1ca0f45d52383p-56, 0x1.2ecafa93e2f56p0},
+ {0x1.6f46ad23182e4p-55, 0x1.306fe0a31b715p0},
+ {0x1.a9ce78e18047cp-55, 0x1.32170fc4cd831p0},
+ {0x1.32721843659a6p-54, 0x1.33c08b26416ffp0},
+ {-0x1.b5cee5c4e4628p-55, 0x1.356c55f929ff1p0},
+ {-0x1.63aeabf42eae2p-54, 0x1.371a7373aa9cbp0},
+ {-0x1.e958d3c9904bdp-54, 0x1.38cae6d05d866p0},
+ {-0x1.5e436d661f5e3p-56, 0x1.3a7db34e59ff7p0},
+ {-0x1.efff8375d29c3p-54, 0x1.3c32dc313a8e5p0},
+ {0x1.ada0911f09ebcp-55, 0x1.3dea64c123422p0},
+ {-0x1.7d023f956f9f3p-54, 0x1.3fa4504ac801cp0},
+ {-0x1.ef3691c309278p-58, 0x1.4160a21f72e2ap0},
+ {-0x1.1c7dde35f7999p-55, 0x1.431f5d950a897p0},
+ {0x1.89b7a04ef80dp-59, 0x1.44e086061892dp0},
+ {0x1.c944bd1648a76p-54, 0x1.46a41ed1d0057p0},
+ {0x1.3c1a3b69062fp-56, 0x1.486a2b5c13cdp0},
+ {0x1.9cb62f3d1be56p-54, 0x1.4a32af0d7d3dep0},
+ {0x1.d4397afec42e2p-56, 0x1.4bfdad5362a27p0},
+ {0x1.8ecdbbc6a7833p-54, 0x1.4dcb299fddd0dp0},
+ {-0x1.4b309d25957e3p-54, 0x1.4f9b2769d2ca7p0},
+ {-0x1.f768569bd93efp-55, 0x1.516daa2cf6642p0},
+ {-0x1.07abe1db13cadp-55, 0x1.5342b569d4f82p0},
+ {-0x1.d689cefede59bp-55, 0x1.551a4ca5d920fp0},
+ {0x1.9bb2c011d93adp-54, 0x1.56f4736b527dap0},
+ {0x1.295e15b9a1de8p-55, 0x1.58d12d497c7fdp0},
+ {0x1.6324c054647adp-54, 0x1.5ab07dd485429p0},
+ {0x1.c4b1b816986a2p-60, 0x1.5c9268a5946b7p0},
+ {0x1.ba6f93080e65ep-54, 0x1.5e76f15ad2148p0},
+ {-0x1.3e2429b56de47p-54, 0x1.605e1b976dc09p0},
+ {-0x1.383c17e40b497p-54, 0x1.6247eb03a5585p0},
+ {-0x1.c483c759d8933p-55, 0x1.6434634ccc32p0},
+ {-0x1.bb60987591c34p-54, 0x1.6623882552225p0},
+ {0x1.038ae44f73e65p-57, 0x1.68155d44ca973p0},
+ {-0x1.bdd3413b26456p-54, 0x1.6a09e667f3bcdp0},
+ {-0x1.2895667ff0b0dp-56, 0x1.6c012750bdabfp0},
+ {-0x1.bbe3a683c88abp-57, 0x1.6dfb23c651a2fp0},
+ {-0x1.83c0f25860ef6p-55, 0x1.6ff7df9519484p0},
+ {-0x1.16e4786887a99p-55, 0x1.71f75e8ec5f74p0},
+ {-0x1.0a8d96c65d53cp-54, 0x1.73f9a48a58174p0},
+ {-0x1.0245957316dd3p-54, 0x1.75feb564267c9p0},
+ {0x1.866b80a02162dp-54, 0x1.780694fde5d3fp0},
+ {-0x1.41577ee04992fp-55, 0x1.7a11473eb0187p0},
+ {0x1.f124cd1164dd6p-54, 0x1.7c1ed0130c132p0},
+ {0x1.05d02ba15797ep-56, 0x1.7e2f336cf4e62p0},
+ {-0x1.27c86626d972bp-54, 0x1.80427543e1a12p0},
+ {-0x1.d4c1dd41532d8p-54, 0x1.82589994cce13p0},
+ {-0x1.8d684a341cdfbp-55, 0x1.8471a4623c7adp0},
+ {-0x1.fc6f89bd4f6bap-54, 0x1.868d99b4492edp0},
+ {0x1.994c2f37cb53ap-54, 0x1.88ac7d98a6699p0},
+ {0x1.6e9f156864b27p-54, 0x1.8ace5422aa0dbp0},
+ {-0x1.0d55e32e9e3aap-56, 0x1.8cf3216b5448cp0},
+ {0x1.5cc13a2e3976cp-55, 0x1.8f1ae99157736p0},
+ {-0x1.dd6792e582524p-54, 0x1.9145b0b91ffc6p0},
+ {-0x1.75fc781b57ebcp-57, 0x1.93737b0cdc5e5p0},
+ {-0x1.64b7c96a5f039p-56, 0x1.95a44cbc8520fp0},
+ {-0x1.d185b7c1b85d1p-54, 0x1.97d829fde4e5p0},
+ {-0x1.173bd91cee632p-54, 0x1.9a0f170ca07bap0},
+ {0x1.c7c46b071f2bep-56, 0x1.9c49182a3f09p0},
+ {0x1.824ca78e64c6ep-56, 0x1.9e86319e32323p0},
+ {-0x1.359495d1cd533p-54, 0x1.a0c667b5de565p0},
+ {0x1.6305c7ddc36abp-54, 0x1.a309bec4a2d33p0},
+ {-0x1.d2f6edb8d41e1p-54, 0x1.a5503b23e255dp0},
+ {0x1.bcb7ecac563c7p-54, 0x1.a799e1330b358p0},
+ {0x1.0fac90ef7fd31p-54, 0x1.a9e6b5579fdbfp0},
+ {-0x1.f9234cae76cdp-55, 0x1.ac36bbfd3f37ap0},
+ {0x1.7a1cd345dcc81p-54, 0x1.ae89f995ad3adp0},
+ {-0x1.bdef54c80e425p-54, 0x1.b0e07298db666p0},
+ {-0x1.2805e3084d708p-57, 0x1.b33a2b84f15fbp0},
+ {-0x1.c71dfbbba6de3p-54, 0x1.b59728de5593ap0},
+ {-0x1.5584f7e54ac3bp-56, 0x1.b7f76f2fb5e47p0},
+ {-0x1.efcd30e54292ep-54, 0x1.ba5b030a1064ap0},
+ {0x1.23dd07a2d9e84p-55, 0x1.bcc1e904bc1d2p0},
+ {-0x1.efdca3f6b9c73p-54, 0x1.bf2c25bd71e09p0},
+ {0x1.11065895048ddp-55, 0x1.c199bdd85529cp0},
+ {0x1.b4537e083c60ap-54, 0x1.c40ab5fffd07ap0},
+ {0x1.2884dff483cadp-54, 0x1.c67f12e57d14bp0},
+ {0x1.1acbc48805c44p-56, 0x1.c8f6d9406e7b5p0},
+ {0x1.503cbd1e949dbp-56, 0x1.cb720dcef9069p0},
+ {-0x1.dd83b53829d72p-55, 0x1.cdf0b555dc3fap0},
+ {-0x1.cbc3743797a9cp-54, 0x1.d072d4a07897cp0},
+ {-0x1.d487b719d8578p-54, 0x1.d2f87080d89f2p0},
+ {0x1.2ed02d75b3707p-55, 0x1.d5818dcfba487p0},
+ {-0x1.11ec18beddfe8p-54, 0x1.d80e316c98398p0},
+ {0x1.c2300696db532p-54, 0x1.da9e603db3285p0},
+ {0x1.2da5778f018c3p-54, 0x1.dd321f301b46p0},
+ {-0x1.1a5cd4f184b5cp-54, 0x1.dfc97337b9b5fp0},
+ {-0x1.7b627817a1496p-54, 0x1.e264614f5a129p0},
+ {0x1.39e8980a9cc8fp-55, 0x1.e502ee78b3ff6p0},
+ {0x1.2d522ca0c8de2p-54, 0x1.e7a51fbc74c83p0},
+ {-0x1.e9c23179c2893p-54, 0x1.ea4afa2a490dap0},
+ {-0x1.c93f3b411ad8cp-54, 0x1.ecf482d8e67f1p0},
+ {0x1.dc7f486a4b6bp-54, 0x1.efa1bee615a27p0},
+ {0x1.3a1a5bf0d8e43p-54, 0x1.f252b376bba97p0},
+ {0x1.9d3e12dd8a18bp-54, 0x1.f50765b6e454p0},
+ {-0x1.dbb12d006350ap-54, 0x1.f7bfdad9cbe14p0},
+ {0x1.74853f3a5931ep-55, 0x1.fa7c1819e90d8p0},
+ {0x1.2eb74966579e7p-57, 0x1.fd3c22b8f71f1p0},
+};
+
+LIBC_NOINLINE LIBC_INLINE double
+exp_denorm(double x, int hi, double r_hi, double lo_u, double err_d,
+ const ExpBase192 *log2_base = nullptr) {
+ using FPBits = typename fputil::FPBits<double>;
+ // Exact subnormal power of 2 when x is an integer in [-1074, -1022].
+ if (log2_base == nullptr && LIBC_UNLIKELY(static_cast<double>(hi) == x))
+ return FPBits(uint64_t(1) << (hi + 1074)).get_val();
+
+ if (auto r_denorm = ziv_test_denorm(hi, r_hi, lo_u - err_d, err_d);
+ LIBC_LIKELY(r_denorm.has_value())) {
+ double res = r_denorm.value();
+ if (LIBC_UNLIKELY(FPBits(res).is_normal()))
+ return res;
+
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_underflow_except_if_required<double>();
+ return res;
+ }
+ return exp_accurate_128(x, log2_base);
+}
+
+LIBC_ALWAYS_INLINE double exp2_eval(double x, DoubleDouble u, int hi,
+ unsigned idx, double err_d = EXP2_ERR_D,
+ const ExpBase192 *log2_base = nullptr) {
+ DoubleDouble exp2_mid = EXP2_MID[idx];
+
+ // Reconstruction of 2^(idx / 2^7) * exp(u) ~ exp2_mid * (1 + u + u^2 * P(u))
+ // where:
+ // exp2_mid = 2^(idx / 2^7) = exp2_mid.hi + exp2_mid.lo
+ // u = dx * log(b) = u.hi + u.lo
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ // Decompose exp2_mid * (1 + u + u^2 * P(u)) ~ r_hi + (r_lo + tail):
+ // r_hi + r_lo = exp2_mid.hi * (1 + u.hi)
+ // tail = exp2_mid.hi * (u.lo + u^2 * P(u))
+ // + exp2_mid.lo * (1 + u.hi)
+ // where the omitted cross term exp2_mid.lo * (u.lo + u^2 * P(u)) is bounded
+ // by 2^(-54 - 18) = 2^-72 and absorbed into err_d.
+
+ // p ~ (exp(u) - 1 - u) / u^2
+ double u2 = u.hi * u.hi;
+ double p = exp_poly(u.hi, u2);
+
+ // r_hi = round(exp2_mid.hi * (1 + u.hi))
+ double r_hi = fputil::multiply_add(exp2_mid.hi, u.hi, exp2_mid.hi);
+ // mid_lo ~ exp2_mid.lo * (1 + u.hi)
+ double mid_lo = fputil::multiply_add(exp2_mid.lo, u.hi, exp2_mid.lo);
+ // poly_lo ~ u.lo + u^2 * P(u)
+ double poly_lo = fputil::multiply_add(u2, p, u.lo);
+
+ // Fold +/-err_d into mid_lo so tail_u, tail_l, and r_lo can be evaluated in
+ // parallel:
+ // mid_u = mid_lo + err_d
+ // mid_l = mid_lo - err_d
+ double mid_u = mid_lo + err_d;
+ double mid_l = mid_lo - err_d;
+ // tail_{u,l} ~ exp2_mid.hi * (u.lo + u^2 * P(u)) + exp2_mid.lo * (1 + u.hi)
+ // +/- err_d
+ double tail_u = fputil::multiply_add(exp2_mid.hi, poly_lo, mid_u);
+ double tail_l = fputil::multiply_add(exp2_mid.hi, poly_lo, mid_l);
+ // r_lo = exp2_mid.hi * (1 + u.hi) - r_hi (exact to within 2^-105)
+ double r_lo = fputil::multiply_add(exp2_mid.hi, u.hi, exp2_mid.hi - r_hi);
+
+ // lo_{u,l} ~ r_lo + tail +/- err_d
+ double lo_u = r_lo + tail_u;
+ double lo_l = r_lo + tail_l;
+#else
+ // p ~ (exp(u) - 1 - u) / u^2
+ double u_sum = u.hi + u.lo;
+ double u2 = u_sum * u_sum;
+ double p = exp_poly(u_sum, u2);
+
+ // Split exp2_mid = m_hi + m_lo with lsb(m_hi) >= 2^-22 and |m_lo| <= 2^-23:
+ double m_hi = (exp2_mid.hi + SPLIT_22) - SPLIT_22; // lsb = 2^-22
+ double m_lo = (exp2_mid.hi - m_hi) + exp2_mid.lo;
+
+ // Since lsb(u.hi) >= 2^-30 and r_hi is in [0.997, 2.0), both m_hi * u.hi and
+ // r_hi = m_hi * (1 + u.hi) are exact:
+ double prod_hi = m_hi * u.hi; // exact, lsb = 2^-52
+ double r_hi = m_hi + prod_hi; // exact
+
+ // poly_lo ~ u.lo + u^2 * P(u)
+ double u2_p = u2 * p;
+ double poly_lo = u.lo + u2_p;
+
+ // r_lo ~ exp2_mid.hi * (u.lo + u^2 * P(u)) + m_lo * (1 + u.hi)
+ double l0 = exp2_mid.hi * poly_lo;
+ double l1 = m_lo * (1.0 + u.hi);
+ double r_lo = l0 + l1;
+
+ // lo_{u,l} = r_lo +/- err_d
+ double lo_u = r_lo + err_d;
+ double lo_l = r_lo - err_d;
+#endif
+ // Ziv's rounding test in the current rounding mode:
+ double upper = r_hi + lo_u;
+ double lower = r_hi + lo_l;
+
+ if (LIBC_LIKELY(upper == lower && hi > -1022))
+ return scale_exp2(upper, hi);
+
+ if (LIBC_UNLIKELY(hi <= -1022))
+ return exp_denorm(x, hi, r_hi, lo_u, err_d, log2_base);
+
+ return exp_accurate_128(x, log2_base);
+}
+
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+} // namespace double_eval
+} // namespace math
} // namespace LIBC_NAMESPACE_DECL
#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_UTILS_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 92d599fba947e..9e9542b891445 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1377,7 +1377,9 @@ add_fp_unittest(
SRCS
exp2_test.cpp
DEPENDS
+ libc.hdr.fenv_macros
libc.src.math.exp2
+ libc.src.__support.math.exp2_double_eval
libc.src.__support.FPUtil.fp_bits
)
diff --git a/libc/test/src/math/exp2_test.cpp b/libc/test/src/math/exp2_test.cpp
index f388a1b57fa4a..ab2ae08af11fd 100644
--- a/libc/test/src/math/exp2_test.cpp
+++ b/libc/test/src/math/exp2_test.cpp
@@ -1,123 +1,188 @@
-//===-- Unittests for 2^x -------------------------------------------------===//
+//===----------------------------------------------------------------------===//
//
// 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 unittests for exp2.
+///
+//===----------------------------------------------------------------------===//
+#include "hdr/fenv_macros.h"
#include "hdr/math_macros.h"
+#include "hdr/stdint_proxy.h"
#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/exp2_double_eval.h"
#include "src/math/exp2.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
#include "utils/MPFRWrapper/MPFRUtils.h"
-#include "hdr/stdint_proxy.h"
-
-using LlvmLibcExp2Test = LIBC_NAMESPACE::testing::FPTest<double>;
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#if defined(LIBC_MATH_HAS_SMALL_TABLES) && \
+ !defined(LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY)
+#define TOLERANCE 2
+#else // LIBC_MATH_HAS_SMALL_TABLE && LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+#define TOLERANCE 1
+#endif // LIBC_MATH_HAS_SMALL_TABLES && !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#define TOLERANCE 0
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+#define ALL_ROUNDING false
+#else
+#define ALL_ROUNDING true
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
using LIBC_NAMESPACE::testing::tlog;
-TEST_F(LlvmLibcExp2Test, SpecialNumbers) {
- EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::exp2(aNaN));
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::exp2(inf));
- EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::exp2(neg_inf));
- EXPECT_FP_EQ_WITH_EXCEPTION(zero, LIBC_NAMESPACE::exp2(-0x1.0p20),
- FE_UNDERFLOW);
- EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::exp2(0x1.0p20), FE_OVERFLOW);
- EXPECT_FP_EQ_ALL_ROUNDING(1.0, LIBC_NAMESPACE::exp2(0.0));
- EXPECT_FP_EQ_ALL_ROUNDING(1.0, LIBC_NAMESPACE::exp2(-0.0));
+LIBC_NOINLINE double exp2_double_eval(double x) {
+ return LIBC_NAMESPACE::math::double_eval::exp2(x);
}
-TEST_F(LlvmLibcExp2Test, TrickyInputs) {
- constexpr int N = 16;
- constexpr uint64_t INPUTS[N] = {
- 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
- 0xc08ff80000000000, // x=-0x1.ff8p+9
- 0xbc971547652b82fe, // x=-0x1.71547652b82fep-54
- };
- for (int i = 0; i < N; ++i) {
- double x = FPBits(INPUTS[i]).get_val();
- EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x,
- LIBC_NAMESPACE::exp2(x), 0.5);
+class Exp2Test : public LIBC_NAMESPACE::testing::FPTest<double> {
+public:
+ void test_special_numbers(double (*func)(double)) {
+ EXPECT_FP_EQ(aNaN, func(aNaN));
+ EXPECT_FP_EQ(inf, func(inf));
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, func(neg_inf));
+ EXPECT_FP_EQ_ALL_ROUNDING(1.0, func(0.0));
+ EXPECT_FP_EQ_ALL_ROUNDING(1.0, func(-0.0));
}
-}
-TEST_F(LlvmLibcExp2Test, 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 = [&](mpfr::RoundingMode rounding_mode) {
- mpfr::ForceRoundingMode __r(rounding_mode);
- if (!__r.success)
- return;
-
- 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 = LIBC_NAMESPACE::exp2(x);
- ++cc;
- if (FPBits(result).is_nan() || FPBits(result).is_inf())
- continue;
-
- ++count;
-
- if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Exp2, x, result,
- 0.5, rounding_mode)) {
- ++fails;
- while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Exp2, x,
- result, tol, rounding_mode)) {
- mx = x;
- mr = result;
-
- if (tol > 1000.0)
- break;
-
- tol *= 2.0;
- }
+ void test_overflow(double (*func)(double)) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(0x1.0p20), FE_OVERFLOW);
+ }
+
+ void test_underflow(double (*func)(double)) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(zero, func(-0x1.0p20), FE_UNDERFLOW);
+ }
+
+ void test_tricky_inputs(double (*func)(double), double ulp_tolerance = 0.5,
+ bool all_rounding = true) {
+ constexpr int N = 15;
+ constexpr uint64_t INPUTS[N] = {
+ 0x3fd79289c6e6a5c0, // x = 0x1.79289c6e6a5cp-2
+ 0x3fd05de80a173ea0, // x = 0x1.05de80a173eap-2
+ 0x3fda7b764e2cf47a, // x = 0x1.a7b764e2cf47ap-2
+ 0x4044c19e5712e377, // x = 0x1.4c19e5712e377p+5
+ 0xc039a74cdab36c28, // x = -0x1.9a74cdab36c28p+4
+ 0xc04757852a4b93aa, // x = -0x1.757852a4b93aap+5
+ 0xc085b3e4e2e3bba9, // x = -0x1.5b3e4e2e3bba9p+9
+ 0xc086960d591aec34, // x = -0x1.6960d591aec34p+9
+ 0xc086232c09d58d91, // x = -0x1.6232c09d58d91p+9
+ 0xc0874910d52d3051, // x = -0x1.74910d52d3051p+9
+ 0xc0867a172ceb0990, // x = -0x1.67a172ceb099p+9
+ 0xc08ff80000000000, // x = -0x1.ff8p+9
+ 0xbc971547652b82fe, // x = -0x1.71547652b82fep-54
+ 0xbf1eb7a4cb841fcc, // x = -0x1.eb7a4cb841fccp-14
+ 0xbf19a61fb925970d, // x = -0x1.9a61fb925970dp-14
+ };
+ for (int i = 0; i < N; ++i) {
+ double x = FPBits(INPUTS[i]).get_val();
+ if (all_rounding) {
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp2, x, func(x),
+ ulp_tolerance);
+ } else {
+ EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, x, func(x), ulp_tolerance);
}
}
- if (fails) {
- tlog << " Exp2 failed: " << fails << "/" << count << "/" << cc
- << " tests.\n";
- tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
- EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, mx, mr, 0.5, rounding_mode);
- }
- };
+ }
+
+ void test_in_range(double (*func)(double), double ulp_tolerance = 0.5,
+ bool all_rounding = true) {
+ 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 = [&](mpfr::RoundingMode rounding_mode) {
+ mpfr::ForceRoundingMode __r(rounding_mode);
+ if (!__r.success)
+ return;
+
+ uint64_t fails = 0;
+ uint64_t count = 0;
+ uint64_t cc = 0;
+ double mx = 0.0, 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 = func(x);
+ ++cc;
+ if (FPBits(result).is_nan() || FPBits(result).is_inf())
+ continue;
+
+ ++count;
+ if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Exp2, x, result,
+ ulp_tolerance, rounding_mode)) {
+ ++fails;
+ while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(
+ mpfr::Operation::Exp2, x, result, tol, rounding_mode)) {
+ mx = x;
+ mr = result;
+
+ if (tol > 1000.0)
+ break;
+
+ tol *= 2.0;
+ }
+ }
+ }
+ if (fails) {
+ tlog << " Exp2 failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol)
+ << ".\n";
+ EXPECT_MPFR_MATCH(mpfr::Operation::Exp2, mx, mr, ulp_tolerance,
+ rounding_mode);
+ }
+ };
- tlog << " Test Rounding To Nearest...\n";
- test(mpfr::RoundingMode::Nearest);
+ tlog << " Test Rounding To Nearest...\n";
+ test(mpfr::RoundingMode::Nearest);
- tlog << " Test Rounding Downward...\n";
- test(mpfr::RoundingMode::Downward);
+ if (all_rounding) {
+ tlog << " Test Rounding Downward...\n";
+ test(mpfr::RoundingMode::Downward);
- tlog << " Test Rounding Upward...\n";
- test(mpfr::RoundingMode::Upward);
+ tlog << " Test Rounding Upward...\n";
+ test(mpfr::RoundingMode::Upward);
- tlog << " Test Rounding Toward Zero...\n";
- test(mpfr::RoundingMode::TowardZero);
-}
+ tlog << " Test Rounding Toward Zero...\n";
+ test(mpfr::RoundingMode::TowardZero);
+ }
+ }
+};
+
+#define LIST_EXP2_TESTS(suffix, func, ulp_tolerance, all_rounding) \
+ using LlvmLibcExp2Test##suffix = Exp2Test; \
+ TEST_F(LlvmLibcExp2Test##suffix, SpecialNumbers) { \
+ test_special_numbers(&func); \
+ } \
+ TEST_F(LlvmLibcExp2Test##suffix, Overflow) { test_overflow(&func); } \
+ TEST_F(LlvmLibcExp2Test##suffix, Underflow) { test_underflow(&func); } \
+ TEST_F(LlvmLibcExp2Test##suffix, TrickyInputs) { \
+ test_tricky_inputs(&func, ulp_tolerance, all_rounding); \
+ } \
+ TEST_F(LlvmLibcExp2Test##suffix, InDoubleRange) { \
+ test_in_range(&func, ulp_tolerance, all_rounding); \
+ } \
+ static_assert(true, "Require semicolon.")
+
+LIST_EXP2_TESTS(Default, LIBC_NAMESPACE::exp2,
+ /*ulp_tolerance=*/TOLERANCE + 0.5,
+ /*all_rounding=*/ALL_ROUNDING);
+LIST_EXP2_TESTS(DoubleEval, exp2_double_eval,
+ /*ulp_tolerance=*/TOLERANCE + 0.5,
+ /*all_rounding=*/ALL_ROUNDING);
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index fe9fafdfab9a3..93caef9f5bf82 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -1480,7 +1480,9 @@ add_fp_unittest(
exp2_test.cpp
DEPENDS
libc.hdr.errno_macros
+ libc.hdr.fenv_macros
libc.src.math.exp2
+ libc.src.__support.math.exp2_double_eval
libc.src.__support.FPUtil.fp_bits
)
diff --git a/libc/test/src/math/smoke/exp2_test.cpp b/libc/test/src/math/smoke/exp2_test.cpp
index 3fd732b790f6b..f968489033571 100644
--- a/libc/test/src/math/smoke/exp2_test.cpp
+++ b/libc/test/src/math/smoke/exp2_test.cpp
@@ -1,69 +1,109 @@
-//===-- Unittests for 2^x -------------------------------------------------===//
+//===----------------------------------------------------------------------===//
//
// 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 exp2.
+///
+//===----------------------------------------------------------------------===//
#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
#include "hdr/math_macros.h"
#include "hdr/stdint_proxy.h"
#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/math/exp2_double_eval.h"
#include "src/math/exp2.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
-using LlvmLibcExp2Test = LIBC_NAMESPACE::testing::FPTest<double>;
+LIBC_NOINLINE double exp2_double_eval(double x) {
+ return LIBC_NAMESPACE::math::double_eval::exp2(x);
+}
-TEST_F(LlvmLibcExp2Test, SpecialNumbers) {
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::exp2(sNaN), FE_INVALID);
- EXPECT_MATH_ERRNO(0);
+class Exp2Test : public LIBC_NAMESPACE::testing::FPTest<double> {
+public:
+ void test_special_numbers(double (*func)(double),
+ bool check_snan_invalid = false) {
+ if (check_snan_invalid) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, func(sNaN), FE_INVALID);
+ } else {
+ EXPECT_FP_EQ(aNaN, func(sNaN));
+ }
+ EXPECT_MATH_ERRNO(0);
- EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::exp2(aNaN));
- EXPECT_MATH_ERRNO(0);
- EXPECT_FP_EQ(inf, LIBC_NAMESPACE::exp2(inf));
- EXPECT_MATH_ERRNO(0);
- EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::exp2(neg_inf));
- EXPECT_MATH_ERRNO(0);
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, func(aNaN));
+ EXPECT_MATH_ERRNO(0);
- EXPECT_FP_EQ_WITH_EXCEPTION(zero, LIBC_NAMESPACE::exp2(-0x1.0p20),
- FE_UNDERFLOW);
- EXPECT_MATH_ERRNO(ERANGE);
- EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::exp2(0x1.0p20), FE_OVERFLOW);
- EXPECT_MATH_ERRNO(ERANGE);
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, func(inf));
+ EXPECT_MATH_ERRNO(0);
- EXPECT_FP_EQ_ALL_ROUNDING(1.0, LIBC_NAMESPACE::exp2(0.0));
- EXPECT_FP_EQ_ALL_ROUNDING(1.0, LIBC_NAMESPACE::exp2(-0.0));
- EXPECT_FP_EQ_ALL_ROUNDING(2.0, LIBC_NAMESPACE::exp2(1.0));
- EXPECT_FP_EQ_ALL_ROUNDING(0.5, LIBC_NAMESPACE::exp2(-1.0));
- EXPECT_FP_EQ_ALL_ROUNDING(4.0, LIBC_NAMESPACE::exp2(2.0));
- EXPECT_FP_EQ_ALL_ROUNDING(0.25, LIBC_NAMESPACE::exp2(-2.0));
-}
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, func(neg_inf));
+ EXPECT_MATH_ERRNO(0);
-#ifdef LIBC_TEST_FTZ_DAZ
+ EXPECT_FP_EQ_ALL_ROUNDING(1.0, func(0.0));
+ EXPECT_MATH_ERRNO(0);
-using namespace LIBC_NAMESPACE::testing;
+ EXPECT_FP_EQ_ALL_ROUNDING(1.0, func(-0.0));
+ EXPECT_MATH_ERRNO(0);
-TEST_F(LlvmLibcExp2Test, FTZMode) {
- ModifyMXCSR mxcsr(FTZ);
+ EXPECT_FP_EQ_ALL_ROUNDING(2.0, func(1.0));
+ EXPECT_FP_EQ_ALL_ROUNDING(0.5, func(-1.0));
+ EXPECT_FP_EQ_ALL_ROUNDING(4.0, func(2.0));
+ EXPECT_FP_EQ_ALL_ROUNDING(0.25, func(-2.0));
+ }
- EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::exp2(min_denormal));
- EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::exp2(max_denormal));
-}
+ void test_overflow(double (*func)(double)) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(0x1.0p20), FE_OVERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+ }
-TEST_F(LlvmLibcExp2Test, DAZMode) {
- ModifyMXCSR mxcsr(DAZ);
+ void test_underflow(double (*func)(double)) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(zero, func(-0x1.0p20), FE_UNDERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+ }
- EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::exp2(min_denormal));
- EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::exp2(max_denormal));
-}
+#ifdef LIBC_TEST_FTZ_DAZ
+ void test_denormals(double (*func)(double)) {
+ EXPECT_FP_EQ(1.0, func(min_denormal));
+ EXPECT_FP_EQ(1.0, func(max_denormal));
+ }
+#endif // LIBC_TEST_FTZ_DAZ
+};
-TEST_F(LlvmLibcExp2Test, FTZDAZMode) {
- ModifyMXCSR mxcsr(FTZ | DAZ);
+#ifdef LIBC_TEST_FTZ_DAZ
+#define LIST_EXP2_FTZ_DAZ_TESTS(suffix, func) \
+ TEST_F(LlvmLibcExp2Test##suffix, FTZMode) { \
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::FTZ); \
+ test_denormals(&func); \
+ } \
+ TEST_F(LlvmLibcExp2Test##suffix, DAZMode) { \
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::DAZ); \
+ test_denormals(&func); \
+ } \
+ TEST_F(LlvmLibcExp2Test##suffix, FTZDAZMode) { \
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::FTZ | \
+ LIBC_NAMESPACE::testing::DAZ); \
+ test_denormals(&func); \
+ }
+#else // !LIBC_TEST_FTZ_DAZ
+#define LIST_EXP2_FTZ_DAZ_TESTS(suffix, func)
+#endif // LIBC_TEST_FTZ_DAZ
- EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::exp2(min_denormal));
- EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::exp2(max_denormal));
-}
+#define LIST_EXP2_TESTS(suffix, func, check_snan_invalid) \
+ using LlvmLibcExp2Test##suffix = Exp2Test; \
+ TEST_F(LlvmLibcExp2Test##suffix, SpecialNumbers) { \
+ test_special_numbers(&func, check_snan_invalid); \
+ } \
+ TEST_F(LlvmLibcExp2Test##suffix, Overflow) { test_overflow(&func); } \
+ TEST_F(LlvmLibcExp2Test##suffix, Underflow) { test_underflow(&func); } \
+ LIST_EXP2_FTZ_DAZ_TESTS(suffix, func) \
+ static_assert(true, "Require semicolon.")
-#endif
+LIST_EXP2_TESTS(Default, LIBC_NAMESPACE::exp2, /*check_snan_invalid=*/true);
+LIST_EXP2_TESTS(DoubleEval, exp2_double_eval, /*check_snan_invalid=*/false);
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index cbe1c92e9e6be..7f3d900173082 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -10350,13 +10350,42 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_exp2_accurate_128",
+ hdrs = ["src/__support/math/exp2_accurate_128.h"],
+ deps = [
+ ":__support_common",
+ ":__support_fputil_dyadic_float",
+ ":__support_fputil_fp_bits",
+ ":__support_frac128",
+ ":__support_frac64",
+ ":__support_macros_attributes",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_uint128",
+ ],
+)
+
libc_support_library(
name = "__support_math_exp_utils",
hdrs = ["src/__support/math/exp_utils.h"],
deps = [
+ ":__support_common",
":__support_cpp_bit",
":__support_cpp_optional",
+ ":__support_fputil_double_double",
+ ":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_rounding_mode",
+ ":__support_macros_attributes",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_macros_properties_architectures",
+ ":__support_macros_properties_cpu_features",
+ ":__support_math_exp2_accurate_128",
+ ":hdr_errno_macros",
+ ":hdr_fenv_macros",
],
)
@@ -10464,30 +10493,51 @@ libc_support_library(
)
libc_support_library(
- name = "__support_math_exp2",
- hdrs = ["src/__support/math/exp2.h"],
+ name = "__support_math_exp_double_utils",
+ hdrs = ["src/__support/math/exp_double_utils.h"],
deps = [
":__support_common",
":__support_cpp_bit",
- ":__support_cpp_optional",
- ":__support_fputil_double_double",
- ":__support_fputil_dyadic_float",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
":__support_fputil_multiply_add",
+ ":__support_macros_attributes",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_math_exp_utils",
+ ":hdr_errno_macros",
+ ],
+)
+
+libc_support_library(
+ name = "__support_math_exp2_double_eval",
+ hdrs = ["src/__support/math/exp2_double_eval.h"],
+ deps = [
+ ":__support_common",
+ ":__support_fputil_double_double",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
":__support_fputil_nearest_integer",
- ":__support_fputil_polyeval",
- ":__support_fputil_rounding_mode",
- ":__support_fputil_triple_double",
- ":__support_integer_literals",
+ ":__support_macros_attributes",
":__support_macros_config",
":__support_macros_optimization",
- ":__support_math_common_constants",
- ":__support_math_exp_constants",
+ ":__support_macros_properties_cpu_features",
+ ":__support_math_exp_double_utils",
":__support_math_exp_utils",
],
)
+libc_support_library(
+ name = "__support_math_exp2",
+ hdrs = ["src/__support/math/exp2.h"],
+ deps = [
+ ":__support_common",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_math_exp2_double_eval",
+ ],
+)
+
libc_support_library(
name = "__support_math_exp2f_double_eval",
hdrs = ["src/__support/math/exp2f_double_eval.h"],
diff --git a/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel
index bb65e582cacd0..1a517f90f585a 100644
--- a/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel
@@ -184,7 +184,12 @@ math_mpfr_test(
],
)
-math_mpfr_test(name = "exp2")
+math_mpfr_test(
+ name = "exp2",
+ deps = [
+ "//libc:__support_math_exp2_double_eval",
+ ],
+)
math_mpfr_test(
name = "exp2f",
diff --git a/utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel
index d7492346e08cf..b19b640506e1e 100644
--- a/utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel
@@ -520,7 +520,12 @@ math_test(
],
)
-math_test(name = "exp2")
+math_test(
+ name = "exp2",
+ deps = [
+ "//libc:__support_math_exp2_double_eval",
+ ],
+)
math_test(
name = "exp2f",
More information about the libc-commits
mailing list