[libc-commits] [libc] [llvm] [libc][math] Implement double precision sinh correctly rounded for all rounding modes. (PR #235097)
via libc-commits
libc-commits at lists.llvm.org
Fri Oct 23 13:58:54 PDT 2026
llvmorg-github-actions[bot] wrote:
<!--LLVM PR SUMMARY COMMENT-->
@llvm/pr-subscribers-backend-amdgpu
Author: lntue
<details>
<summary>Changes</summary>
- Extend even/odd polynomial evalutions of exp(x) +- exp(-x) from single precision to double-double and 128-bit fixed point.
- Optimize range reduction using the trict from previous exp implementation, combining with the new fp-bitwise operations, and reapply it to exp and exp2.
Assisted-by: Gemini for analysis and refactoring.
---
Patch is 85.30 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/235097.diff
36 Files Affected:
- (modified) libc/config/baremetal/aarch64/entrypoints.txt (+1)
- (modified) libc/config/baremetal/arm/entrypoints.txt (+1)
- (modified) libc/config/baremetal/riscv/entrypoints.txt (+1)
- (modified) libc/config/darwin/aarch64/entrypoints.txt (+1)
- (modified) libc/config/freebsd/x86_64/entrypoints.txt (+1)
- (modified) libc/config/gpu/amdgpu/entrypoints.txt (+1)
- (modified) libc/config/gpu/nvptx/entrypoints.txt (+1)
- (modified) libc/config/linux/aarch64/entrypoints.txt (+1)
- (modified) libc/config/linux/arm/entrypoints.txt (+1)
- (modified) libc/config/linux/riscv/entrypoints.txt (+1)
- (modified) libc/config/linux/x86_64/entrypoints.txt (+1)
- (modified) libc/config/windows/entrypoints.txt (+1)
- (modified) libc/docs/headers/math/index.md (+1-1)
- (modified) libc/include/math.yaml (+6)
- (modified) libc/shared/math.h (+1)
- (added) libc/shared/math/sinh.h (+28)
- (modified) libc/src/__support/FPUtil/double_double.h (+18)
- (modified) libc/src/__support/math/CMakeLists.txt (+73-2)
- (modified) libc/src/__support/math/exp2_double_eval.h (+15-8)
- (modified) libc/src/__support/math/exp_double_eval.h (+17-8)
- (modified) libc/src/__support/math/exp_utils.h (+65)
- (added) libc/src/__support/math/sinh.h (+34)
- (added) libc/src/__support/math/sinh_double_eval.h (+237)
- (added) libc/src/__support/math/sinhcosh_accurate_128.h (+309)
- (added) libc/src/__support/math/sinhcosh_utils.h (+467)
- (modified) libc/src/math/generic/CMakeLists.txt (+10)
- (added) libc/src/math/generic/sinh.cpp (+21)
- (modified) libc/test/shared/CMakeLists.txt (+1)
- (modified) libc/test/shared/shared_math_test.cpp (+1)
- (modified) libc/test/src/math/CMakeLists.txt (+16)
- (added) libc/test/src/math/sinh_test.cpp (+187)
- (modified) libc/test/src/math/smoke/CMakeLists.txt (+15)
- (added) libc/test/src/math/smoke/sinh_test.cpp (+123)
- (modified) utils/bazel/llvm-project-overlay/libc/BUILD.bazel (+78-2)
- (modified) utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel (+8)
- (modified) utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel (+8)
``````````diff
diff --git a/libc/config/baremetal/aarch64/entrypoints.txt b/libc/config/baremetal/aarch64/entrypoints.txt
index 7ec212d8709692..db6033e9a2c04c 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/entrypoints.txt
@@ -596,6 +596,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sinpif
libc.src.math.sqrt
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index 779d9fc23b28e8..e4964aa02a6ec5 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -608,6 +608,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sinpif
libc.src.math.sqrt
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index b847a2bb1f47f0..d9baaaeba8a10e 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -599,6 +599,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sinpif
libc.src.math.sqrt
diff --git a/libc/config/darwin/aarch64/entrypoints.txt b/libc/config/darwin/aarch64/entrypoints.txt
index 5600b82ad3410b..cabd03e2442dfb 100644
--- a/libc/config/darwin/aarch64/entrypoints.txt
+++ b/libc/config/darwin/aarch64/entrypoints.txt
@@ -417,6 +417,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sinpif
libc.src.math.sqrt
diff --git a/libc/config/freebsd/x86_64/entrypoints.txt b/libc/config/freebsd/x86_64/entrypoints.txt
index 8f6c178301c2f2..6d567b01571e84 100644
--- a/libc/config/freebsd/x86_64/entrypoints.txt
+++ b/libc/config/freebsd/x86_64/entrypoints.txt
@@ -305,6 +305,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincosf
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sqrt
libc.src.math.sqrtf
diff --git a/libc/config/gpu/amdgpu/entrypoints.txt b/libc/config/gpu/amdgpu/entrypoints.txt
index 2eaef608e9a7ac..196b6a3fd20266 100644
--- a/libc/config/gpu/amdgpu/entrypoints.txt
+++ b/libc/config/gpu/amdgpu/entrypoints.txt
@@ -534,6 +534,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sinpif
libc.src.math.sqrt
diff --git a/libc/config/gpu/nvptx/entrypoints.txt b/libc/config/gpu/nvptx/entrypoints.txt
index 79d0c733ec3c01..96305e3d1293f1 100644
--- a/libc/config/gpu/nvptx/entrypoints.txt
+++ b/libc/config/gpu/nvptx/entrypoints.txt
@@ -531,6 +531,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sinpif
libc.src.math.sqrt
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index ae4d060ad2e65c..f78eb4d3c8a372 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -834,6 +834,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sinpif
libc.src.math.sqrt
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index e35ca9c7211ab3..3dd2eede1952d4 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -611,6 +611,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sqrt
libc.src.math.sqrtf
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index f8906e3adbfd30..689eeb8f6e3d69 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -895,6 +895,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sinpif
libc.src.math.sqrt
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5d67548617ea5b..a42172e2cce6fa 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -917,6 +917,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincos
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sinpif
libc.src.math.sqrt
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index 1712b2ad26fda4..f38000fcddb7ac 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -338,6 +338,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sincosf
libc.src.math.sincosf
libc.src.math.sinf
+ libc.src.math.sinh
libc.src.math.sinhf
libc.src.math.sqrt
libc.src.math.sqrtf
diff --git a/libc/docs/headers/math/index.md b/libc/docs/headers/math/index.md
index 30aba90c60a9bc..a4596408ee5648 100644
--- a/libc/docs/headers/math/index.md
+++ b/libc/docs/headers/math/index.md
@@ -211,7 +211,7 @@ log
| rsqrt | {{ check }} | | | {{ check }} | | | 7.12.7.9 | F.10.4.9 |
| sin | {{ check }} | {{ check }} | | {{ check }} | | | 7.12.4.6 | F.10.1.6 |
| sincos | {{ check }} | {{ check }} | | | | | | |
-| sinh | {{ check }} | | | {{ check }} | | | 7.12.5.5 | F.10.2.5 |
+| sinh | {{ check }} | {{ check }} | | {{ check }} | | | 7.12.5.5 | F.10.2.5 |
| sinpi | {{ check }} | | | {{ check }} | | | 7.12.4.13 | F.10.1.13 |
| sqrt | {{ check }} | {{ check }} | {{ check }} | {{ check }} | {{ check }} | {{ check }} | 7.12.7.10 | F.10.4.10 |
| tan | {{ check }} | {{ check }} | | {{ check }} | | {{ check }} | 7.12.4.7 | F.10.1.7 |
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index 006233e579f86c..89e0ab4f0dac60 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -2653,6 +2653,12 @@ functions:
- type: float
- type: float *
- type: float *
+ - name: sinh
+ standards:
+ - stdc
+ return_type: double
+ arguments:
+ - type: double
- name: sinhf
standards:
- stdc
diff --git a/libc/shared/math.h b/libc/shared/math.h
index eb72bebcb53bfa..04f85e000bd0b6 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -467,6 +467,7 @@
#include "math/sincosf.h"
#include "math/sinf.h"
#include "math/sinf16.h"
+#include "math/sinh.h"
#include "math/sinhf.h"
#include "math/sinhf16.h"
#include "math/sinpif.h"
diff --git a/libc/shared/math/sinh.h b/libc/shared/math/sinh.h
new file mode 100644
index 00000000000000..8e7859fdcf665f
--- /dev/null
+++ b/libc/shared/math/sinh.h
@@ -0,0 +1,28 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+//
+//===----------------------------------------------------------------------===//
+///
+/// \file
+/// Shared sinh function.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SHARED_MATH_SINH_H
+#define LLVM_LIBC_SHARED_MATH_SINH_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/sinh.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::sinh;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_SINH_H
diff --git a/libc/src/__support/FPUtil/double_double.h b/libc/src/__support/FPUtil/double_double.h
index 6e41bb9b2a4bed..4d93c4769b1fb2 100644
--- a/libc/src/__support/FPUtil/double_double.h
+++ b/libc/src/__support/FPUtil/double_double.h
@@ -67,6 +67,24 @@ LIBC_INLINE constexpr NumberPair<T> exact_add(T a, T b) {
return r;
}
+template <bool FAST2SUM = true, typename T = double>
+LIBC_INLINE constexpr NumberPair<T> exact_sub(T a, T b) {
+ NumberPair<T> r{0.0, 0.0};
+ if constexpr (FAST2SUM) {
+ r.hi = a - b;
+ T t = a - r.hi;
+ r.lo = t - b;
+ } else {
+ r.hi = a - b;
+ T t1 = r.hi - a;
+ T t2 = r.hi - t1;
+ T t3 = a - t2;
+ T t4 = b + t1;
+ r.lo = t3 - t4;
+ }
+ return r;
+}
+
// Following the above analysis of FastTwoSum,
// If `lsb(a.hi) >= ulp(b.hi)`, then the errors:
// err = (a.hi + a.lo + b.hi + b.lo) - (r.hi - r.lo) is bounded by:
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index c0751d7e15f102..f55bcf412bbffd 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -4293,10 +4293,11 @@ add_header_library(
DEPENDS
.exp_double_utils
.exp_utils
+ libc.src.__support.CPP.bit
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.fp_bitwise
libc.src.__support.FPUtil.multiply_add
- libc.src.__support.FPUtil.nearest_integer
libc.src.__support.common
libc.src.__support.macros.attributes
libc.src.__support.macros.config
@@ -4322,10 +4323,11 @@ add_header_library(
DEPENDS
.exp_double_utils
.exp_utils
+ libc.src.__support.CPP.bit
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.fp_bitwise
libc.src.__support.FPUtil.multiply_add
- libc.src.__support.FPUtil.nearest_integer
libc.src.__support.common
libc.src.__support.macros.attributes
libc.src.__support.macros.config
@@ -5734,6 +5736,75 @@ add_header_library(
libc.src.__support.macros.properties.types
)
+add_header_library(
+ sinhcosh_accurate_128
+ HDRS
+ sinhcosh_accurate_128.h
+ DEPENDS
+ .exp2_accurate_128
+ 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(
+ sinhcosh_utils
+ HDRS
+ sinhcosh_utils.h
+ DEPENDS
+ .exp_utils
+ .sinhcosh_accurate_128
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.fp_bitwise
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.nearest_integer
+ 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(
+ sinh_double_eval
+ HDRS
+ sinh_double_eval.h
+ DEPENDS
+ .sinhcosh_utils
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.fp_bitwise
+ 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.cpu_features
+)
+
+add_header_library(
+ sinh
+ HDRS
+ sinh.h
+ DEPENDS
+ .sinh_double_eval
+ libc.src.__support.common
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+)
+
add_header_library(
sinhf
HDRS
diff --git a/libc/src/__support/math/exp2_double_eval.h b/libc/src/__support/math/exp2_double_eval.h
index 8e2e372b58fbac..f45e20e30fe15d 100644
--- a/libc/src/__support/math/exp2_double_eval.h
+++ b/libc/src/__support/math/exp2_double_eval.h
@@ -14,10 +14,11 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_DOUBLE_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP2_DOUBLE_EVAL_H
+#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/fp_bitwise.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"
@@ -63,17 +64,23 @@ LIBC_ALWAYS_INLINE double exp2(double x) {
// Range reduction:
// Let x = k * 2^-N + dx
// in which:
- // k = round(x * 2^N)
- // hi = k >> N
+ // k = round(x * 2^N)
+ // hi = k >> N
// idx = k & (2^N - 1)
// then:
// 2^x = 2^hi * 2^(idx / 2^N) * 2^dx.
+ //
+ // Using ROUND_BIAS, ROUND_SHIFT, and ROUND_MASK (see exp_utils.h), we compute
+ // k_biased = fma(x, EXP_SCALE, ROUND_BIAS) and extract in parallel:
+ // kd = (k_biased & ROUND_MASK) - ROUND_SHIFT
+ // hi = (int32_t)(k_bits >> 19) >> (13 + N)
+ // idx = (k_bits >> 32) & (2^N - 1).
+ double k_biased = fputil::multiply_add(x, EXP_SCALE, ROUND_BIAS);
+ double kd = fputil::fp_and(k_biased, ROUND_MASK) - ROUND_SHIFT;
+ uint64_t k_bits = cpp::bit_cast<uint64_t>(k_biased);
- 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);
+ int hi = static_cast<int32_t>(k_bits >> 19) >> (13 + EXP_N_BITS);
+ unsigned idx = static_cast<unsigned>(k_bits >> 32) & EXP_MASK;
double dx = fputil::multiply_add(kd, -EXP_SCALE_INV, x); // exact
diff --git a/libc/src/__support/math/exp_double_eval.h b/libc/src/__support/math/exp_double_eval.h
index ef75ee37c5be62..b47736d0321f90 100644
--- a/libc/src/__support/math/exp_double_eval.h
+++ b/libc/src/__support/math/exp_double_eval.h
@@ -14,10 +14,11 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_EVAL_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_EXP_DOUBLE_EVAL_H
+#include "src/__support/CPP/bit.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/double_double.h"
+#include "src/__support/FPUtil/fp_bitwise.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"
@@ -70,17 +71,25 @@ LIBC_ALWAYS_INLINE double exp(double x) {
// Range reduction:
// Let x = k * (log(2) / 2^N) + u
// in which:
- // k = round(x * 2^N * log2(e))
- // hi = k >> N
+ // k = round(x * 2^N * log2(e))
+ // hi = k >> N
// idx = k & (2^N - 1)
// then:
// exp(x) = 2^hi * 2^(idx / 2^N) * exp(u).
- double x_scaled = x * (EXP_SCALE * LOG2_E.hi);
- double kd = fputil::nearest_integer(x_scaled);
- int k = static_cast<int>(kd);
+ //
+ // Using ROUND_BIAS, ROUND_SHIFT, and ROUND_MASK (see exp_utils.h), we compute
+ // k_biased = fma(x, LOG2_E_SCALE, ROUND_BIAS) and extract in parallel:
+ // kd = (k_biased & ROUND_MASK) - ROUND_SHIFT
+ // hi = (int32_t)(k_bits >> 19) >> (13 + N)
+ // idx = (k_bits >> 32) & (2^N - 1).
+ constexpr double LOG2_E_SCALE = EXP_SCALE * LOG2_E.hi;
+
+ double k_biased = fputil::multiply_add(x, LOG2_E_SCALE, ROUND_BIAS);
+ double kd = fputil::fp_and(k_biased, ROUND_MASK) - ROUND_SHIFT;
+ uint64_t k_bits = cpp::bit_cast<uint64_t>(k_biased);
- int hi = k >> EXP_N_BITS;
- unsigned idx = static_cast<unsigned>(k & EXP_MASK);
+ int hi = static_cast<int32_t>(k_bits >> 19) >> (13 + EXP_N_BITS);
+ unsigned idx = static_cast<unsigned>(k_bits >> 32) & EXP_MASK;
#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
double t1 = fputil::multiply_add(kd, MLOG_2_EXP2_HI, x); // exact
diff --git a/libc/src/__support/math/exp_utils.h b/libc/src/__support/math/exp_utils.h
index fa909c824efc5f..a3af63eabdaefa 100644
--- a/libc/src/__support/math/exp_utils.h
+++ b/libc/src/__support/math/exp_utils.h
@@ -161,6 +161,71 @@ LIBC_NOINLINE LIBC_INLINE double set_exceptional(double x) {
return x + FPBits::inf().get_val();
}
+// Fast range-reduction rounding constants for exp, exp2, sinh, and cosh:
+//
+// Instead of calling fputil::nearest_integer(x * SCALE) and
+// static_cast<int>(kd) (which has a serial dependency chain:
+// multiply -> roundsd/frintn -> cvttsd2si/fcvtzs),
+// we fuse the scaling and rounding bias into a single FMA:
+// k_biased = fma(x, SCALE, ROUND_BIAS).
+//
+// Derivation of ROUND_SHIFT, ROUND_BIAS, and ROUND_MASK:
+// 1) Why we cannot use 1.5 * 2^52:
+// Adding 1.5 * 2^52 (ulp = 1) relies on the hardware FMA rounding mode to
+// round z = x * SCALE to an integer. In directed rounding modes
+// (FE_DOWNWARD, FE_UPWARD, FE_TOWARDZERO), that rounds toward -inf/+inf/0
+// instead of nearest, doubling the reduced argument |u| from <= 0.5 * 2^-N
+// to <= 2^-N.
+// To round to nearest in all rounding modes, we keep F = 52 - M > 0
+// fractional bits in the 52-bit mantissa using:
+// ROUND_SHIFT = 1.5 * 2^M,
+// add +0.5 explicitly in ROUND_BIAS, and mask out the F fractional bits.
+//
+// 2) Choosing M = 20 (20 integer bits + 32 fractional bits):
+// For all non-overflowing/non-underflowing inputs of exp, exp2, sinh, and
+// cosh, the scaled input z = x * SCALE satisfies |z| < 2^18. To keep
+// z + 1.5 * 2^M inside the single binade:
+// [1.25 * 2^M, 1.75 * 2^M) subset [2^M, 2^(M + 1))
+// without carrying across 2^M, we need 2^18 <= 2^(M - 2), i.e., M >= 20.
+// Choosing the minimum M = 20 maximizes the fractional precision at
+// F = 52 - 20 = 32 fractional bits and aligns the binary point of k right
+// at the 32-bit word boundary (bit 32):
+// ROUND_SHIFT = 1.5 * 2^20 = 0x1.8p+20, ulp = 2^(20 - 52) = 2^-32.
+// In k_biased = z + ROUND_BIAS:
+// - The low 32 bits [31:0] hold the 32 fractional bits in [0, 1).
+// - Bits [50:32] hold the 19-bit signed two's-complement integer k (with
+// |k| < 2^18 and sign bit at bit 50).
+//
+// 3) Deriving the bias (+0.5 + 2^-32):
+// Masking out the low 32 fractional bits with:
+// ROUND_MASK = 0xffff'ffff'0000'0000
+// truncates toward -inf (floor). Adding +0.5 turns floor into nearest:
+// floor(z + 0.5) = round(z).
+// However, the FMA rounds z + ROUND_BIAS at ulp = 2^-32 before the mask.
+// With only +0.5, a tiny negative input (-2^-32 < z < 0) in FE_DOWNWARD or
+// FE_TOWARDZERO would round z + (1.5 * 2^20 + 0.5) down by 1 ULP to:
+// 1.5 * 2^20 + 0.5 - 2^-32,
+// which truncates to k = -1 instead of k = 0.
+// Adding +2^-32 (1 ULP at 2^20, i.e., mantissa bit 0) gives:
+// ROUND_BIAS = 1.5 * 2^20 + 0.5 + 2^-32 = 0x1.8000080000001p+20,
+// ensuring k_biased >= 1.5 * 2^20 (so k = 0 and kd = 0.0) for all
+// |x| < 2^-9 across all rounding modes.
+//
+// 4) Note on rounding errors:
+// Rounding the FMA...
[truncated]
``````````
</details>
https://github.com/llvm/llvm-project/pull/235097
More information about the libc-commits
mailing list