[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
Sun Oct 25 19:16:56 PDT 2026
https://github.com/lntue updated https://github.com/llvm/llvm-project/pull/235097
>From ba9c47cddca55882e6e64a8e03a359ce23dd8f71 Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 23 Oct 2026 20:48:33 +0000
Subject: [PATCH 1/3] [libc][math] Implement double precision sinh correctly
rounded for all rounding modes.
- 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.
---
libc/config/baremetal/aarch64/entrypoints.txt | 1 +
libc/config/baremetal/arm/entrypoints.txt | 1 +
libc/config/baremetal/riscv/entrypoints.txt | 1 +
libc/config/darwin/aarch64/entrypoints.txt | 1 +
libc/config/freebsd/x86_64/entrypoints.txt | 1 +
libc/config/gpu/amdgpu/entrypoints.txt | 1 +
libc/config/gpu/nvptx/entrypoints.txt | 1 +
libc/config/linux/aarch64/entrypoints.txt | 1 +
libc/config/linux/arm/entrypoints.txt | 1 +
libc/config/linux/riscv/entrypoints.txt | 1 +
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/config/windows/entrypoints.txt | 1 +
libc/docs/headers/math/index.md | 2 +-
libc/include/math.yaml | 6 +
libc/shared/math.h | 1 +
libc/shared/math/sinh.h | 28 ++
libc/src/__support/FPUtil/double_double.h | 18 +
libc/src/__support/math/CMakeLists.txt | 75 ++-
libc/src/__support/math/exp2_double_eval.h | 23 +-
libc/src/__support/math/exp_double_eval.h | 25 +-
libc/src/__support/math/exp_utils.h | 65 +++
libc/src/__support/math/sinh.h | 34 ++
libc/src/__support/math/sinh_double_eval.h | 237 +++++++++
.../__support/math/sinhcosh_accurate_128.h | 309 ++++++++++++
libc/src/__support/math/sinhcosh_utils.h | 467 ++++++++++++++++++
libc/src/math/generic/CMakeLists.txt | 10 +
libc/src/math/generic/sinh.cpp | 21 +
libc/test/shared/CMakeLists.txt | 1 +
libc/test/shared/shared_math_test.cpp | 1 +
libc/test/src/math/CMakeLists.txt | 16 +
libc/test/src/math/sinh_test.cpp | 187 +++++++
libc/test/src/math/smoke/CMakeLists.txt | 15 +
libc/test/src/math/smoke/sinh_test.cpp | 123 +++++
.../llvm-project-overlay/libc/BUILD.bazel | 80 ++-
.../libc/test/src/math/BUILD.bazel | 8 +
.../libc/test/src/math/smoke/BUILD.bazel | 8 +
36 files changed, 1751 insertions(+), 21 deletions(-)
create mode 100644 libc/shared/math/sinh.h
create mode 100644 libc/src/__support/math/sinh.h
create mode 100644 libc/src/__support/math/sinh_double_eval.h
create mode 100644 libc/src/__support/math/sinhcosh_accurate_128.h
create mode 100644 libc/src/__support/math/sinhcosh_utils.h
create mode 100644 libc/src/math/generic/sinh.cpp
create mode 100644 libc/test/src/math/sinh_test.cpp
create mode 100644 libc/test/src/math/smoke/sinh_test.cpp
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 at ulp(k_biased) = 2^-32 before truncating the 32
+// fractional bits gives:
+// k = floor(z + 0.5 + 2^-32 + delta), |delta| < 2^-32.
+// Since k (extracted from bits [50:32]) and:
+// kd = (k_biased & ROUND_MASK) - ROUND_SHIFT
+// are still equal integers, the range-reduction identity remains exact; the
+// extra rounding error only shifts the half-integer rounding boundary by at
+// most 2^-31, widening the reduced argument |u| by at most a factor of
+// (1 + 2^-30) and increasing the polynomial approximation error by < 2^-90.
+LIBC_INLINE_VAR constexpr double ROUND_BIAS = 0x1.8000080000001p+20;
+LIBC_INLINE_VAR constexpr double ROUND_SHIFT = 0x1.8p+20;
+LIBC_BIT_CAST_CONSTEXPR_VAR double ROUND_MASK =
+ cpp::bit_cast<double>(0xffff'ffff'0000'0000ULL);
+
#if !defined(LIBC_MATH_HAS_SKIP_ACCURATE_PASS) || \
!defined(LIBC_MATH_HAS_SMALL_TABLES)
diff --git a/libc/src/__support/math/sinh.h b/libc/src/__support/math/sinh.h
new file mode 100644
index 00000000000000..a3252dfe862c19
--- /dev/null
+++ b/libc/src/__support/math/sinh.h
@@ -0,0 +1,34 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 sinh.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINH_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_SINH_H
+
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+
+#include "src/__support/math/sinh_double_eval.h"
+#define LIBC_MATH_SINH_IMPL double_eval
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+using LIBC_MATH_SINH_IMPL::sinh;
+
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#undef LIBC_MATH_SINH_IMPL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINH_H
diff --git a/libc/src/__support/math/sinh_double_eval.h b/libc/src/__support/math/sinh_double_eval.h
new file mode 100644
index 00000000000000..d0de9b4aaf60aa
--- /dev/null
+++ b/libc/src/__support/math/sinh_double_eval.h
@@ -0,0 +1,237 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 sinh.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINH_DOUBLE_EVAL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_SINH_DOUBLE_EVAL_H
+
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/CPP/bit.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/fp_bitwise.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/cpu_features.h"
+#include "src/__support/math/sinhcosh_utils.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace double_eval {
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// Leading 11-bit split of c_0 ~ 1/6 for exact cubic reconstruction:
+// SINH_C0_HI_11 = round(1/6, 11, RD) = 0x1.554p-3
+// 1/6 - SINH_C0_HI_11 = 0x1.5555555555555p-15.
+LIBC_INLINE_VAR constexpr double SINH_C0_HI_11 = 0x1.554p-3;
+
+// Polynomial approximation for ((sinh(x) - x) / x^3 - 0x1.554p-3) on
+// [-2^-2, 2^-2] as a degree-5 even polynomial in x^2:
+// sinh(x) - x ~ x^3 * (SINH_C0_HI_11 + P_rem(x^2))
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > display = hexadecimal!;
+// > prec = 500;
+// > P = fpminimax((sinh(sqrt(x))/sqrt(x) - 1)/x - 0x1.554p-3, 5, [|D...|],
+// [2^(-100), 2^(-4)]);
+// > for i from 0 to 5 do print(coeff(P, i), ",");
+// > dirtyinfnorm(sinh(x) - (x + x^3 * (0x1.554p-3 + subst(P, x, x^2))),
+// [-2^(-2), 2^(-2)]);
+// 0x1.0315...p-75 < 2^-74.
+LIBC_INLINE_VAR constexpr double SINH_SMALL_P_COEFFS[6] = {
+ 0x1.5555555555555p-15, // c_0 - 0x1.554p-3 ~ 1/6 - 0x1.554p-3
+ 0x1.1111111111111p-7, // c_1 ~ 1/120
+ 0x1.a01a01a01a38ap-13, // c_2 ~ 1/5040
+ 0x1.71de3a52826f3p-19, // c_3 ~ 1/362880
+ 0x1.ae6471bbce233p-26, // c_4 ~ 1/39916800
+ 0x1.60eca3b25ae29p-33, // c_5 ~ 1/6227020800
+};
+
+// Evaluate P_rem(x^2) ~ (sinh(x) - x) / x^3 - 0x1.554p-3 using Estrin's
+// scheme:
+LIBC_ALWAYS_INLINE double sinh_small_poly(double x2) {
+ double x4 = x2 * x2;
+ double x8 = x4 * x4;
+ double c0 =
+ fputil::multiply_add(x2, SINH_SMALL_P_COEFFS[1], SINH_SMALL_P_COEFFS[0]);
+ double c1 =
+ fputil::multiply_add(x2, SINH_SMALL_P_COEFFS[3], SINH_SMALL_P_COEFFS[2]);
+ double c2 =
+ fputil::multiply_add(x2, SINH_SMALL_P_COEFFS[5], SINH_SMALL_P_COEFFS[4]);
+ double c01 = fputil::multiply_add(x4, c1, c0);
+ return fputil::multiply_add(x8, c2, c01);
+}
+
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// Evaluate sinh(x) for |x| < SMALL_BOUND (2^-26 + ulp(2^-26) when
+// LIBC_MATH_HAS_SKIP_ACCURATE_PASS is defined, or 2^-2 otherwise):
+LIBC_INLINE double sinh_small_eval(double x, uint64_t x_abs) {
+ // Tiny and subnormal range: |x| <= 2^-26.
+ // For 0 < |x| <= 2^-26:
+ // 0 < (sinh(|x|) - |x|) / |x| ~ x^2 / 6 < 2^-54 < 2^-1 * ulp(x),
+ // so sinh(x) rounds to:
+ // - nextafter(x, +inf) when x > 0 and rounding mode is FE_UPWARD,
+ // - nextafter(x, -inf) when x < 0 and rounding mode is FE_DOWNWARD,
+ // - x otherwise, with FE_INEXACT raised for all x != 0 and FE_UNDERFLOW +
+ // errno = ERANGE raised when |x| < 2^-1022.
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ if (LIBC_UNLIKELY(x_abs <= 0x3e50'0000'0000'0000ULL))
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ {
+ if (LIBC_UNLIKELY(x_abs == 0))
+ return x;
+ if (LIBC_UNLIKELY(x_abs < 0x0010'0000'0000'0000ULL)) {
+ fputil::set_errno_if_required(ERANGE);
+ return fputil::multiply_add(x, 0x1.0p-55, x);
+ }
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ return fputil::multiply_add(x, 0x1.0p-55, x);
+#else
+ // Avoid spurious FE_UNDERFLOW on x * 2^-55 when 2^-1022 <= |x| < 2^-967:
+ double x_scaled = fputil::multiply_add(x, 0x1.0p55, x);
+ return x_scaled * 0x1.0p-55;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ }
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // 2^-26 < |x| < 2^-2:
+ // Split x = x_hi + x_lo into 14-bit x_hi and 39-bit x_lo so that:
+ // x2_hi = x_hi * x_hi (28 bits)
+ // x3_hi = x2_hi * x_hi (42 bits)
+ // y_hi = x3_hi * SINH_C0_HI_11 (42 + 11 = 53 bits)
+ // are all exact in double precision across all rounding modes (with or
+ // without FMA).
+ LIBC_BIT_CAST_CONSTEXPR_VAR double HI_MASK_14 =
+ cpp::bit_cast<double>(0xffff'ff80'0000'0000ULL);
+ double x_hi = fputil::fp_and(x, HI_MASK_14);
+ double x_lo = x - x_hi;
+ double x2 = x * x;
+ double x2_hi = x_hi * x_hi; // exact (28 bits)
+ double x3_hi = x2_hi * x_hi; // exact (42 bits)
+ double y_hi = x3_hi * SINH_C0_HI_11; // exact (53 bits)
+
+ // Cubic remainder:
+ // x3_lo ~ x^3 - x_hi^3 = x_lo * (x^2 + x * x_hi + x_hi^2)
+ double x2_mid = fputil::multiply_add(x, x_hi, x2);
+ double x3_sum = x2_mid + x2_hi;
+ double x3_lo = x_lo * x3_sum;
+
+ // Polynomial remainder:
+ // x^3 * (SINH_C0_HI_11 + P_rem(x^2)) = y_hi + y_lo
+ // where:
+ // y_lo ~ x_hi^3 * P_rem(x^2) + x3_lo * (SINH_C0_HI_11 + P_rem(x^2))
+ double p_rem = sinh_small_poly(x2);
+ double p_full = SINH_C0_HI_11 + p_rem;
+ double y_lo_t = x3_lo * p_full;
+ double y_lo = fputil::multiply_add(x3_hi, p_rem, y_lo_t);
+
+ // Leading sum r_hi + s_err = x + y_hi (exact FastTwoSum in all rounding modes
+ // since |y_hi| < 2^-6 * |x|):
+ double r_hi = x + y_hi;
+ double r_t = x - r_hi;
+ double s_err = r_t + y_hi;
+
+ // Error bound for Ziv's rounding test:
+ // The error of y_hi + y_lo is bounded by |y_hi| * 2^-57, and the rounding
+ // error of s_err + y_lo is bounded by |x| * 2^-103.
+ double x_err = x * 0x1.0p-103;
+ double err = fputil::multiply_add(y_hi, 0x1.0p-57, x_err);
+
+ // Low tail bounds lo_{u,l} ~ s_err + (y_lo +/- err):
+ double y_lo_u = y_lo + err;
+ double y_lo_l = y_lo - err;
+ double lo_u = s_err + y_lo_u;
+ double lo_l = s_err + y_lo_l;
+
+ // 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))
+ return upper;
+
+ return sinh_small_accurate_128(x);
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+}
+
+// Handle overflow, inf, and nan for sinh(x) when |x| >= 0x1.633ce8fb9f87ep+9:
+LIBC_NOINLINE LIBC_INLINE double sinh_exceptional(double x) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ if (xbits.is_nan())
+ return x + 1.0;
+
+ if (xbits.is_inf())
+ return x;
+
+ 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 (xbits.is_neg()) {
+ if (LIBC_UNLIKELY(rounding == FE_UPWARD || rounding == FE_TOWARDZERO))
+ return FPBits::max_normal(Sign::NEG).get_val();
+ } else {
+ if (LIBC_UNLIKELY(rounding == FE_DOWNWARD || rounding == FE_TOWARDZERO))
+ return FPBits::max_normal(Sign::POS).get_val();
+ }
+#endif // !LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+
+ return FPBits::inf(xbits.sign()).get_val();
+}
+
+LIBC_ALWAYS_INLINE double sinh(double x) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ uint64_t x_abs = xbits.abs().uintval();
+
+ // Bound constants (with biased exponent and mantissa bits):
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // - Small bound: 2^-26 + ulp(2^-26) (exp_pm_eval handles |x| > 2^-26 directly
+ // with k = 0 and zero cancellation when |x| < log(2) / 2)
+ constexpr uint64_t SMALL_BOUND = 0x3e50'0000'0000'0001ULL;
+#else
+ // - Small bound: 2^-2 = 0x1.0p-2
+ constexpr uint64_t SMALL_BOUND = 0x3fd0'0000'0000'0000ULL;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // - Overflow bound:
+ // > round(asinh(2^1023 * (2 - 2^-52)), D, RU) = 0x1.633ce8fb9f87ep+9
+ // > round(asinh(2^1023 * (2 - 2^-52)), D, RD) = 0x1.633ce8fb9f87dp+9
+ // > round(sinh(0x1.633ce8fb9f87dp+9), D, RN) = 0x1.fffffffffffffp+1023
+ // > round(sinh(0x1.633ce8fb9f87ep+9), D, RN) = +inf
+ constexpr uint64_t OVERFLOW_BOUND = 0x4086'33ce'8fb9'f87eULL;
+
+ if (LIBC_UNLIKELY(x_abs < SMALL_BOUND || x_abs >= OVERFLOW_BOUND)) {
+ if (x_abs < SMALL_BOUND)
+ return sinh_small_eval(x, x_abs);
+
+ return sinh_exceptional(x);
+ }
+
+ return exp_pm_eval(x, /*is_sinh=*/true);
+}
+
+} // namespace double_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINH_DOUBLE_EVAL_H
diff --git a/libc/src/__support/math/sinhcosh_accurate_128.h b/libc/src/__support/math/sinhcosh_accurate_128.h
new file mode 100644
index 00000000000000..58fdbf1250de6d
--- /dev/null
+++ b/libc/src/__support/math/sinhcosh_accurate_128.h
@@ -0,0 +1,309 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 paths for double-precision sinh and cosh functions.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_ACCURATE_128_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_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/frac64.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/math/exp2_accurate_128.h"
+#include "src/__support/uint128.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace double_eval {
+
+// Polynomial approximation for (sinh(x) - x) / x^3 on [-2^-2, 2^-2] as a
+// degree-9 polynomial in x^2 in 128-bit fixed-point format:
+// P(x^2) ~ (sinh(x) - x) / x^3
+// sinh(x) = x + x^3 * P(x^2).
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > display = hexadecimal!;
+// > prec = 500;
+// > P = fpminimax((sinh(sqrt(x))/sqrt(x) - 1)/x, 9, [|128...|],
+// [2^(-150), 2^(-4)], fixed);
+// > for i from 0 to 9 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(sinh(x) - (x + x^3 * subst(P, x, x^2)), [-2^-2, 2^-2]);
+// 0x1.1cb4...p-142 < 2^-141.
+LIBC_INLINE_VAR constexpr Frac128 SINH_SMALL_POLY_128[10] = {
+ // lsb = 2^-128
+ Frac128({0xaaaaaaaaaaaaaaabULL, 0x2aaaaaaaaaaaaaaaULL}), // c_0
+ Frac128({0x2222222222222115ULL, 0x0222222222222222ULL}), // c_1
+ Frac128({0x00d00d00d00de2b5ULL, 0x000d00d00d00d00dULL}), // c_2
+ Frac128({0x671f558390d8d543ULL, 0x00002e3bc74aad8eULL}), // c_3
+ Frac128({0x138e3f9d219de9cbULL, 0x0000006b99159fd5ULL}), // c_4
+ Frac128({0x43684be86ba485efULL, 0x00000000b092309dULL}), // c_5
+ Frac128({0x399dc020c9d4a0d9ULL, 0x0000000000d73f9fULL}), // c_6
+ Frac128({0x3b819e88de5f1579ULL, 0x000000000000ca96ULL}), // c_7
+ Frac128({0xa4d89e719113ac86ULL, 0x0000000000000097ULL}), // c_8
+ Frac128({0x5c7ba5e9742ae90bULL, 0x0000000000000000ULL}), // c_9
+};
+
+// 128-bit accurate path for sinh(x) when 2^-26 < |x| < 2^-2:
+LIBC_NOINLINE LIBC_INLINE double sinh_small_accurate_128(double x) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ int exp_x = xbits.get_exponent(); // in [-26, -3]
+
+ // Normalize |x|'s 53-bit mantissa to m_x in [1, 2) with lsb = 2^-63:
+ // |x| = 2^exp_x * m_x.
+ uint64_t mant_x = (xbits.get_mantissa() | (1ULL << 52)) << 11; // lsb = 2^-63
+
+ // Normalized mantissa of |x| and m_x^2:
+ // m_x_frac = m_x / 2 (lsb = 2^-127)
+ // m_x2 = mant_x^2 * 2^-126 (lsb = 2^-126)
+ Frac128 m_x_frac({0, mant_x}); // lsb = 2^-127
+ UInt<128> m_x2 = UInt<64>(mant_x).ful_mul(UInt<64>(mant_x)); // lsb = 2^-126
+
+ // Since Frac128(m_x2.val) represents m_x^2 / 4 (lsb = 2^-128) and m_x_frac
+ // represents m_x / 2 (lsb = 2^-127), their Frac128 product (>> 128)
+ // represents m_x^3 / 8 with lsb = 2^-125:
+ Frac128 m_x3 = Frac128(m_x2.val) * m_x_frac; // lsb = 2^-125
+
+ // Scale m_x^2 by 2^(2 * exp_x) to obtain x^2 = 2^(2 * exp_x) * m_x^2 with
+ // lsb = 2^-128. Since exp_x is in [-26, -3], shift = -2 * exp_x - 2 is in
+ // [4, 50]:
+ unsigned shift = static_cast<unsigned>(-2 * exp_x - 2); // in [4, 50]
+ Frac128 x2_frac = Frac128(m_x2.val) >> shift; // lsb = 2^-128
+
+ // Evaluate P(x^2) ~ (sinh(x) - x) / x^3 in Horner form:
+ Frac128 p = SINH_SMALL_POLY_128[9]; // lsb = 2^-128
+ LIBC_LOOP_NOUNROLL
+ for (int i = 8; i >= 0; --i)
+ p = p * x2_frac + SINH_SMALL_POLY_128[i]; // lsb = 2^-128
+
+ // Reconstruction of |sinh(x)| = |x| + |x|^3 * P(x^2):
+ // Since y_poly = m_x3 * p represents m_x^3 * P(x^2) with lsb = 2^-125,
+ // shifting right by shift = -2 * exp_x - 2 aligns:
+ // 2^(3 * exp_x) * m_x^3 * P(x^2)
+ // to the same scale as m_x_frac (lsb = 2^(exp_x - 127)):
+ Frac128 y_poly = m_x3 * p; // lsb = 2^-125
+ Frac128 y_m = y_poly >> shift; // lsb = 2^(exp_x - 127)
+
+ // Add |x| and |x|^3 * P(x^2), normalizing on carry out:
+ UInt<128> x_m = m_x_frac; // lsb = 2^(exp_x - 127)
+ int exp_res = exp_x - 127;
+ if (LIBC_UNLIKELY(x_m.add_overflow(y_m) != 0)) {
+ x_m = (x_m >> 1) | (UInt<128>(1) << 127);
+ ++exp_res;
+ }
+
+ DFloat128 r;
+ r.sign = xbits.sign();
+ r.exponent = exp_res;
+ r.mantissa = x_m;
+ return static_cast<double>(r);
+}
+
+// Polynomial approximations for sinh(u * log(2)) and cosh(u * log(2)) - 1 on
+// u in [0, 2^-5] in terms of v = 32 * u in [0, 1] and w = v^2 in [0, 1] in
+// 128-bit fixed-point format (lsb = 2^-128):
+// sinh(u * log(2)) = sinh(2^-5 * v * log(2)) ~ v * P_odd(w)
+// cosh(u * log(2)) - 1 = cosh(2^-5 * v * log(2)) - 1 ~ w * P_even(w).
+//
+// For k >= 4, the coefficients c_k of P_odd(w) and P_even(w) satisfy
+// c_k < 2^-64, so their high 64-bit words val[1] are identically zero and are
+// stored as 64-bit low words in EXP2_SINH_W_LO64 and EXP2_COSH_W_LO64.
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > prec = 500;
+// > P_odd = fpminimax(sinh(2^-5 * sqrt(x) * log(2)) / sqrt(x), 6, [|128...|],
+// [2^(-150), 1], fixed);
+// > for i from 0 to 6 do {
+// c = coeff(P_odd, 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}), // odd_", i);
+// };
+// > dirtyinfnorm(sinh(2^-5 * sqrt(x) * log(2)) - sqrt(x) * P_odd, [0, 1]);
+// 0x1.6276...p-132 < 2^-131.
+LIBC_INLINE_VAR constexpr Frac128 EXP2_SINH_W_128[4] = {
+ // lsb = 2^-128
+ Frac128({0x5e4f1d9cc01f97b5ULL, 0x058b90bfbe8e7bcdULL}), // odd_0
+ Frac128({0xbf8b33a762bb32caULL, 0x00001c6b08d704a0ULL}), // odd_1
+ Frac128({0x14ce6220e2fed2eaULL, 0x000000002bb0ffcfULL}), // odd_2
+ Frac128({0xbfc588b0c686b280ULL, 0x0000000000001ffcULL}), // odd_3
+};
+
+LIBC_INLINE_VAR constexpr Frac64 EXP2_SINH_W_LO64[3] = {
+ // lsb = 2^-128 (val[1] == 0)
+ Frac64(0x0da929e9caf3e033ULL), // odd_4
+ Frac64(0x000003d1958e6b49ULL), // odd_5
+ Frac64(0x0000000000c0b0b4ULL), // odd_6
+};
+
+// > P_even = fpminimax((cosh(2^-5 * sqrt(x) * log(2)) - 1) / x, 6, [|128...|],
+// [2^(-150), 1], fixed);
+// > for i from 0 to 6 do {
+// c = coeff(P_even, 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}), // even_", i);
+// };
+// > dirtyinfnorm((cosh(2^-5 * sqrt(x) * log(2)) - 1) - x * P_even, [0, 1]);
+// 0x1.a7d4...p-130 < 2^-129.
+LIBC_INLINE_VAR constexpr Frac128 EXP2_COSH_W_128[4] = {
+ // lsb = 2^-128
+ Frac128({0x54378b583764b9b0ULL, 0x000f5fdeffc162c7ULL}), // even_0
+ Frac128({0x9cee539977c16a7dULL, 0x000000276556df74ULL}), // even_1
+ Frac128({0x5f0d8f0ede96150bULL, 0x0000000000286122ULL}), // even_2
+ Frac128({0x2c0223a5c823fdddULL, 0x0000000000000016ULL}), // even_3
+};
+
+LIBC_INLINE_VAR constexpr Frac64 EXP2_COSH_W_LO64[3] = {
+ // lsb = 2^-128 (val[1] == 0)
+ Frac64(0x0007933d4562e2e0ULL), // even_4
+ Frac64(0x00000001c3bd65faULL), // even_5
+ Frac64(0x0000000000004bf6ULL), // even_6
+};
+
+// Compute:
+// sinh(x) = (e^x - e^-x) / 2 (when is_sinh is true)
+// cosh(x) = (e^x + e^-x) / 2 (when is_sinh is false)
+// in 128-bit precision using direct base-2 192-bit fixed-point range reduction
+// for 2^-2 <= |x| < 0x1.633ce8fb9f87ep+9 (sinh) or
+// 2^-27 < |x| < 0x1.633ce8fb9f87ep+9 (cosh).
+LIBC_NOINLINE LIBC_INLINE double exp_pm_accurate_128(double x, bool is_sinh) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ int exp_x = xbits.get_exponent(); // in [-2, 9] for sinh, [-27, 9] for cosh
+
+ // Normalize |x|'s 53-bit mantissa to [2^63, 2^64):
+ uint64_t mant_x = (xbits.get_mantissa() | (1ULL << 52)) << 11;
+
+ // 1-limb x 3-limb product for |z| = |x| * log2(e):
+ // mant_x has lsb = 2^(exp_x - 63)
+ // LOG2_E_192.log2_b has lsb = 2^-191
+ // p = mant_x * LOG2_E_192.log2_b has lsb = 2^(exp_x - 254)
+ UInt<256> p = LOG2_E_192.log2_b.ful_mul(UInt<64>(mant_x));
+
+ // Nearest-integer range reduction:
+ // k_p = round(|z| * 2^4), u = |z| - k_p * 2^-4 in [-2^-5, 2^-5]
+ // v = 32 * |u| in [0, 1) (lsb = 2^-128)
+ // w = v^2 in [0, 1) (lsb = 2^-128).
+ int k_p = 0;
+ bool u_is_neg = false;
+ Frac128 v;
+ if (LIBC_LIKELY(exp_x >= -6)) {
+ // Since p.val[3] has lsb = 2^(exp_x - 62), the 2^-5 bit of |z| is at bit
+ // shr = 57 - exp_x (in [48, 63]) of p.val[3]:
+ int shr = 57 - exp_x;
+ uint64_t k_twice = p.val[3] >> shr; // floor(|z| * 2^5), lsb = 2^-5
+ k_p = static_cast<int>((k_twice + 1) >> 1);
+
+ // Shift p left by shl = 64 - shr = exp_x + 7 (in [1, 16]) so that the top
+ // two limbs hold the fractional bits of 32 * (|z| - k_p * 2^-4) with
+ // lsb = 2^-128:
+ int shl = exp_x + 7;
+ uint64_t v_lo = (p.val[2] << shl) | (p.val[1] >> shr);
+ uint64_t v_hi = (p.val[3] << shl) | (p.val[2] >> shr);
+ Frac128 u_bits({v_lo, v_hi});
+
+ // When u < 0 (k_twice is odd), one's-complement negation (~u_bits) gives
+ // v = 32 * |u| in [0, 1) to within 2^-128:
+ u_is_neg = (k_twice & 1ULL) != 0;
+ v = u_is_neg ? ~u_bits : u_bits; // lsb = 2^-128
+ } else {
+ // For cosh with exp_x in [-27, -7], |z| < 2^-5 so k_p = 0, u = |z| >= 0,
+ // and v = 32 * |z| is obtained by shifting the top 128 bits of p right by
+ // (-7 - exp_x):
+ size_t shift_v = static_cast<size_t>(-7 - exp_x);
+ Frac128 p_hi({p.val[2], p.val[3]});
+ v = p_hi >> shift_v; // lsb = 2^-128
+ }
+ Frac128 w = v * v; // lsb = 2^-128
+
+ // Decompose k_p = 16 * hi_p + idx_p and k_m = -k_p = 16 * hi_m + idx_m:
+ // (e^|x| -/+ e^-|x|) / 2 = 2^hi_p * (m_p * 2^u -/+ m_m * 2^-u)
+ // where:
+ // m_p = 2^(idx_p / 16) / 2 in [2^-1, 1) (lsb = 2^-127)
+ // m_m = 2^(hi_m - hi_p) * 2^(idx_m / 16) / 2 < 2^-1 (lsb = 2^-127).
+ int k_m = -k_p;
+ int hi_p = k_p >> 4;
+ unsigned idx_p = static_cast<unsigned>(k_p & 0x0f);
+ int hi_m = k_m >> 4;
+ unsigned idx_m = static_cast<unsigned>(k_m & 0x0f);
+
+ Frac128 m_p = EXP2_MID_FRAC128[idx_p] >> 1; // lsb = 2^-127
+ unsigned shift_m = static_cast<unsigned>(hi_p - hi_m + 1);
+ Frac128 m_m(0);
+ if (LIBC_LIKELY(shift_m < 128))
+ m_m = EXP2_MID_FRAC128[idx_m] >> shift_m; // lsb = 2^-127
+
+ // Split 2^(+/- u) = cosh(u * log(2)) +/- sinh(u * log(2)):
+ // m_p * 2^u -/+ m_m * 2^-u
+ // = s * (1 + w * P_even(w)) + sign(u) * d * v * P_odd(w)
+ // where s = m_p -/+ m_m and d = m_p +/- m_m are both in (2^-2, 2) with
+ // lsb = 2^-127.
+ Frac128 s = is_sinh ? (m_p - m_m) : (m_p + m_m); // lsb = 2^-127
+ Frac128 d = is_sinh ? (m_p + m_m) : (m_p - m_m); // lsb = 2^-127
+
+ Frac128 d_v = d * v; // lsb = 2^-127
+ Frac128 s_w = s * w; // lsb = 2^-127
+
+ // High-degree Horner steps (k = 6, 5, 4 -> 3) using 64-bit multiplications
+ // since the high 64-bit words val[1] of c_4, c_5, c_6 are zero:
+ Frac64 w_hi = w.to_frac64(); // lsb = 2^-64
+ Frac64 odd_lo = EXP2_SINH_W_LO64[2];
+ Frac64 even_lo = EXP2_COSH_W_LO64[2];
+ for (int i = 1; i >= 0; --i) {
+ odd_lo = odd_lo * w_hi + EXP2_SINH_W_LO64[i];
+ even_lo = even_lo * w_hi + EXP2_COSH_W_LO64[i];
+ }
+ odd_lo = odd_lo * w_hi + Frac64(EXP2_SINH_W_128[3].val[0]);
+ even_lo = even_lo * w_hi + Frac64(EXP2_COSH_W_128[3].val[0]);
+
+ // Low-degree Horner steps (k = 2, 1, 0) in 128-bit fixed-point:
+ Frac128 p_odd({odd_lo.val[0], EXP2_SINH_W_128[3].val[1]}); // lsb = 2^-128
+ Frac128 p_even({even_lo.val[0], EXP2_COSH_W_128[3].val[1]}); // lsb = 2^-128
+ LIBC_LOOP_NOUNROLL
+ for (int i = 2; i >= 0; --i) {
+ p_odd = p_odd * w + EXP2_SINH_W_128[i];
+ p_even = p_even * w + EXP2_COSH_W_128[i];
+ }
+
+ // Reconstruct:
+ // r_m = s + s_w * P_even(w) + sign(u) * d_v * P_odd(w)
+ // in (2^-2, 2) with lsb = 2^-127:
+ Frac128 even_part = s + s_w * p_even;
+ Frac128 odd_part = d_v * p_odd;
+ Frac128 r_m = u_is_neg ? (even_part - odd_part)
+ : (even_part + odd_part); // lsb = 2^-127
+
+ // Normalize r_m and round 2^hi_p * r_m (with lsb = 2^(hi_p - 127)) to double
+ // in the current rounding mode:
+ Sign sign = is_sinh ? xbits.sign() : Sign::POS;
+ DFloat128 r(sign, hi_p - 127, r_m);
+ return static_cast<double>(r);
+}
+
+} // namespace double_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_ACCURATE_128_H
diff --git a/libc/src/__support/math/sinhcosh_utils.h b/libc/src/__support/math/sinhcosh_utils.h
new file mode 100644
index 00000000000000..f93d5b56a9d181
--- /dev/null
+++ b/libc/src/__support/math/sinhcosh_utils.h
@@ -0,0 +1,467 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 for double-precision sinh and cosh functions.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_UTILS_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_UTILS_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"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/exp_utils.h"
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#include "src/__support/math/sinhcosh_accurate_128.h"
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+namespace double_eval {
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// Splitting -log(2) ~ hi + lo with 35-bit-precision hi part so that
+// kd * MLOG_2_HI is exact for |kd| < 2^18:
+// > hi = round(-log(2), 35, RN);
+// > lo = round(-log(2) - hi, D, RN);
+LIBC_INLINE_VAR constexpr double MLOG_2_HI = -0x1.62e42fefcp-1;
+LIBC_INLINE_VAR constexpr double MLOG_2_LO = 0x1.c610ca86c3899p-37;
+
+// Degree-4 minimax polynomials in u^2 for (cosh(u) - 1) / u^2 and
+// (sinh(u) - u) / u^3 on [-log(2) / 2, log(2) / 2]:
+// cosh(u) - 1 ~ u^2 * P_even(u^2)
+// sinh(u) - u ~ u^3 * P_odd(u^2).
+//
+// Generated by Sollya with:
+// > display = hexadecimal!;
+// > prec = 500;
+// > B = (log(2) / 2)^2;
+// > P_even = fpminimax((cosh(sqrt(x)) - 1) / x, 4, [|D...|], [2^(-100), B]);
+// > for i from 0 to 4 do print(coeff(P_even, i), ",");
+// > dirtyinfnorm((cosh(sqrt(x)) - 1) - x * P_even, [0, B]);
+// 0x1.c244...p-57 < 2^-56.
+LIBC_INLINE_VAR constexpr double EXP_EVEN_N0_COEFFS[5] = {
+ 0x1.0000000000001p-1, // ~ 1/2!
+ 0x1.5555555553ce4p-5, // ~ 1/4!
+ 0x1.6c16c178a71cep-10, // ~ 1/6!
+ 0x1.a019b8c6ec68dp-16, // ~ 1/8!
+ 0x1.2891c21cf0c03p-22, // ~ 1/10!
+};
+
+// > P_odd = fpminimax((sinh(sqrt(x)) / sqrt(x) - 1) / x, 4, [|D...|],
+// [2^(-100), B]);
+// > for i from 0 to 4 do print(coeff(P_odd, i), ",");
+// > dirtyinfnorm((sinh(sqrt(x)) - sqrt(x)) - x * sqrt(x) * P_odd, [0, B]);
+// 0x1.4882...p-62 < 2^-61.
+LIBC_INLINE_VAR constexpr double EXP_ODD_N0_COEFFS[5] = {
+ 0x1.5555555555556p-3, // ~ 1/3!
+ 0x1.111111111070cp-7, // ~ 1/5!
+ 0x1.a01a01a901652p-13, // ~ 1/7!
+ 0x1.71de0a3fd7069p-19, // ~ 1/9!
+ 0x1.af3effc3463e9p-26, // ~ 1/11!
+};
+
+// Tableless (N = 0) evaluation of:
+// sinh(x) = (e^x - e^-x) / 2 (when is_sinh is true)
+// cosh(x) = (e^x + e^-x) / 2 (when is_sinh is false, with x >= 0)
+// for 2^-26 < |x| < 0x1.633ce8fb9f87ep+9 (sinh) or
+// 2^-27 < x < 0x1.633ce8fb9f87ep+9 (cosh).
+//
+// Range reduction:
+// Let |x| = k * log(2) + u, where:
+// k = round(|x| * log2(e))
+// |u| <= log(2) / 2
+// then:
+// e^|x| = 2^k * e^u
+// e^-|x| = 2^-k * e^-u = 2^k * 2^(-2 * k) * e^-u.
+//
+// Factoring out 2^(k - 1) and separating e^(+/- u) = cosh(u) +/- sinh(u):
+// (e^|x| -/+ e^-|x|) / 2 = 2^(k - 1) * (e^u + m_m * e^-u)
+// = 2^(k - 1) * (s * cosh(u) + d * sinh(u))
+// where:
+// m_m = (-/+ 1) * 2^(-2 * k)
+// s = 1 + m_m
+// d = 1 - m_m.
+//
+// Key properties across the domain:
+// - For |x| < 0x1.633ce8fb9f87ep+9, k is in [0, 1025] so k - 1 is in
+// [-1, 1024]. When k = 1025, u < 0 keeps:
+// r = e^u + m_m * e^-u < 1.0,
+// so 2^1024 * r stays in the normal range and scale_exp2(r, k - 1) handles
+// the entire domain without a separate top-octave branch.
+//
+// - For |x| < log(2) / 2, k = 0 gives u = |x| (exact, t2 = 0) and m_m = -/+ 1:
+// * For sinh (m_m = -1): s = 0 and d = 2, reducing the evaluation to:
+// 2^-1 * 2 * (u + u^3 * P_odd(u^2)) = x + x^3 * P_odd(x^2)
+// with no cancellation down to |x| = 2^-26.
+// * For cosh (m_m = +1): s = 2 and d = 0, reducing to:
+// 1 + u^2 * P_even(u^2).
+//
+// - For 1 <= k <= 26, 2^(-2 * k) is in [2^-52, 2^-2], so both:
+// s = 1.0 + m_m
+// d = 1.0 - m_m
+// are exact in all rounding modes. For k >= 27:
+// |m_m| <= 2^-54 < 0.5 * ulp(1.0),
+// so clamping m_m to 0.0 avoids a 1-ULP rounding error in directed rounding
+// modes.
+LIBC_MATH_SMALL_TABLES_INLINE double exp_pm_eval(double x, bool is_sinh) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ double x_a = xbits.abs().get_val();
+
+ double x_scaled = x_a * LOG2_E.hi;
+ double kd = fputil::nearest_integer(x_scaled);
+ int k = static_cast<int>(kd);
+
+ // Secondary scale factor:
+ // m_m = (-/+ 1) * 2^(-2 * k) for k < 27
+ // m_m = 0.0 for k >= 27 (where 2^(-2 * k) <= 2^-54).
+ //
+ // Subtracting (1 << 11) = 0x800 (for sinh) sets bit 11 of the 12-bit
+ // two's-complement exponent field, flipping the sign bit of m_m when shifted
+ // left by 52 in scale_exp2:
+ int sign_shift = static_cast<int>(is_sinh) << 11;
+ int exp_diff = -2 * k - sign_shift;
+ double m_m = (k < 27) ? scale_exp2(1.0, exp_diff) : 0.0;
+
+ // Reduced argument:
+ // u = (x_a + kd * MLOG_2_HI) + kd * MLOG_2_LO.
+ // Since MLOG_2_HI has 35 significant bits and |kd| <= 1025,
+ // t1 = fma(kd, MLOG_2_HI, x_a)
+ // is exact in all rounding modes with or without FMA.
+ double t1 = fputil::multiply_add(kd, MLOG_2_HI, x_a); // exact
+ double t2 = kd * MLOG_2_LO;
+ double u = t1 + t2;
+ double u2 = u * u;
+ double u4 = u2 * u2;
+
+ // Evaluate degree-4 even and odd minimax polynomials P_even(u^2) and
+ // P_odd(u^2) in parallel using Estrin's scheme:
+ // cosh(u) ~ 1 + u^2 * P_even(u^2)
+ // sinh(u) ~ t1 + (t2 + u^3 * P_odd(u^2))
+
+ // First Estrin layer: c_0 + u^2 * c_1 and c_2 + u^2 * c_3
+ double e01 =
+ fputil::multiply_add(u2, EXP_EVEN_N0_COEFFS[1], EXP_EVEN_N0_COEFFS[0]);
+ double o01 =
+ fputil::multiply_add(u2, EXP_ODD_N0_COEFFS[1], EXP_ODD_N0_COEFFS[0]);
+ double e23 =
+ fputil::multiply_add(u2, EXP_EVEN_N0_COEFFS[3], EXP_EVEN_N0_COEFFS[2]);
+ double o23 =
+ fputil::multiply_add(u2, EXP_ODD_N0_COEFFS[3], EXP_ODD_N0_COEFFS[2]);
+
+ // Second Estrin layer: (c_2 + u^2 * c_3) + u^4 * c_4
+ double e24 = fputil::multiply_add(u4, EXP_EVEN_N0_COEFFS[4], e23);
+ double o24 = fputil::multiply_add(u4, EXP_ODD_N0_COEFFS[4], o23);
+
+ // Final polynomial values P_even(u^2) and P_odd(u^2):
+ double p_even = fputil::multiply_add(u4, e24, e01);
+ double p_odd = fputil::multiply_add(u4, o24, o01);
+
+ // Since m_m is either 0 or +/- 2^(-2 * k) with 0 <= k <= 26, both s and d are
+ // exact in all rounding modes:
+ double s = 1.0 + m_m; // exact
+ double d = 1.0 - m_m; // exact
+
+ // Apply sign(x) to the exact weights s and d before accumulating the
+ // polynomial so every subsequent rounding operation rounds in the active
+ // rounding direction for negative x as well (for cosh, the caller passes
+ // x >= 0 so sign_x == +0.0):
+ double sign_x = fputil::fp_xor(x, x_a);
+ s = fputil::fp_xor(s, sign_x);
+ d = fputil::fp_xor(d, sign_x);
+
+ // Accumulate from smallest to largest magnitude before adding the leading
+ // constant s. Adding s last avoids a double rounding at the leading ULP
+ // position of s when |d * t1| (< 2^-1) and |s * u^2 * P_even(u^2)| (< 2^-4)
+ // are added:
+ double d_t2 = d * t2;
+ double s_u2 = s * u2;
+ double d_u = d * u;
+ double d_u3 = d_u * u2;
+
+ // even_lo ~ s * u^2 * P_even(u^2) + d * t2
+ double even_lo = fputil::multiply_add(s_u2, p_even, d_t2);
+ // poly_lo ~ d * u^3 * P_odd(u^2) + even_lo
+ double poly_lo = fputil::multiply_add(d_u3, p_odd, even_lo);
+ // poly ~ d * t1 + poly_lo
+ double poly = fputil::multiply_add(d, t1, poly_lo);
+
+ // r ~ s * cosh(u) + d * sinh(u)
+ double r = s + poly;
+ return scale_exp2(r, k - 1);
+}
+
+#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// Fixed-point splitting constants (ulp = 2^-22, 2^-29) for exact multiplication
+// without Veltkamp splitting:
+LIBC_INLINE_VAR constexpr double SPLIT_22 = 0x1.8p30; // 1.5 * 2^(-22 + 52)
+LIBC_INLINE_VAR constexpr double SPLIT_29 = 0x1.8p23; // 1.5 * 2^(-29 + 52)
+
+// Evaluate:
+// sinh(x) = (e^x - e^-x) / 2 (when is_sinh is true)
+// cosh(x) = (e^x + e^-x) / 2 (when is_sinh is false, with x >= 0)
+// in double-double precision for 2^-2 <= |x| < 0x1.633ce8fb9f87ep+9 (sinh) or
+// 2^-27 < x < 0x1.633ce8fb9f87ep+9 (cosh), falling back to
+// exp_pm_accurate_128(x, is_sinh) when Ziv's rounding test fails.
+//
+// Range reduction:
+// Let |x| = k_p * (log(2) / 2^7) + u
+// in which:
+// k_p = round(|x| * 2^7 * log2(e)), k_m = -k_p
+// hi_p = k_p >> 7, idx_p = k_p & 0x7f
+// hi_m = k_m >> 7, idx_m = k_m & 0x7f
+// |u| <= 2^-8 * log(2)
+// then:
+// e^|x| = 2^hi_p * 2^(idx_p / 2^7) * e^u
+// e^-|x| = 2^hi_m * 2^(idx_m / 2^7) * e^-u
+// = 2^hi_p * (2^(hi_m - hi_p) * 2^(idx_m / 2^7)) * e^-u.
+//
+// Factoring out 2^(hi_p - 1):
+// (e^|x| -/+ e^-|x|) / 2 = 2^(hi_p - 1) * (m_p * e^u + m_m * e^-u)
+// where:
+// m_p = 2^(idx_p / 2^7)
+// m_m = (-/+ 1) * 2^(hi_m - hi_p) * 2^(idx_m / 2^7).
+//
+// For |x| < 0x1.633ce8fb9f87ep+9, hi_p is in [0, 1025] so hi_p - 1 is in
+// [-1, 1024]. When hi_p == 1025, idx_p == 0 and u < 0, which keeps:
+// r = m_p * e^u + m_m * e^-u < 1.0
+// (biased exponent 1022). Adding hi_p - 1 = 1024 in scale_exp2(r, hi_p - 1)
+// produces biased exponent 2046, handling the top octave
+// [1024 * log(2), 1025 * log(2)) without a separate branch.
+LIBC_MATH_SMALL_TABLES_INLINE double exp_pm_eval(double x, bool is_sinh) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ double x_a = xbits.abs().get_val();
+
+ constexpr double LOG2_E_SCALE = EXP_SCALE * LOG2_E.hi;
+ // Using ROUND_BIAS, ROUND_SHIFT, and ROUND_MASK (see exp_utils.h), we compute
+ // k_biased = fma(x_a, LOG2_E_SCALE, ROUND_BIAS) and extract in parallel:
+ // kd = (k_biased & ROUND_MASK) - ROUND_SHIFT
+ // k_p = (k_bits >> 32) & 0x3'ffff.
+ double k_biased = fputil::multiply_add(x_a, LOG2_E_SCALE, ROUND_BIAS);
+ double kd = fputil::fp_and(k_biased, ROUND_MASK) - ROUND_SHIFT;
+ int k_p =
+ static_cast<int>(cpp::bit_cast<uint64_t>(k_biased) >> 32) & 0x3'ffff;
+ int k_m = -k_p;
+
+ int hi_p = k_p >> EXP_N_BITS;
+ unsigned idx_p = static_cast<unsigned>(k_p & EXP_MASK);
+
+ int hi_m = k_m >> EXP_N_BITS;
+ unsigned idx_m = static_cast<unsigned>(k_m & EXP_MASK);
+
+ // Secondary exponent offset:
+ // exp_diff = hi_m - hi_p
+ // clamped at -100 (where e^-|x| <= 2^-99 * e^|x| << SINH_ERR_D) to avoid
+ // underflow in m_m:
+ int exp_diff = hi_m - hi_p;
+ int clamped_diff = (exp_diff > -100) ? exp_diff : -100;
+
+ // Subtracting (1 << 11) = 0x800 (for sinh) sets bit 11 of the 12-bit
+ // two's-complement exponent field, flipping the sign bit of m_m when shifted
+ // left by 52 in scale_exp2:
+ int sign_shift = static_cast<int>(is_sinh) << 11;
+ int scale_m = clamped_diff - sign_shift;
+
+ DoubleDouble m_p = EXP2_MID[idx_p];
+ // Scale m_m = EXP2_MID[idx_m] by (-/+ 1) * 2^clamped_diff:
+ // - Since m_m.hi is in [1, 2), we use scale_exp2(m_m.hi, scale_m) directly.
+ // - Since EXP2_MID[0].lo == 0.0, we scale m_m.lo by multiplying by:
+ // s_m = scale_exp2(1.0, scale_m).
+ double s_m = scale_exp2(1.0, scale_m);
+ DoubleDouble m_m = EXP2_MID[idx_m];
+ m_m.lo *= s_m;
+ m_m.hi = scale_exp2(m_m.hi, scale_m);
+
+ // Reduced argument:
+ // u = (x_a + kd * MLOG_2_EXP2_HI) + kd * MLOG_2_EXP2_MID.
+ // Since MLOG_2_EXP2_HI has 35 significant bits and |kd| < 2^18,
+ // t1 = fma(kd, MLOG_2_EXP2_HI, x_a)
+ // is exact in all rounding modes with or without FMA.
+ double t1 = fputil::multiply_add(kd, MLOG_2_EXP2_HI, x_a); // exact
+ double u = fputil::multiply_add(kd, MLOG_2_EXP2_MID, t1);
+ double u2 = u * u;
+
+ // Separate e^(+/-u) into even and odd parts:
+ // e^u ~ (1 + u^2 * P_even(u^2)) + (u + u^3 * P_odd(u^2))
+ // e^-u ~ (1 + u^2 * P_even(u^2)) - (u + u^3 * P_odd(u^2))
+ // where:
+ // P_even(u^2) = p_0 + p_2 * u^2
+ // P_odd(u^2) = p_1 + p_3 * u^2
+ // Then:
+ // m_p * e^u + m_m * e^-u
+ // ~ (m_p + m_m) * (1 + u^2 * P_even(u^2))
+ // + (m_p - m_m) * (u + u^3 * P_odd(u^2)).
+ double p_even = fputil::multiply_add(u2, EXP_P_COEFFS[2], EXP_P_COEFFS[0]);
+ double p_odd = fputil::multiply_add(u2, EXP_P_COEFFS[3], EXP_P_COEFFS[1]);
+
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ // Decompose m_p * e^u + m_m * e^-u ~ r_hi + r_lo:
+ // s.hi + s.lo = m_p.hi + m_m.hi (exact FastTwoSum)
+ // d.hi + d.lo = m_p.hi - m_m.hi (exact FastTwoSum)
+ // r_hi + r_err = s.hi + d.hi * t1
+ // r_lo ~ r_err + (s.lo + (m_p.lo - m_m.lo + d.lo) * t1)
+ // + ((m_p.lo + m_m.lo) + s.hi * u^2 * P_even(u^2))
+ // + d.hi * (kd * MLOG_2_EXP2_MID + u^3 * P_odd(u^2)).
+ //
+ // For |x| >= 2^-2, k_p >= 46 so hi_p >= 0 and hi_m <= -1, hence
+ // m_p.hi >= 1.0 > |m_m.hi|. By Sterbenz's lemma, both (s.hi - m_p.hi) and
+ // (m_p.hi - d.hi) are exact in all rounding modes.
+ DoubleDouble s = fputil::exact_add(m_p.hi, m_m.hi);
+ DoubleDouble d = fputil::exact_sub(m_p.hi, m_m.hi);
+
+ // Leading term r_hi = round(s.hi + d.hi * t1) and its rounding error r_err:
+ // Since |d.hi * t1| <= 2^-7 * log(2) < 2^-7 and s.hi > 2^-2, (s.hi - r_hi) is
+ // exact by Sterbenz's lemma.
+ double r_hi = fputil::multiply_add(d.hi, t1, s.hi);
+ double r_t = s.hi - r_hi; // exact
+ double r_err = fputil::multiply_add(d.hi, t1, r_t);
+
+ // Odd polynomial and low reduced-argument correction:
+ // odd_lo ~ (u - t1) + u^3 * P_odd(u^2) = kd * MLOG_2_EXP2_MID + u^3 * p_odd
+ double t2 = kd * MLOG_2_EXP2_MID;
+ double u3 = u * u2;
+ double odd_lo = fputil::multiply_add(u3, p_odd, t2);
+
+ // Combine table low parts and polynomial terms in parallel with the
+ // FastTwoSum residuals (s.lo, d.lo) and r_hi rounding error (r_err):
+ // mp_lo ~ m_p.lo * (1 + t1)
+ // mm_lo ~ m_m.lo * (1 - t1)
+ // even_lo ~ mp_lo + (s.hi * u^2) * P_even(u^2)
+ // poly_lo ~ even_lo + d.hi * odd_lo
+ double mp_lo = fputil::multiply_add(m_p.lo, t1, m_p.lo);
+ double mm_lo = fputil::multiply_add(-m_m.lo, t1, m_m.lo);
+ double s_u2 = s.hi * u2;
+ double even_lo = fputil::multiply_add(s_u2, p_even, mp_lo);
+ double poly_lo = fputil::multiply_add(d.hi, odd_lo, even_lo);
+
+ // Total linear error from FastTwoSum residuals and r_hi rounding:
+ // sd_lo ~ s.lo + d.lo * t1
+ // err_tot ~ r_err + sd_lo
+ double sd_lo = fputil::multiply_add(d.lo, t1, s.lo);
+ double err_tot = r_err + sd_lo;
+
+ // Apply sign(x) to r_hi, mm_lo, poly_lo, and err_tot before the final Ziv
+ // test additions so that (r_hi + lo_u) and (r_hi + lo_l) round in the active
+ // rounding direction for negative inputs (for cosh, the caller passes x >= 0
+ // so sign_x == +0.0):
+ double sign_x = fputil::fp_xor(x, x_a);
+ r_hi = fputil::fp_xor(r_hi, sign_x);
+ mm_lo = fputil::fp_xor(mm_lo, sign_x);
+ poly_lo = fputil::fp_xor(poly_lo, sign_x);
+ err_tot = fputil::fp_xor(err_tot, sign_x);
+
+ // Fold +/-SINH_ERR_D into mm_lo off the critical path so tail_{u,l} and
+ // err_tot finish in parallel:
+ constexpr double SINH_ERR_D = 0x1.4p-64;
+ double mm_lo_u = mm_lo + SINH_ERR_D;
+ double mm_lo_l = mm_lo - SINH_ERR_D;
+ double tail_u = poly_lo + mm_lo_u;
+ double tail_l = poly_lo + mm_lo_l;
+
+ double lo_u = err_tot + tail_u;
+ double lo_l = err_tot + tail_l;
+#else // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ // Fixed-point grid splitting for exact leading reconstruction in all rounding
+ // modes without FMA:
+ // lsb(mp_hi), lsb(mm_hi) >= 2^-22 => lsb(s_hi), lsb(d_hi) >= 2^-22
+ // lsb(u_hi) >= 2^-29
+ // Since s_hi = mp_hi + mm_hi > 2^-2, both s_hi, d_hi, and
+ // r_hi = s_hi + d_hi * u_hi (with lsb >= 2^-51)
+ // are exact in double.
+
+ // Split u = u_hi + u_lo with lsb(u_hi) >= 2^-29 and |u_lo| <= 2^-30:
+ double u_hi_t = t1 + SPLIT_29;
+ double u_hi = u_hi_t - SPLIT_29; // exact, lsb = 2^-29
+ double u_lo_t = t1 - u_hi; // exact
+ double u_lo = fputil::multiply_add(kd, MLOG_2_EXP2_MID, u_lo_t);
+
+ // Split m_p = mp_hi + mp_lo and m_m = mm_hi + mm_lo with lsb >= 2^-22:
+ double mp_hi_t = m_p.hi + SPLIT_22;
+ double mp_hi = mp_hi_t - SPLIT_22; // exact, lsb = 2^-22
+ double mp_lo_t = m_p.hi - mp_hi; // exact
+ double mp_lo = mp_lo_t + m_p.lo;
+
+ double mm_hi_t = m_m.hi + SPLIT_22;
+ double mm_hi = mm_hi_t - SPLIT_22; // exact, lsb = 2^-22
+ double mm_lo_t = m_m.hi - mm_hi; // exact
+ double mm_lo = mm_lo_t + m_m.lo;
+
+ // Exact leading sum, difference, and linear term r_hi = s_hi + d_hi * u_hi:
+ double s_hi = mp_hi + mm_hi; // exact, lsb = 2^-22
+ double d_hi = mp_hi - mm_hi; // exact, lsb = 2^-22
+ double r_hi = fputil::multiply_add(d_hi, u_hi, s_hi); // exact, lsb = 2^-51
+
+ // Low parts of the split table values and their linear contribution:
+ // s_lo = mp_lo + mm_lo, d_lo = mp_lo - mm_lo
+ // m_lo ~ s_lo + d_lo * u_hi
+ double s_lo = mp_lo + mm_lo;
+ double d_lo = mp_lo - mm_lo;
+ double m_lo = fputil::multiply_add(d_lo, u_hi, s_lo);
+
+ // Full-precision leading weights for the quadratic and cubic terms:
+ double s_m_hi = m_p.hi + m_m.hi;
+ double d_m_hi = m_p.hi - m_m.hi;
+
+ // Apply sign(x) before the final polynomial accumulation so that (r_hi +
+ // lo_u) and (r_hi + lo_l) round in the active rounding direction for negative
+ // x:
+ double sign_x = fputil::fp_xor(x, x_a);
+ r_hi = fputil::fp_xor(r_hi, sign_x);
+ s_m_hi = fputil::fp_xor(s_m_hi, sign_x);
+ d_m_hi = fputil::fp_xor(d_m_hi, sign_x);
+ m_lo = fputil::fp_xor(m_lo, sign_x);
+
+ // Combine even and odd polynomial terms with m_lo:
+ // even_lo ~ m_lo + s_m_hi * u^2 * P_even(u^2)
+ // odd_lo ~ u_lo + u^3 * P_odd(u^2)
+ // r_lo ~ even_lo + d_m_hi * odd_lo
+ double u2_peven = u2 * p_even;
+ double even_lo = fputil::multiply_add(s_m_hi, u2_peven, m_lo);
+ double u3 = u * u2;
+ double odd_lo = fputil::multiply_add(u3, p_odd, u_lo);
+ double r_lo = fputil::multiply_add(d_m_hi, odd_lo, even_lo);
+
+ // Total approximation and rounding error is bounded by SINH_ERR_D =
+ // 0x1.4p-64:
+ constexpr double SINH_ERR_D = 0x1.4p-64;
+ double lo_u = r_lo + SINH_ERR_D;
+ double lo_l = r_lo - SINH_ERR_D;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+
+ // 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))
+ return scale_exp2(upper, hi_p - 1);
+
+ return exp_pm_accurate_128(x, is_sinh);
+}
+
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+} // namespace double_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_SINHCOSH_UTILS_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 13cb65ff37e0b5..6b025c4fa60620 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -3753,6 +3753,16 @@ add_entrypoint_object(
libc.src.__support.math.coshf16
)
+add_entrypoint_object(
+ sinh
+ SRCS
+ sinh.cpp
+ HDRS
+ ../sinh.h
+ DEPENDS
+ libc.src.__support.math.sinh
+)
+
add_entrypoint_object(
sinhf
SRCS
diff --git a/libc/src/math/generic/sinh.cpp b/libc/src/math/generic/sinh.cpp
new file mode 100644
index 00000000000000..0f63d3191918e5
--- /dev/null
+++ b/libc/src/math/generic/sinh.cpp
@@ -0,0 +1,21 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 sinh(x) function.
+///
+//===----------------------------------------------------------------------===//
+
+#include "src/math/sinh.h"
+#include "src/__support/math/sinh.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(double, sinh, (double x)) { return math::sinh(x); }
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index cde49cd6e3c73d..80c3eb1e1a8f20 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -473,6 +473,7 @@ else()
libc.src.__support.math.sincosf
libc.src.__support.math.sinf
libc.src.__support.math.sinf16
+ libc.src.__support.math.sinh
libc.src.__support.math.sinhf
libc.src.__support.math.sinhf16
libc.src.__support.math.sqrtf128
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index d0581b5dababa9..1d70b4d3578203 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -351,6 +351,7 @@ TEST(LlvmLibcSharedMathTest, AllDouble) {
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::log2(1.0));
EXPECT_FP_EQ(1.0, LIBC_NAMESPACE::shared::pow(0.0, 0.0));
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::sin(0.0));
+ EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::sinh(0.0));
EXPECT_FP_EQ(1.0, cos);
EXPECT_FP_EQ(0.0, sin);
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::sqrt(0.0));
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index d88afe6370e27b..394f347782807b 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2595,6 +2595,22 @@ add_fp_unittest(
libc.src.math.coshf16
)
+add_fp_unittest(
+ sinh_test
+ NEED_MPFR
+ SUITE
+ libc-math-tests
+ SRCS
+ sinh_test.cpp
+ DEPENDS
+ libc.hdr.fenv_macros
+ libc.src.math.sinh
+ libc.src.__support.math.sinh_double_eval
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.optimization
+)
+
add_fp_unittest(
sinhf_test
NEED_MPFR
diff --git a/libc/test/src/math/sinh_test.cpp b/libc/test/src/math/sinh_test.cpp
new file mode 100644
index 00000000000000..5cada6e299c8b0
--- /dev/null
+++ b/libc/test/src/math/sinh_test.cpp
@@ -0,0 +1,187 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 sinh.
+///
+//===----------------------------------------------------------------------===//
+
+#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/sinh_double_eval.h"
+#include "src/math/sinh.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#ifndef LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+#define TOLERANCE 2
+#else // LIBC_MATH_HAS_ASSUME_ROUND_NEAREST_ONLY
+#define TOLERANCE 1
+#endif // !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;
+
+LIBC_NOINLINE double sinh_double_eval(double x) {
+ return LIBC_NAMESPACE::math::double_eval::sinh(x);
+}
+
+class SinhTest : 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(neg_inf, func(neg_inf));
+ EXPECT_FP_EQ_ALL_ROUNDING(0.0, func(0.0));
+ EXPECT_FP_EQ_ALL_ROUNDING(-0.0, func(-0.0));
+ }
+
+ void test_overflow(double (*func)(double)) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, func(0x1.0p20), FE_OVERFLOW);
+ EXPECT_FP_EQ_WITH_EXCEPTION(neg_inf, func(-0x1.0p20), FE_OVERFLOW);
+ }
+
+ void test_tricky_inputs(double (*func)(double), double ulp_tolerance = 0.5,
+ bool all_rounding = true) {
+ constexpr int N = 16;
+ constexpr uint64_t INPUTS[N] = {
+ 0x3e50'0000'0000'0000ULL, // x = 0x1.0p-26
+ 0x3e50'0000'0000'0001ULL, // x = 0x1.0000000000001p-26
+ 0x3e6b'cee7'0ebe'7ec9ULL, // x = 0x1.bcee70ebe7ec9p-25
+ 0x3f18'a62b'e5b5'f93cULL, // x = 0x1.8a62be5b5f93cp-14
+ 0x3f8f'ffff'ffff'ffffULL, // x = 0x1.fffffffffffffp-7
+ 0x3fb3'30b5'a711'0cf0ULL, // x = 0x1.330b5a7110cfp-4
+ 0x3fcf'ffff'ffff'ffffULL, // x = 0x1.fffffffffffffp-3
+ 0x3fd0'0000'0000'0000ULL, // x = 0x1.0p-2
+ 0x3fd3'b6db'f12c'e919ULL, // x = 0x1.3b6dbf12ce919p-2
+ 0x3fe6'2e42'fefa'39efULL, // x = 0x1.62e42fefa39efp-1
+ 0x3ff7'4385'446d'71c3ULL, // x = 0x1.74385446d71c3p+0
+ 0x401c'ac90'6a4f'1a46ULL, // x = 0x1.cac906a4f1a46p+2
+ 0x4032'57b9'12e2'b431ULL, // x = 0x1.257b912e2b431p+4
+ 0x4086'2e42'fefa'39efULL, // x = 0x1.62e42fefa39efp+9
+ 0x4086'33ce'8fb9'f87dULL, // x = 0x1.633ce8fb9f87dp+9
+ 0x4086'33ce'8fb9'f87eULL, // x = 0x1.633ce8fb9f87ep+9
+ };
+ for (int i = 0; i < N; ++i) {
+ double x = FPBits(INPUTS[i]).get_val();
+ if (all_rounding) {
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinh, x, func(x),
+ ulp_tolerance);
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinh, -x, func(-x),
+ ulp_tolerance);
+ } else {
+ EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, x, func(x), ulp_tolerance);
+ EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, -x, func(-x), ulp_tolerance);
+ }
+ }
+ }
+
+ 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(0x1.0p-30).uintval();
+ constexpr uint64_t STOP = FPBits(710.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::Sinh, x, result,
+ ulp_tolerance, rounding_mode)) {
+ ++fails;
+ while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(
+ mpfr::Operation::Sinh, x, result, tol, rounding_mode)) {
+ mx = x;
+ mr = result;
+
+ if (tol > 1000.0)
+ break;
+
+ tol *= 2.0;
+ }
+ }
+ }
+ if (fails) {
+ tlog << " Sinh failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol)
+ << ".\n";
+ EXPECT_MPFR_MATCH(mpfr::Operation::Sinh, mx, mr, ulp_tolerance,
+ rounding_mode);
+ }
+ };
+
+ tlog << " Test Rounding To Nearest...\n";
+ test(mpfr::RoundingMode::Nearest);
+
+ if (all_rounding) {
+ tlog << " Test Rounding Downward...\n";
+ test(mpfr::RoundingMode::Downward);
+
+ tlog << " Test Rounding Upward...\n";
+ test(mpfr::RoundingMode::Upward);
+
+ tlog << " Test Rounding Toward Zero...\n";
+ test(mpfr::RoundingMode::TowardZero);
+ }
+ }
+};
+
+#define LIST_SINH_TESTS(suffix, func, ulp_tolerance, all_rounding) \
+ using LlvmLibcSinhTest##suffix = SinhTest; \
+ TEST_F(LlvmLibcSinhTest##suffix, SpecialNumbers) { \
+ test_special_numbers(&func); \
+ } \
+ TEST_F(LlvmLibcSinhTest##suffix, Overflow) { test_overflow(&func); } \
+ TEST_F(LlvmLibcSinhTest##suffix, TrickyInputs) { \
+ test_tricky_inputs(&func, ulp_tolerance, all_rounding); \
+ } \
+ TEST_F(LlvmLibcSinhTest##suffix, InDoubleRange) { \
+ test_in_range(&func, ulp_tolerance, all_rounding); \
+ } \
+ static_assert(true, "Require semicolon.")
+
+LIST_SINH_TESTS(Default, LIBC_NAMESPACE::sinh,
+ /*ulp_tolerance=*/TOLERANCE + 0.5,
+ /*all_rounding=*/ALL_ROUNDING);
+LIST_SINH_TESTS(DoubleEval, sinh_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 76caeef863e342..448add6cf8a73f 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -4758,6 +4758,21 @@ add_fp_unittest(
libc.src.__support.FPUtil.cast
)
+add_fp_unittest(
+ sinh_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ sinh_test.cpp
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.math.sinh
+ libc.src.__support.math.sinh_double_eval
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.macros.attributes
+)
+
add_fp_unittest(
sinhf_test
SUITE
diff --git a/libc/test/src/math/smoke/sinh_test.cpp b/libc/test/src/math/smoke/sinh_test.cpp
new file mode 100644
index 00000000000000..e7f8439973ff96
--- /dev/null
+++ b/libc/test/src/math/smoke/sinh_test.cpp
@@ -0,0 +1,123 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 sinh.
+///
+//===----------------------------------------------------------------------===//
+
+#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/sinh_double_eval.h"
+#include "src/math/sinh.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+LIBC_NOINLINE double sinh_double_eval(double x) {
+ return LIBC_NAMESPACE::math::double_eval::sinh(x);
+}
+
+class SinhTest : 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_ALL_ROUNDING(aNaN, func(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, func(inf));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(neg_inf, func(neg_inf));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(0.0, func(0.0));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(-0.0, func(-0.0));
+ EXPECT_MATH_ERRNO(0);
+ }
+
+ void test_overflow(double (*func)(double)) {
+ // x = 0x1.fffffffffffffp+1023 (max_normal)
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf, func(FPBits(uint64_t(0x7fef'ffff'ffff'ffff)).get_val()),
+ FE_OVERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ // x = 0x1.633ce8fb9f87ep+9 (smallest positive input that overflows)
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf, func(FPBits(uint64_t(0x4086'33ce'8fb9'f87e)).get_val()),
+ FE_OVERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ // x = -0x1.633ce8fb9f87ep+9 (largest negative input that overflows)
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ neg_inf, func(FPBits(uint64_t(0xc086'33ce'8fb9'f87e)).get_val()),
+ FE_OVERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+ }
+
+ void test_underflow(double (*func)(double)) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(min_denormal, func(min_denormal), FE_UNDERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(-min_denormal, func(-min_denormal),
+ FE_UNDERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+ }
+
+#ifdef LIBC_TEST_FTZ_DAZ
+ void test_denormals(double (*func)(double)) {
+ EXPECT_FP_EQ(0.0, func(min_denormal));
+ EXPECT_FP_EQ(0.0, func(max_denormal));
+ }
+#endif // LIBC_TEST_FTZ_DAZ
+};
+
+#ifdef LIBC_TEST_FTZ_DAZ
+#define LIST_SINH_FTZ_DAZ_TESTS(suffix, func) \
+ TEST_F(LlvmLibcSinhTest##suffix, FTZMode) { \
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::FTZ); \
+ test_denormals(&func); \
+ } \
+ TEST_F(LlvmLibcSinhTest##suffix, DAZMode) { \
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::DAZ); \
+ test_denormals(&func); \
+ } \
+ TEST_F(LlvmLibcSinhTest##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_SINH_FTZ_DAZ_TESTS(suffix, func)
+#endif // LIBC_TEST_FTZ_DAZ
+
+#define LIST_SINH_TESTS(suffix, func, check_snan_invalid) \
+ using LlvmLibcSinhTest##suffix = SinhTest; \
+ TEST_F(LlvmLibcSinhTest##suffix, SpecialNumbers) { \
+ test_special_numbers(&func, check_snan_invalid); \
+ } \
+ TEST_F(LlvmLibcSinhTest##suffix, Overflow) { test_overflow(&func); } \
+ TEST_F(LlvmLibcSinhTest##suffix, Underflow) { test_underflow(&func); } \
+ LIST_SINH_FTZ_DAZ_TESTS(suffix, func) \
+ static_assert(true, "Require semicolon.")
+
+LIST_SINH_TESTS(Default, LIBC_NAMESPACE::sinh, /*check_snan_invalid=*/true);
+LIST_SINH_TESTS(DoubleEval, sinh_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 210e2e1eb4a09a..d967412e293f40 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -10762,10 +10762,11 @@ libc_support_library(
hdrs = ["src/__support/math/exp2_double_eval.h"],
deps = [
":__support_common",
+ ":__support_cpp_bit",
":__support_fputil_double_double",
":__support_fputil_fp_bits",
+ ":__support_fputil_fp_bitwise",
":__support_fputil_multiply_add",
- ":__support_fputil_nearest_integer",
":__support_macros_attributes",
":__support_macros_config",
":__support_macros_optimization",
@@ -10791,10 +10792,11 @@ libc_support_library(
hdrs = ["src/__support/math/exp_double_eval.h"],
deps = [
":__support_common",
+ ":__support_cpp_bit",
":__support_fputil_double_double",
":__support_fputil_fp_bits",
+ ":__support_fputil_fp_bitwise",
":__support_fputil_multiply_add",
- ":__support_fputil_nearest_integer",
":__support_macros_attributes",
":__support_macros_config",
":__support_macros_optimization",
@@ -11523,6 +11525,75 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_sinhcosh_accurate_128",
+ hdrs = ["src/__support/math/sinhcosh_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_math_exp2_accurate_128",
+ ":__support_uint128",
+ ],
+)
+
+libc_support_library(
+ name = "__support_math_sinhcosh_utils",
+ hdrs = ["src/__support/math/sinhcosh_utils.h"],
+ deps = [
+ ":__support_common",
+ ":__support_cpp_bit",
+ ":__support_fputil_double_double",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_fp_bitwise",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_nearest_integer",
+ ":__support_macros_attributes",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_macros_properties_cpu_features",
+ ":__support_math_exp_utils",
+ ":__support_math_sinhcosh_accurate_128",
+ ],
+)
+
+libc_support_library(
+ name = "__support_math_sinh_double_eval",
+ hdrs = ["src/__support/math/sinh_double_eval.h"],
+ deps = [
+ ":__support_common",
+ ":__support_cpp_bit",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_fp_bitwise",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_rounding_mode",
+ ":__support_macros_attributes",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_macros_properties_cpu_features",
+ ":__support_math_sinhcosh_utils",
+ ":hdr_errno_macros",
+ ":hdr_fenv_macros",
+ ],
+)
+
+libc_support_library(
+ name = "__support_math_sinh",
+ hdrs = ["src/__support/math/sinh.h"],
+ deps = [
+ ":__support_common",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_math_sinh_double_eval",
+ ],
+)
+
libc_support_library(
name = "__support_math_sinhf",
hdrs = ["src/__support/math/sinhf.h"],
@@ -15409,6 +15480,11 @@ libc_math_function(
additional_deps = [":__support_math_sincosf"],
)
+libc_math_function(
+ name = "sinh",
+ additional_deps = [":__support_math_sinh"],
+)
+
libc_math_function(
name = "sinhf",
additional_deps = [":__support_math_sinhf"],
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 c39dc612851536..89289849d4ba73 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
@@ -936,6 +936,14 @@ math_mpfr_test(
],
)
+math_mpfr_test(
+ name = "sinh",
+ deps = [
+ "//libc:__support_macros_attributes",
+ "//libc:__support_math_sinh_double_eval",
+ ],
+)
+
math_mpfr_test(
name = "sinhf",
deps = [
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 238b2852e4cc6c..e5bd15412d0611 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
@@ -2148,6 +2148,14 @@ math_test(name = "sincos")
math_test(name = "sincosf")
+math_test(
+ name = "sinh",
+ deps = [
+ "//libc:__support_macros_attributes",
+ "//libc:__support_math_sinh_double_eval",
+ ],
+)
+
math_test(
name = "sinhf",
deps = [
>From 3e903d2b64ae5b07f14a3917025ccca7a794a33f Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Fri, 23 Oct 2026 21:09:58 +0000
Subject: [PATCH 2/3] Also use rounding trick for small table version.
---
libc/src/__support/math/CMakeLists.txt | 1 -
libc/src/__support/math/sinhcosh_utils.h | 11 +++++++----
utils/bazel/llvm-project-overlay/libc/BUILD.bazel | 1 -
3 files changed, 7 insertions(+), 6 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index f55bcf412bbffd..29b085e1d0c163 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -5765,7 +5765,6 @@ add_header_library(
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
diff --git a/libc/src/__support/math/sinhcosh_utils.h b/libc/src/__support/math/sinhcosh_utils.h
index f93d5b56a9d181..2782095cc5fd3d 100644
--- a/libc/src/__support/math/sinhcosh_utils.h
+++ b/libc/src/__support/math/sinhcosh_utils.h
@@ -19,7 +19,6 @@
#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"
@@ -126,9 +125,13 @@ LIBC_MATH_SMALL_TABLES_INLINE double exp_pm_eval(double x, bool is_sinh) {
FPBits xbits(x);
double x_a = xbits.abs().get_val();
- double x_scaled = x_a * 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_a, LOG2_E.hi, ROUND_BIAS) and extract in parallel:
+ // kd = (k_biased & ROUND_MASK) - ROUND_SHIFT
+ // k = (k_bits >> 32) & 0x3'ffff.
+ double k_biased = fputil::multiply_add(x_a, LOG2_E.hi, ROUND_BIAS);
+ double kd = fputil::fp_and(k_biased, ROUND_MASK) - ROUND_SHIFT;
+ int k = static_cast<int>(cpp::bit_cast<uint64_t>(k_biased) >> 32) & 0x3'ffff;
// Secondary scale factor:
// m_m = (-/+ 1) * 2^(-2 * k) for k < 27
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index d967412e293f40..299993bbebba70 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -11552,7 +11552,6 @@ libc_support_library(
":__support_fputil_fp_bits",
":__support_fputil_fp_bitwise",
":__support_fputil_multiply_add",
- ":__support_fputil_nearest_integer",
":__support_macros_attributes",
":__support_macros_config",
":__support_macros_optimization",
>From 4748fc2e52c0b348a2ffb82ccdb93488c211717f Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Mon, 26 Oct 2026 02:16:02 +0000
Subject: [PATCH 3/3] Use default initialization for NumberPair<T> in
double_double.h.
---
libc/src/__support/FPUtil/double_double.h | 10 +++++-----
1 file changed, 5 insertions(+), 5 deletions(-)
diff --git a/libc/src/__support/FPUtil/double_double.h b/libc/src/__support/FPUtil/double_double.h
index 4d93c4769b1fb2..45ce8a10acdac4 100644
--- a/libc/src/__support/FPUtil/double_double.h
+++ b/libc/src/__support/FPUtil/double_double.h
@@ -51,7 +51,7 @@ using FloatFloat = NumberPair<float>;
// <https://www.arith2025.org/proceedings/215900a141.pdf>.
template <bool FAST2SUM = true, typename T = double>
LIBC_INLINE constexpr NumberPair<T> exact_add(T a, T b) {
- NumberPair<T> r{0.0, 0.0};
+ NumberPair<T> r;
if constexpr (FAST2SUM) {
r.hi = a + b;
T t = r.hi - a;
@@ -69,7 +69,7 @@ LIBC_INLINE constexpr NumberPair<T> exact_add(T a, T b) {
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};
+ NumberPair<T> r;
if constexpr (FAST2SUM) {
r.hi = a - b;
T t = a - r.hi;
@@ -113,7 +113,7 @@ LIBC_INLINE constexpr NumberPair<T> add(const NumberPair<T> &a, T b) {
// Default splitting constant = 2^ceil(prec(double)/2) + 1 = 2^27 + 1.
template <typename T = double, size_t N = DefaultSplit<T>::VALUE>
LIBC_INLINE constexpr NumberPair<T> split(T a) {
- NumberPair<T> r{0.0, 0.0};
+ NumberPair<T> r;
// CN = 2^N.
constexpr T CN = static_cast<T>(1 << N);
constexpr T C = CN + T(1);
@@ -129,7 +129,7 @@ template <typename T = double, size_t SPLIT_B = DefaultSplit<T>::VALUE>
LIBC_INLINE constexpr NumberPair<T> exact_mult(const NumberPair<T> &as, T a,
T b) {
NumberPair<T> bs = split<T, SPLIT_B>(b);
- NumberPair<T> r{0.0, 0.0};
+ NumberPair<T> r;
r.hi = a * b;
T t1 = as.hi * bs.hi - r.hi;
@@ -169,7 +169,7 @@ template <> struct TargetHasFmaInstruction<double> {
// then a * b = r.hi + r.lo is exact for all rounding modes.
template <typename T = double, size_t SPLIT_B = DefaultSplit<T>::VALUE>
LIBC_INLINE LIBC_CONSTEXPR NumberPair<T> exact_mult(T a, T b) {
- NumberPair<T> r{0.0, 0.0};
+ NumberPair<T> r;
if constexpr (TargetHasFmaInstruction<T>::VALUE) {
r.hi = a * b;
More information about the libc-commits
mailing list