[libc-commits] [libc] [libc][math] Integer-only, statically rounded implementation of exp (PR #225803)
Hoàng Minh Thiên via libc-commits
libc-commits at lists.llvm.org
Sun Oct 25 09:24:48 PDT 2026
https://github.com/hmthien050209 updated https://github.com/llvm/llvm-project/pull/225803
>From 35548113ba69bd3232c4d10c3cccc4802d2e65e5 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Ho=C3=A0ng=20Minh=20Thi=C3=AAn?=
<hoangminhthien05022009 at gmail.com>
Date: Sun, 25 Oct 2026 23:23:24 +0700
Subject: [PATCH] feat: implemented integer-only exp evaluation, complying with
the new structure
---
libc/shared/math/static_rounding/exp.h | 30 +++
libc/shared/static_rounding_math.h | 1 +
libc/src/__support/math/CMakeLists.txt | 33 +++
libc/src/__support/math/exp.h | 8 +
libc/src/__support/math/exp_integer_eval.h | 270 ++++++++++++++++++++
libc/src/__support/math/exp_integer_utils.h | 261 +++++++++++++++++++
libc/src/__support/math/expf_integer_eval.h | 2 +-
libc/test/src/math/exp_test.cpp | 28 +-
libc/test/src/math/smoke/exp_test.cpp | 27 +-
9 files changed, 649 insertions(+), 11 deletions(-)
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/src/__support/math/exp_integer_utils.h
diff --git a/libc/shared/math/static_rounding/exp.h b/libc/shared/math/static_rounding/exp.h
new file mode 100644
index 00000000000000..d09f270c791bad
--- /dev/null
+++ b/libc/shared/math/static_rounding/exp.h
@@ -0,0 +1,30 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 {
+namespace static_rounding {
+
+using math::static_rounding::exp;
+
+} // namespace static_rounding
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_STATIC_ROUNDING_EXPF_H
diff --git a/libc/shared/static_rounding_math.h b/libc/shared/static_rounding_math.h
index 9fc00ade2be962..bd639869aa44af 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/exp2f.h"
#include "shared/math/static_rounding/expf.h"
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 400e5e2d6aeaa8..5cf6ec830002a3 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -4339,11 +4339,44 @@ add_header_library(
exp.h
DEPENDS
.exp_double_eval
+ .exp_integer_eval
libc.src.__support.common
libc.src.__support.macros.config
libc.src.__support.macros.optimization
)
+add_header_library(
+ exp_integer_eval
+ HDRS
+ exp_integer_eval.h
+ DEPENDS
+ .exp_integer_utils
+ libc.hdr.fenv_macros
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.frac64
+ libc.src.__support.frac128
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+)
+
+add_header_library(
+ exp_integer_utils
+ HDRS
+ exp_integer_utils.h
+ DEPENDS
+ libc.hdr.fenv_macros
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.frac64
+ libc.src.__support.frac128
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
+)
+
add_header_library(
exp2f_double_eval
HDRS
diff --git a/libc/src/__support/math/exp.h b/libc/src/__support/math/exp.h
index b2bcb2fd3dedc7..30d1dccbccc2fc 100644
--- a/libc/src/__support/math/exp.h
+++ b/libc/src/__support/math/exp.h
@@ -21,6 +21,14 @@
#include "src/__support/math/exp_double_eval.h"
#define LIBC_MATH_EXP_IMPL double_eval
+#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)
+
+#include "src/__support/math/exp_integer_eval.h"
+#define LIBC_MATH_EXP_IMPL integer_eval
+#endif
+
namespace LIBC_NAMESPACE_DECL {
namespace math {
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 00000000000000..c34d8271cfcc99
--- /dev/null
+++ b/libc/src/__support/math/exp_integer_eval.h
@@ -0,0 +1,270 @@
+//===----------------------------------------------------------------------===//
+//
+// 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
+/// Integer-only implementation of exp.
+///
+//===----------------------------------------------------------------------===//
+
+#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/frac128.h"
+#include "src/__support/frac64.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/exp_integer_utils.h"
+
+// TODO: share *_integer_utils.h with exp2 the same way we did with expf/exp2f.
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace static_rounding {
+
+LIBC_ALWAYS_INLINE uint64_t exp_bits(uint64_t x_u,
+ [[maybe_unused]] int rounding) {
+ using FPBits = fputil::FPBits<double>;
+
+ uint64_t x_abs = x_u & 0x7fff'ffff'ffff'ffff;
+
+ if (x_abs > FPBits::inf().uintval())
+ return x_u | (uint64_t(1) << (FPBits::FRACTION_LEN - 1));
+
+ if (LIBC_UNLIKELY(
+ x_u >= 0xc087'4910'd52d'3052 ||
+ (x_u < 0xbca0'0000'0000'0000 && x_u >= 0x4086'2e42'fefa'39f0) ||
+ x_u <= 0x3ca0'0000'0000'0000)) {
+ // |x| <= 2^-53
+ if (x_abs <= 0x3ca0'0000'0000'0000) {
+ uint64_t one = FPBits::one().uintval();
+ if (x_abs == 0)
+ return one;
+ bool is_neg = (x_u >> 63) != 0;
+ if (x_abs == 0x3ca0'0000'0000'0000) {
+ if (rounding == FE_TONEAREST || (rounding == FE_UPWARD && !is_neg) ||
+ ((rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) && is_neg))
+ return is_neg ? one - 1 : one + 1;
+ }
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (rounding == FE_UPWARD && !is_neg)
+ return one + 1;
+ if ((rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO) && is_neg)
+ return one - 1;
+#endif
+ return one;
+ }
+
+ // x <= log(2^-1075) || x >= 0x1.6232bdd7abcd3p+9 or inf/nan.
+
+ // x <= log(2^-1075) or -inf/nan
+ if (x_u >= 0xc087'4910'd52d'3052) {
+ // exp(-Inf) = 0
+ if (x_abs == FPBits::inf().uintval())
+ return FPBits::zero().uintval();
+
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (rounding == FE_UPWARD)
+ return FPBits::min_subnormal().uintval();
+#endif
+ return FPBits::zero().uintval();
+ }
+
+ // x >= round(log(MAX_NORMAL), D, RU) = 0x1.62e42fefa39fp+9 or +inf/nan
+ // x is finite
+ if (x_u < 0x7ff0'0000'0000'0000) {
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO)
+ return FPBits::max_normal().uintval();
+#endif
+ }
+
+ // x is +inf or nan
+ return FPBits::inf().uintval();
+ }
+
+ // Main calculations
+
+ bool is_neg = (x_u >> 63) != 0;
+ uint32_t x_e = static_cast<uint32_t>(x_abs >> FPBits::FRACTION_LEN);
+ uint64_t x_s = x_abs & FPBits::FRACTION_MASK;
+
+ // The idea of the algorithm below is that:
+ // Let:
+ // - y = x / ln(2) --> exp(x) = 2^y.
+ // - y = hi + mid + low
+ //
+ // With:
+ // - hi is an integer
+ // - mid * 2^4 is an integer
+ // - lo are the remainder bits
+ //
+ // Then:
+ // exp(x) = 2^hi * 2^mid * 2^lo
+ //
+ // With this formula:
+ // - Multiplying by 2^hi is exact and cheap, via adding into the exponent
+ // field
+ // - 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.
+
+ // Range reduction
+
+ // binary64 recall:
+ // 1 sign bit
+ // 11 exp bits
+ // 52 sig bits
+
+ uint64_t x_s_shifted = (x_s << 11) | (1ull << 63);
+ // LSB(x_s_frac) = 2^-52
+ Frac64 x_s_frac(x_s_shifted);
+
+ Frac64 x_ln2 = x_s_frac * LOG2_E_F64;
+ uint64_t x_ln2_bits = x_ln2.val[0];
+
+ uint64_t k = 0;
+ uint64_t frac_bits = 0;
+ int shift = 0;
+
+ int x_e_unbiased = static_cast<int>(x_e) - FPBits::EXP_BIAS;
+ if (x_e_unbiased >= -1) {
+ shift = 62 - x_e_unbiased;
+ k = x_ln2_bits >> shift;
+ frac_bits = x_ln2_bits << (64 - shift);
+ } else {
+ k = 0;
+ shift = -x_e_unbiased - 2;
+ frac_bits = (shift < 64) ? (x_ln2_bits >> shift) : 0;
+ }
+
+ // As Frac64/Frac128 can't store the sign, we need to handle the sign
+ // separately:
+ // - Both branches are computing floor(x * log2(e)).
+ // - For negative x, we round up to the next multiple of 2^52, then clear
+ // the last 52 bits.
+ // - For positive x, we round down (just clear) the last 52 bits.
+ //
+ // Then, l2y_r_hi is the remainder of x * log2(e) after removing the
+ // integer part, which is used to look up EXP_MID and compute 2^lo.
+ //
+ // e_y is the biased exponent field, positioned for the final bit assembly.
+ uint64_t l2y_r_hi;
+ uint64_t e_y;
+
+ if (LIBC_UNLIKELY(is_neg)) {
+ if (frac_bits != 0) {
+ k = k + 1;
+ l2y_r_hi = ~frac_bits + 1; // 1 - r
+ } else {
+ l2y_r_hi = 0;
+ }
+ e_y = (FPBits::EXP_BIAS << 20) - static_cast<uint32_t>(k << 20);
+ } else {
+ l2y_r_hi = frac_bits;
+ e_y = (FPBits::EXP_BIAS << 20) + static_cast<uint32_t>(k << 20);
+ }
+
+ 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 top 4 fractional bits for the LUT index.
+ uint16_t x_mid = static_cast<uint16_t>((l2y_r_hi >> 60) & 0xf);
+
+ // The remaining 60 bits are the polynomial input.
+ uint64_t x_lo = l2y_r_hi & ((1ull << 60) - 1);
+
+ // Fast path: 64-bit calculations first
+
+ // Don't << 4, as the polynomial approximation is correct in range [0, 1/16],
+ // and x_lo is already in that range.
+ // LSB(x_lo_frac) = 2^-64
+ Frac64 x_lo_frac(x_lo);
+
+ Frac64 p = x_lo_frac * fputil::polyeval(x_lo_frac, EXP_64_COEFFS[0],
+ EXP_64_COEFFS[1], EXP_64_COEFFS[2],
+ EXP_64_COEFFS[3], EXP_64_COEFFS[4],
+ EXP_64_COEFFS[5]);
+
+ // With:
+ // - p = 2^lo - 1 --> 2^lo = p + 1
+ // - 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_val = EXP_MID[x_mid].to_frac64();
+ Frac64 result = fputil::multiply_add(mid_val, p, mid_val);
+
+ uint64_t result_bits = result.val[0] << 1;
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return exp_handle_rounding(result, is_neg, d, e_y, rounding);
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+ // Rounding test
+ constexpr uint32_t LAST_BITS = 12;
+ constexpr uint32_t ROUNDING_ERROR = (1ull << 11);
+ uint32_t result_last_bits =
+ static_cast<uint32_t>(result_bits & ((1u << LAST_BITS) - 1));
+ bool is_hard;
+ if (rounding == FE_TONEAREST) {
+ uint32_t rounded_lo =
+ (result_last_bits + (1u << (LAST_BITS - 1)) - ROUNDING_ERROR) >>
+ LAST_BITS;
+ uint32_t rounded_hi =
+ (result_last_bits + (1u << (LAST_BITS - 1)) + ROUNDING_ERROR) >>
+ LAST_BITS;
+ is_hard = rounded_lo != rounded_hi;
+ } else {
+ is_hard = result_last_bits <= ROUNDING_ERROR ||
+ result_last_bits >= (1u << LAST_BITS) - ROUNDING_ERROR;
+ }
+
+ if (LIBC_LIKELY(!is_hard)) {
+ return exp_handle_rounding(result, is_neg, d, e_y, rounding);
+ }
+
+ // Dial back to 128-bit path for hard-to-round cases
+ return exp_accurate_path(x_s_shifted, x_e_unbiased, is_neg, rounding);
+}
+
+LIBC_ALWAYS_INLINE double exp(double x, [[maybe_unused]] int rounding) {
+ return cpp::bit_cast<double>(exp_bits(cpp::bit_cast<uint64_t>(x), rounding));
+}
+
+} // namespace static_rounding
+
+namespace integer_eval {
+
+LIBC_ALWAYS_INLINE double exp(double x) {
+ return 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/__support/math/exp_integer_utils.h b/libc/src/__support/math/exp_integer_utils.h
new file mode 100644
index 00000000000000..bcb1e0de4ac5ce
--- /dev/null
+++ b/libc/src/__support/math/exp_integer_utils.h
@@ -0,0 +1,261 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Shared utilities and constants for integer-only double-precision exponential
+/// functions.
+///
+//===----------------------------------------------------------------------===//
+//
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/frac128.h"
+#include "src/__support/frac64.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/types.h"
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_UTILS_H
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+namespace static_rounding {
+
+// All constants below are generated via Sollya
+
+// > display = hexadecimal;
+// > print(2+round(1/log(2), 128, RN));
+// LSB(LOG2_E_F128) = 2^-127
+LIBC_INLINE_VAR constexpr Frac128 LOG2_E_F128 =
+ Frac128({0xbe87'fed0'691d'3e89, 0xb8aa'3b29'5c17'f0bb});
+
+// > display = hexadecimal;
+// > print(2+round(1/log(2), 64, RN));
+// LSB(LOG2_E_F64) = 2^-63
+LIBC_INLINE_VAR constexpr Frac64 LOG2_E_F64 = Frac64(0xb8aa'3b29'5c17'f0bc);
+
+// 2^x for x from 0 to (0b1111/2^4) = 15/16
+// > for i from 0 to 15 do {
+// print(1+round(2^(i/16), 128, RN));
+// };
+// LSB(EXP_MID4[i]) = 2^-127
+LIBC_INLINE_VAR constexpr Frac128 EXP_MID[] = {
+ Frac128({0x0000'0000'0000'0000ULL, 0x8000'0000'0000'0000ULL}),
+ Frac128({0xc5c9'5b8c'2154'c1b2ULL, 0x85aa'c367'cc48'7b14ULL}),
+ Frac128({0xfbe4'6287'58a5'3c90ULL, 0x8b95'c1e3'ea8b'd6e6ULL}),
+ Frac128({0x0fd6'd8e0'ae5a'c9d8ULL, 0x91c3'd373'ab11'c336ULL}),
+ Frac128({0x46ad'2318'2e42'f6f6ULL, 0x9837'f051'8db8'a96fULL}),
+ Frac128({0xa091'1f09'ebb9'fdd1ULL, 0x9ef5'3260'91a1'11adULL}),
+ Frac128({0x1cbd'7f62'1710'701bULL, 0xa5fe'd6a9'b151'38eaULL}),
+ Frac128({0x4980'a8c8'f59a'2ec4ULL, 0xad58'3eea'42a1'4ac6ULL}),
+ Frac128({0x597d'89b3'754a'be9fULL, 0xb504'f333'f9de'6484ULL}),
+ Frac128({0xa881'1fb6'6d0f'af7aULL, 0xbd08'a39f'580c'36beULL}),
+ Frac128({0x3e2a'd0c9'64dd'9f37ULL, 0xc567'2a11'5506'daddULL}),
+ Frac128({0xe235'838f'95f2'c6edULL, 0xce24'8c15'1f84'80e3ULL}),
+ Frac128({0x39a6'8bb9'902d'3fdeULL, 0xd744'fcca'd69d'6af4ULL}),
+ Frac128({0x0658'9504'8dd3'33caULL, 0xe0cc'deec'2a94'e111ULL}),
+ Frac128({0xd02d'75b3'706e'54fbULL, 0xeac0'c6e7'dd24'392eULL}),
+ Frac128({0x7b9d'0c7a'ed98'0fc3ULL, 0xf525'7d15'2486'cc2cULL}),
+};
+
+// 128-bit polynomial approximation of 2^x coefficients generated with Sollya:
+// > P = fpminimax(2^x, 12, [|1, 128...|], [0, 1/16], absolute, fixed);
+// Store the fractional part of the coefficients below
+// > dirtyinfnorm(2^x - P(x), [0, 1/16]);
+// 0x1.b328...p-117
+// LSB(EXPF_COEFFS[i]) = 2^-128
+LIBC_INLINE_VAR constexpr Frac128 EXP_COEFFS[] = {
+ // degree-0 = 1, add back afterwards to reduce calc ops
+ Frac128({0xc9e3'b398'033f'0902ULL, 0xb172'17f7'd1cf'79abULL}),
+ Frac128({0xde2d'60e0'866c'6365ULL, 0x3d7f'7bff'058b'1d50ULL}),
+ Frac128({0x99d3'ad04'd47d'0efeULL, 0x0e35'846b'8250'5fc5ULL}),
+ Frac128({0x399a'b423'1644'1b55ULL, 0x0276'556d'f749'cee5ULL}),
+ Frac128({0x405f'ebef'0f76'011aULL, 0x0057'61ff'9e29'9cc4ULL}),
+ Frac128({0x1ac4'9999'c0ef'b9abULL, 0x000a'1848'97c3'63c4ULL}),
+ Frac128({0xe6f8'52f9'914f'2d9aULL, 0x0000'ffe5'fe2c'4573ULL}),
+ Frac128({0x6056'28a6'5061'b43fULL, 0x0000'162c'0223'a816ULL}),
+ Frac128({0x6b99'3efe'2193'e54cULL, 0x0000'01b5'253d'0671ULL}),
+ Frac128({0x9c76'f49c'ed65'743aULL, 0x0000'001e'4cf8'0b70ULL}),
+ Frac128({0x1ac1'3330'6c4e'3d33ULL, 0x0000'0001'e8ae'6938ULL}),
+ Frac128({0x219f'9904'1f29'5f13ULL, 0x0000'0000'1cd9'af73ULL}),
+};
+
+// 64-bit polynomial approximation of 2^x coefficients generated with Sollya:
+// > P = fpminimax(2^x, 6, [|1, 64...|], [0, 1/16], absolute, fixed);
+// Store the fractional part of the coefficients below
+// > dirtyinfnorm(2^x - P(x), [0, 1/16]);
+// 0x1.2ac9...p-57
+// This is different from EXPF_COEFFS: EXPF_COEFFS is for approximating for x in
+// range of [0, 1]
+// LSB(EXP_COEFFS[i]) = 2^-64
+LIBC_INLINE_VAR constexpr Frac64 EXP_64_COEFFS[] = {
+ // degree-0 = 1, add back afterwards to reduce calc ops
+ Frac64(0xb172'17f7'd1cd'3e3cULL), Frac64(0x3d7f'7bff'0838'4d4aULL),
+ Frac64(0x0e35'846a'6f46'd462ULL), Frac64(0x0276'55a0'd536'1dd9ULL),
+ Frac64(0x0057'5d3d'4b45'0056ULL), Frac64(0x000a'504b'13fe'e008ULL),
+};
+
+// Round the fractional result and combine it with its exponent.
+template <typename TFrac, typename TUInt>
+LIBC_INLINE uint64_t exp_handle_rounding(TFrac result_frac, bool is_neg, int d,
+ TUInt e_y,
+ [[maybe_unused]] int rounding) {
+ static_assert(cpp::is_same<TFrac, Frac64>::value ||
+ cpp::is_same<TFrac, Frac128>::value);
+ constexpr bool IS_FAST_PATH = cpp::is_same<TFrac, Frac64>::value;
+
+ uint32_t shift_length = 11;
+ uint64_t leading_one = 0;
+
+ // subnormal
+ if (LIBC_UNLIKELY(is_neg && d >= 0)) {
+ e_y = 0;
+ leading_one = uint64_t(1) << (52 - d);
+
+ // Truncate the last 2 bits to avoid undefined behavior when shifting by 64
+ // bits
+ if (d >= 51) {
+ d -= 2;
+ if constexpr (IS_FAST_PATH)
+ result_frac.val[0] >>= 2;
+ else
+ result_frac.val[1] >>= 2;
+ }
+
+ shift_length += d + 1;
+ }
+
+ // Get the bits, discarding the leading bit
+ auto frac_bits = [&]() -> uint64_t {
+ if constexpr (IS_FAST_PATH)
+ return result_frac.val[0] << 1;
+ else
+ return (result_frac.val[1] << 1) | (result_frac.val[0] >> 63);
+ };
+
+#ifdef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ TUInt result = static_cast<TUInt>(frac_bits() >> shift_length) +
+ static_cast<TUInt>(leading_one);
+ bool round_bit = (result & 1) != 0;
+ result >>= 1;
+ bool sticky = (frac_bits() & ((uint64_t(1) << shift_length) - 1)) != 0;
+ if constexpr (!IS_FAST_PATH)
+ sticky = sticky || (result_frac.val[0] != 0);
+ if (round_bit && (sticky || (result & 1)))
+ ++result;
+ result += static_cast<TUInt>(e_y) << 32;
+
+ return result;
+#else // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+ if (rounding == FE_TONEAREST) {
+ TUInt result = static_cast<TUInt>(frac_bits() >> shift_length) +
+ static_cast<TUInt>(leading_one);
+ bool round_bit = (result & 1) != 0;
+ result >>= 1;
+ bool sticky = (frac_bits() & ((uint64_t(1) << shift_length) - 1)) != 0;
+ if constexpr (!IS_FAST_PATH)
+ sticky = sticky || (result_frac.val[0] != 0);
+ if (round_bit && (sticky || (result & 1)))
+ ++result;
+ result += static_cast<TUInt>(e_y) << 32;
+
+ return result;
+ }
+
+ TUInt should_round_up = 0;
+
+ if (LIBC_UNLIKELY(rounding == FE_UPWARD)) {
+ uint64_t round_up_mask = (uint64_t(1) << (shift_length + 1)) - 1;
+ bool has_remainder = (frac_bits() & round_up_mask) != 0;
+ if constexpr (!IS_FAST_PATH) {
+ has_remainder = has_remainder || (result_frac.val[0] != 0);
+ }
+ should_round_up = static_cast<TUInt>(has_remainder);
+ }
+
+ TUInt result = (static_cast<TUInt>(frac_bits() >> (shift_length + 1)) +
+ should_round_up + (leading_one >> 1));
+ result += static_cast<TUInt>(e_y) << 32;
+
+ return result;
+#endif // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+}
+
+LIBC_ALWAYS_INLINE uint64_t exp_accurate_path(uint64_t x_s_shifted,
+ int x_e_unbiased, bool is_neg,
+ int rounding) {
+ using FPBits = fputil::FPBits<double>;
+
+ // Recalculate everything in 128-bit precision, with the same idea as the
+ // 64-bit path.
+ //
+ // TODO: any way not to recompute?
+
+ Frac128 x_s_frac({0, x_s_shifted});
+ Frac128 x_ln2 = x_s_frac * LOG2_E_F128;
+
+ uint64_t k = 0;
+ Frac128 l2y_r;
+ if (x_e_unbiased >= -1) {
+ int shift = 62 - x_e_unbiased;
+ k = (x_ln2 >> (64 + shift)).val[0];
+ l2y_r = x_ln2 << (64 - shift);
+ } else {
+ int shift = -x_e_unbiased - 2;
+ l2y_r = (shift < 128) ? (x_ln2 >> shift) : Frac128{};
+ }
+
+ if (LIBC_UNLIKELY(is_neg)) {
+ if (l2y_r.val[0] != 0 || l2y_r.val[1] != 0) {
+ ++k;
+ l2y_r = ~l2y_r + Frac128(1);
+ }
+ }
+
+ uint64_t e_y;
+ if (is_neg)
+ e_y = (FPBits::EXP_BIAS << 20) - static_cast<uint32_t>(k << 20);
+ else
+ e_y = (FPBits::EXP_BIAS << 20) + static_cast<uint32_t>(k << 20);
+
+ int d = static_cast<int>(k) - FPBits::EXP_BIAS;
+
+ if (LIBC_UNLIKELY(is_neg && d >= 53))
+ return 0.0;
+
+ uint16_t x_mid = static_cast<uint16_t>((l2y_r.val[1] >> 60) & 0xf);
+
+ Frac128 x_lo_frac = l2y_r;
+ x_lo_frac.val[1] &= 0x0fff'ffff'ffff'ffffULL;
+
+ Frac128 p =
+ x_lo_frac * fputil::polyeval(x_lo_frac, EXP_COEFFS[0], EXP_COEFFS[1],
+ EXP_COEFFS[2], EXP_COEFFS[3], EXP_COEFFS[4],
+ EXP_COEFFS[5], EXP_COEFFS[6], EXP_COEFFS[7],
+ EXP_COEFFS[8], EXP_COEFFS[9], EXP_COEFFS[10],
+ EXP_COEFFS[11]);
+
+ Frac128 mid_val = EXP_MID[x_mid];
+ Frac128 result = fputil::multiply_add(mid_val, p, mid_val);
+
+ return exp_handle_rounding(result, is_neg, d, e_y, rounding);
+}
+
+} // namespace static_rounding
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXP_INTEGER_UTILS_H
diff --git a/libc/src/__support/math/expf_integer_eval.h b/libc/src/__support/math/expf_integer_eval.h
index 89b61870c4e98f..9d9c687a0473bf 100644
--- a/libc/src/__support/math/expf_integer_eval.h
+++ b/libc/src/__support/math/expf_integer_eval.h
@@ -30,7 +30,7 @@ namespace static_rounding {
// Constants generated by Sollya with:
// > display = hexadecimal;
-// > round(log2(exp(1)) * 2^63, 64, RN);
+// > print(2+round(1/log(2), 64, RN));
// LSB(LOG2_E) = 2^-63
LIBC_INLINE_VAR constexpr Frac64 LOG2_E = Frac64(0xb8aa'3b29'5c17'f0bc);
diff --git a/libc/test/src/math/exp_test.cpp b/libc/test/src/math/exp_test.cpp
index dd764aa5b27dc0..44d3616cdba71a 100644
--- a/libc/test/src/math/exp_test.cpp
+++ b/libc/test/src/math/exp_test.cpp
@@ -18,6 +18,7 @@
#include "src/__support/macros/attributes.h"
#include "src/__support/macros/optimization.h"
#include "src/__support/math/exp_double_eval.h"
+#include "src/__support/math/exp_integer_eval.h"
#include "src/math/exp.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
@@ -168,13 +169,18 @@ class ExpTest : public LIBC_NAMESPACE::testing::FPTest<double> {
}
};
-#define LIST_EXP_TESTS(suffix, func, ulp_tolerance, all_rounding) \
+#define LIST_EXP_TESTS(suffix, func, ulp_tolerance, all_rounding, \
+ check_exception) \
using LlvmLibcExpTest##suffix = ExpTest; \
TEST_F(LlvmLibcExpTest##suffix, SpecialNumbers) { \
test_special_numbers(&func); \
} \
- TEST_F(LlvmLibcExpTest##suffix, Overflow) { test_overflow(&func); } \
- TEST_F(LlvmLibcExpTest##suffix, Underflow) { test_underflow(&func); } \
+ TEST_F(LlvmLibcExpTest##suffix, OverflowUnderflow) { \
+ if (check_exception) { \
+ test_overflow(&func); \
+ test_underflow(&func); \
+ } \
+ } \
TEST_F(LlvmLibcExpTest##suffix, TrickyInputs) { \
test_tricky_inputs(&func, ulp_tolerance, all_rounding); \
} \
@@ -185,7 +191,19 @@ class ExpTest : public LIBC_NAMESPACE::testing::FPTest<double> {
LIST_EXP_TESTS(Default, LIBC_NAMESPACE::exp,
/*ulp_tolerance=*/TOLERANCE + 0.5,
- /*all_rounding=*/ALL_ROUNDING);
+ /*all_rounding=*/ALL_ROUNDING, /*check_exception=*/true);
LIST_EXP_TESTS(DoubleEval, exp_double_eval,
/*ulp_tolerance=*/TOLERANCE + 0.5,
- /*all_rounding=*/ALL_ROUNDING);
+ /*all_rounding=*/ALL_ROUNDING, /*check_exception=*/true);
+LIST_EXP_TESTS(IntegerEval, LIBC_NAMESPACE::math::integer_eval::exp,
+ /*ulp_tolerance=*/1.0, /*all_rounding=*/false,
+ /*check_exception=*/false);
+
+static double exp_static_rounding(double x) {
+ return LIBC_NAMESPACE::math::static_rounding::exp(
+ x, LIBC_NAMESPACE::fputil::quick_get_round());
+}
+
+LIST_EXP_TESTS(StaticRounding, exp_static_rounding,
+ /*ulp_tolerance=*/1.0, /*all_rounding=*/true,
+ /*check_exception=*/false);
diff --git a/libc/test/src/math/smoke/exp_test.cpp b/libc/test/src/math/smoke/exp_test.cpp
index 608d1084ab692b..060fca864e2206 100644
--- a/libc/test/src/math/smoke/exp_test.cpp
+++ b/libc/test/src/math/smoke/exp_test.cpp
@@ -18,6 +18,7 @@
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/macros/attributes.h"
#include "src/__support/math/exp_double_eval.h"
+#include "src/__support/math/exp_integer_eval.h"
#include "src/math/exp.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
@@ -90,15 +91,31 @@ class ExpTest : public LIBC_NAMESPACE::testing::FPTest<double> {
#define LIST_EXP_FTZ_DAZ_TESTS(suffix, func)
#endif // LIBC_TEST_FTZ_DAZ
-#define LIST_EXP_TESTS(suffix, func, check_snan_invalid) \
+#define LIST_EXP_TESTS(suffix, func, check_snan_invalid, check_exception) \
using LlvmLibcExpTest##suffix = ExpTest; \
TEST_F(LlvmLibcExpTest##suffix, SpecialNumbers) { \
test_special_numbers(&func, check_snan_invalid); \
} \
- TEST_F(LlvmLibcExpTest##suffix, Overflow) { test_overflow(&func); } \
- TEST_F(LlvmLibcExpTest##suffix, Underflow) { test_underflow(&func); } \
+ TEST_F(LlvmLibcExpTest##suffix, OverflowUnderflow) { \
+ if (check_exception) { \
+ test_overflow(&func); \
+ test_underflow(&func); \
+ } \
+ } \
LIST_EXP_FTZ_DAZ_TESTS(suffix, func) \
static_assert(true, "Require semicolon.")
-LIST_EXP_TESTS(Default, LIBC_NAMESPACE::exp, /*check_snan_invalid=*/true);
-LIST_EXP_TESTS(DoubleEval, exp_double_eval, /*check_snan_invalid=*/false);
+LIST_EXP_TESTS(Default, LIBC_NAMESPACE::exp, /*check_snan_invalid=*/true,
+ /*check_exception=*/true);
+LIST_EXP_TESTS(DoubleEval, exp_double_eval, /*check_snan_invalid=*/false,
+ /*check_exception=*/true);
+LIST_EXP_TESTS(IntegerEval, LIBC_NAMESPACE::math::integer_eval::exp,
+ /*check_snan_invalid=*/false, /*check_exception=*/false);
+
+static double exp_static_rounding(double x) {
+ return LIBC_NAMESPACE::math::static_rounding::exp(
+ x, LIBC_NAMESPACE::fputil::quick_get_round());
+}
+
+LIST_EXP_TESTS(StaticRounding, exp_static_rounding,
+ /*check_snan_invalid=*/false, /*check_exception=*/false);
More information about the libc-commits
mailing list