[llvm-branch-commits] [libc] [llvm] [libc][math] Implement double precision tanh correctly rounded to all rounding modes. (PR #235618)
via llvm-branch-commits
llvm-branch-commits at lists.llvm.org
Sun Oct 25 23:53:34 PDT 2026
https://github.com/lntue created https://github.com/llvm/llvm-project/pull/235618
- Evaluate `tanh(x) = 1 - (2 * m) / (e^u + m)` for `|x| >= 2^-3` using shared exp range reduction and lookup tables from `exp_utils.h` and `exp2_accurate_128.h`
- For `2^-27 < |x| < 2^-3`, use a degree-5 minimax polynomial in `x^2` with exact leading cubic split in the fast pass, and a degree-11 minimax polynomial in Frac128 in the accurate pass.
- Add a no-LUT (N = 0) sinh(x) / cosh(x) evaluation when `LIBC_MATH_HAS_SKIP_ACCURATE_PASS` is enabled.
Assisted-by: Gemini for analysis and refactoring.
>From 16bd677a6410199e8a1d0ef901ebba2166fbf5bb Mon Sep 17 00:00:00 2001
From: Tue Ly <lntue.h at gmail.com>
Date: Mon, 26 Oct 2026 06:14:32 +0000
Subject: [PATCH] [libc][math] Implement double precision tanh correctly
rounded to all rounding modes.
---
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/tanh.h | 28 +
libc/src/__support/math/CMakeLists.txt | 51 ++
libc/src/__support/math/tanh.h | 34 ++
libc/src/__support/math/tanh_accurate_128.h | 221 +++++++
libc/src/__support/math/tanh_double_eval.h | 559 ++++++++++++++++++
libc/src/math/generic/CMakeLists.txt | 10 +
libc/src/math/generic/tanh.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/smoke/CMakeLists.txt | 16 +
libc/test/src/math/smoke/tanh_test.cpp | 146 +++++
libc/test/src/math/tanh_test.cpp | 187 ++++++
.../llvm-project-overlay/libc/BUILD.bazel | 56 ++
.../libc/test/src/math/BUILD.bazel | 8 +
.../libc/test/src/math/smoke/BUILD.bazel | 9 +
31 files changed, 1384 insertions(+), 1 deletion(-)
create mode 100644 libc/shared/math/tanh.h
create mode 100644 libc/src/__support/math/tanh.h
create mode 100644 libc/src/__support/math/tanh_accurate_128.h
create mode 100644 libc/src/__support/math/tanh_double_eval.h
create mode 100644 libc/src/math/generic/tanh.cpp
create mode 100644 libc/test/src/math/smoke/tanh_test.cpp
create mode 100644 libc/test/src/math/tanh_test.cpp
diff --git a/libc/config/baremetal/aarch64/entrypoints.txt b/libc/config/baremetal/aarch64/entrypoints.txt
index 6ac24afcfbbcdc..ab95344a60e1df 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/entrypoints.txt
@@ -606,6 +606,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.totalorder
libc.src.math.totalorderf
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index a71e3fc8cb1c8a..ceaaca42e22b9f 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -618,6 +618,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.totalorder
libc.src.math.totalorderf
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index 1a481461473787..37bb2c7b4d6605 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -609,6 +609,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.totalorder
libc.src.math.totalorderf
diff --git a/libc/config/darwin/aarch64/entrypoints.txt b/libc/config/darwin/aarch64/entrypoints.txt
index 20ab7d6e180528..930d23f60c5f6d 100644
--- a/libc/config/darwin/aarch64/entrypoints.txt
+++ b/libc/config/darwin/aarch64/entrypoints.txt
@@ -427,6 +427,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.tanpif
libc.src.math.totalorder
diff --git a/libc/config/freebsd/x86_64/entrypoints.txt b/libc/config/freebsd/x86_64/entrypoints.txt
index 4c256333d9228b..8c243cae40f229 100644
--- a/libc/config/freebsd/x86_64/entrypoints.txt
+++ b/libc/config/freebsd/x86_64/entrypoints.txt
@@ -314,6 +314,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.trunc
libc.src.math.truncf
diff --git a/libc/config/gpu/amdgpu/entrypoints.txt b/libc/config/gpu/amdgpu/entrypoints.txt
index 880327cba1cbb0..777e451f179d78 100644
--- a/libc/config/gpu/amdgpu/entrypoints.txt
+++ b/libc/config/gpu/amdgpu/entrypoints.txt
@@ -544,6 +544,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.tanpif
libc.src.math.tgamma
diff --git a/libc/config/gpu/nvptx/entrypoints.txt b/libc/config/gpu/nvptx/entrypoints.txt
index 1c8a02af59a836..7b81ea166c2489 100644
--- a/libc/config/gpu/nvptx/entrypoints.txt
+++ b/libc/config/gpu/nvptx/entrypoints.txt
@@ -541,6 +541,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.tanpif
libc.src.math.tgamma
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index 79626a1ad727d7..c5490a86378404 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -844,6 +844,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.tanpif
libc.src.math.totalorder
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index 2d8bbf275517e9..b381d5128e83b9 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -620,6 +620,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.tanpif
libc.src.math.totalorder
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index a5f834f8e2c62d..c59040976df840 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -905,6 +905,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.tanpif
libc.src.math.totalorder
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 110dca3a73ab72..0b09a0da19941e 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -927,6 +927,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.tanpif
libc.src.math.totalorder
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index 4abd600c61d1dc..5e70ff6c08d6d1 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -347,6 +347,7 @@ set(TARGET_LIBM_ENTRYPOINTS
libc.src.math.sqrtl
libc.src.math.tan
libc.src.math.tanf
+ libc.src.math.tanh
libc.src.math.tanhf
libc.src.math.tanpif
libc.src.math.trunc
diff --git a/libc/docs/headers/math/index.md b/libc/docs/headers/math/index.md
index e0393075cd900b..c863dce0b22748 100644
--- a/libc/docs/headers/math/index.md
+++ b/libc/docs/headers/math/index.md
@@ -215,7 +215,7 @@ log
| 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 |
-| tanh | {{ check }} | | | {{ check }} | | | 7.12.5.6 | F.10.2.6 |
+| tanh | {{ check }} | {{ check }} | | {{ check }} | | | 7.12.5.6 | F.10.2.6 |
| tanpi | {{ check }} | | | {{ check }} | | | 7.12.4.14 | F.10.1.14 |
| tgamma | | | | | | | 7.12.8.4 | F.10.5.4 |
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index 09668300f9ea88..b4fe177793e7b1 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -2742,6 +2742,12 @@ functions:
arguments:
- type: _Float16
guard: LIBC_TYPES_HAS_FLOAT16
+ - name: tanh
+ standards:
+ - stdc
+ return_type: double
+ arguments:
+ - type: double
- name: tanhf
standards:
- stdc
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 1f676b8da2e4f8..2308cc0430f067 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -483,6 +483,7 @@
#include "math/tanbf16.h"
#include "math/tanf.h"
#include "math/tanf16.h"
+#include "math/tanh.h"
#include "math/tanhf.h"
#include "math/tanhf16.h"
#include "math/tanpif.h"
diff --git a/libc/shared/math/tanh.h b/libc/shared/math/tanh.h
new file mode 100644
index 00000000000000..69c0d71a7ac15c
--- /dev/null
+++ b/libc/shared/math/tanh.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 tanh function.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SHARED_MATH_TANH_H
+#define LLVM_LIBC_SHARED_MATH_TANH_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/tanh.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::tanh;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_TANH_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 390220bd5a4b88..a7119e631ec286 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -6026,6 +6026,57 @@ add_header_library(
libc.include.llvm-libc-macros.float16_macros
)
+add_header_library(
+ tanh_accurate_128
+ HDRS
+ tanh_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(
+ tanh_double_eval
+ HDRS
+ tanh_double_eval.h
+ DEPENDS
+ .exp_utils
+ .tanh_accurate_128
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.CPP.bit
+ libc.src.__support.FPUtil.double_double
+ libc.src.__support.FPUtil.except_value_utils
+ 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.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(
+ tanh
+ HDRS
+ tanh.h
+ DEPENDS
+ .tanh_double_eval
+ libc.src.__support.common
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+)
+
add_header_library(
tanhf
HDRS
diff --git a/libc/src/__support/math/tanh.h b/libc/src/__support/math/tanh.h
new file mode 100644
index 00000000000000..7f63daf38bb25b
--- /dev/null
+++ b/libc/src/__support/math/tanh.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 tanh.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_TANH_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_TANH_H
+
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+
+#include "src/__support/math/tanh_double_eval.h"
+#define LIBC_MATH_TANH_IMPL double_eval
+
+namespace LIBC_NAMESPACE_DECL {
+namespace math {
+
+using LIBC_MATH_TANH_IMPL::tanh;
+
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#undef LIBC_MATH_TANH_IMPL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_TANH_H
diff --git a/libc/src/__support/math/tanh_accurate_128.h b/libc/src/__support/math/tanh_accurate_128.h
new file mode 100644
index 00000000000000..de493c133cc226
--- /dev/null
+++ b/libc/src/__support/math/tanh_accurate_128.h
@@ -0,0 +1,221 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 tanh function.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_TANH_ACCURATE_128_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_TANH_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 (x - tanh(x)) / x^3 on [-2^-3, 2^-3] in terms
+// of w = 64 * x^2 in [0, 1] as a degree-12 alternating polynomial in 128-bit
+// fixed-point format (lsb = 2^-128):
+// P(w) = c_0 - c_1 * w + c_2 * w^2 - ... + c_12 * w^12
+// tanh(x) ~ x - x^3 * P(64 * x^2).
+//
+// For k >= 9, the scaled coefficients c_k 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 TANH_SMALL_W_LO64.
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > display = hexadecimal!;
+// > prec = 500;
+// > P = fpminimax((1 - tanh(2^-3 * sqrt(x)) / (2^-3 * sqrt(x))) / (2^-6 * x),
+// 12, [|128...|], [2^(-150), 1], fixed);
+// > for i from 0 to 12 do {
+// c = abs(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(tanh(x) - (x - x^3 * subst(P, x, 64 * x^2)), [-2^-3, 2^-3]);
+// 0x1.7978...p-131 < 2^-130.
+LIBC_INLINE_VAR constexpr Frac128 TANH_SMALL_W_128[9] = {
+ // lsb = 2^-128
+ Frac128({0x5555555555555501ULL, 0x5555555555555555ULL}), // c_0
+ Frac128({0x8888888888881951ULL, 0x0088888888888888ULL}), // c_1
+ Frac128({0xdd0dd0dd0db889b5ULL, 0x0000dd0dd0dd0dd0ULL}), // c_2
+ Frac128({0x0f9f32d6b98d67dfULL, 0x000001664f4882c1ULL}), // c_3
+ Frac128({0xd84791eb52525c45ULL, 0x0000000244dc6abcULL}), // c_4
+ Frac128({0xc2afbbf61d062422ULL, 0x0000000003ada7a1ULL}), // c_5
+ Frac128({0x914add75e7e579bbULL, 0x000000000005f68dULL}), // c_6
+ Frac128({0xc12401ab85fc9a4aULL, 0x00000000000009aaULL}), // c_7
+ Frac128({0xabebb9964eabe8bdULL, 0x000000000000000fULL}), // c_8
+};
+
+LIBC_INLINE_VAR constexpr Frac64 TANH_SMALL_W_LO64[4] = {
+ // lsb = 2^-128 (val[1] == 0)
+ Frac64(0x1967e174ab55fb50ULL), // c_9
+ Frac64(0x00292f9d250de4b0ULL), // c_10
+ Frac64(0x000042b85e8bcc87ULL), // c_11
+ Frac64(0x00000067e2f87d8eULL), // c_12
+};
+
+// 128-bit accurate path for tanh(x) when 2^-27 < |x| < 2^-3:
+LIBC_NOINLINE LIBC_INLINE double tanh_small_accurate_128(double x) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ int exp_x = xbits.get_exponent(); // in [-27, -4]
+
+ // 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 / 4 by 2^(2 * exp_x + 8) to obtain w = 64 * x^2 in (0, 1) with
+ // lsb = 2^-128. Since exp_x is in [-27, -4], shift_w = -2 * exp_x - 8 is in
+ // [0, 46]:
+ unsigned shift_w = static_cast<unsigned>(-2 * exp_x - 8); // in [0, 46]
+ Frac128 w = Frac128(m_x2.val) >> shift_w; // lsb = 2^-128
+
+ // High-degree alternating Horner steps (k = 12, 11, 10, 9 -> 8) using 64-bit
+ // multiplications since val[1] == 0 for c_9..c_12:
+ Frac64 w_hi = w.to_frac64(); // lsb = 2^-64
+ Frac64 p_lo = TANH_SMALL_W_LO64[3];
+ for (int i = 2; i >= 0; --i)
+ p_lo = TANH_SMALL_W_LO64[i] - p_lo * w_hi;
+ p_lo = Frac64(TANH_SMALL_W_128[8].val[0]) - p_lo * w_hi;
+
+ // Low-degree alternating Horner steps (k = 7..0) in 128-bit fixed-point:
+ Frac128 p({p_lo.val[0], TANH_SMALL_W_128[8].val[1]}); // lsb = 2^-128
+ LIBC_LOOP_NOUNROLL
+ for (int i = 7; i >= 0; --i)
+ p = TANH_SMALL_W_128[i] - p * w; // lsb = 2^-128
+
+ // Reconstruction of |tanh(x)| = |x| - |x|^3 * P(64 * x^2):
+ // Since y_poly = m_x3 * p represents m_x^3 * P(w) with lsb = 2^-125,
+ // shifting right by shift = -2 * exp_x - 2 (in [6, 52]) aligns:
+ // 2^(3 * exp_x) * m_x^3 * P(w)
+ // to the same scale as m_x_frac (lsb = 2^(exp_x - 127)):
+ unsigned shift = static_cast<unsigned>(-2 * exp_x - 2); // in [6, 52]
+ Frac128 y_poly = m_x3 * p; // lsb = 2^-125
+ Frac128 y_m = y_poly >> shift; // lsb = 2^(exp_x - 127)
+
+ // Subtract |x|^3 * P(w) from |x| and round to double in the current rounding
+ // mode (DFloat128 normalizes x_m if the leading bit shifts to bit 126):
+ UInt<128> x_m = UInt<128>(m_x_frac) - UInt<128>(y_m); // lsb = 2^(exp_x - 127)
+ DFloat128 r(xbits.sign(), exp_x - 127, x_m);
+ return static_cast<double>(r);
+}
+
+// Compute tanh(x) = (e^(2x) - 1) / (e^(2x) + 1) in 128-bit precision using
+// direct base-2 192-bit fixed-point range reduction for
+// 2^-3 <= |x| < 0x1.30fc1931f09cap+4:
+//
+// Writing z = 2 * |x| * log2(e) = hi_p + idx_p * 2^-4 + u with hi_p in [0, 55],
+// idx_p in [0, 15], and u in [0, 2^-4), we have:
+// e^(2 * |x|) = 2^hi_p * m_final
+// where m_final = exp2_eval_frac128(idx_p, u) = 2^(idx_p * 2^-4 + u) in [1, 2).
+// Dividing numerator and denominator by 2^hi_p gives:
+// tanh(|x|) = (m_final - 2^-hi_p) / (m_final + 2^-hi_p)
+// where 2^-hi_p is a single bit in Frac128.
+LIBC_NOINLINE LIBC_INLINE double tanh_accurate_128(double x) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ int exp_x = xbits.get_exponent(); // in [-3, 4]
+
+ // 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 = 2 * |x| * log2(e):
+ // 2 * |x| = 2^(exp_x - 62) * mant_x
+ // LOG2_E_192.log2_b has lsb = 2^-191
+ // p = mant_x * LOG2_E_192.log2_b has lsb = 2^(exp_x - 253).
+ UInt<256> p = LOG2_E_192.log2_b.ful_mul(UInt<64>(mant_x));
+
+ // Since p.val[3] has lsb = 2^(exp_x - 61), the 2^-4 bit of z is at bit
+ // shr = 57 - exp_x (in [53, 60]) of p.val[3]:
+ int shr = 57 - exp_x; // in [53, 60]
+ int shl = exp_x + 7; // in [4, 11], with shr + shl = 64
+ int k_p = static_cast<int>(p.val[3] >> shr); // floor(z * 2^4)
+ uint64_t u_lo = (p.val[2] << shl) | (p.val[1] >> shr);
+ uint64_t u_hi = (p.val[3] << shl) | (p.val[2] >> shr);
+ // Scale (u_lo, u_hi) from lsb = 2^-132 to u in [0, 2^-4) with lsb = 2^-128:
+ Frac128 u = Frac128({u_lo, u_hi}) >> 4; // lsb = 2^-128
+
+ // Decompose k_p = 16 * hi_p + idx_p and evaluate
+ // m_final = 2^(idx_p * 2^-4 + u) in [1, 2) with lsb = 2^-127:
+ int hi_p = k_p >> 4; // in [0, 55]
+ unsigned idx_p = static_cast<unsigned>(k_p & 0x0f);
+ Frac128 m_final = exp2_eval_frac128(idx_p, u); // lsb = 2^-127
+
+ // Scale m_final and 2^-hi_p to lsb = 2^-126 so that both
+ // num = m_final - 2^-hi_p in (0.28, 2.0)
+ // den = m_final + 2^-hi_p in (1.0, 3.0)
+ // fit in Frac128 without overflow:
+ Frac128 m_half = m_final >> 1; // lsb = 2^-126
+ Frac128 one_neg_hi({0, 1ULL << (62 - hi_p)}); // 2^-hi_p, lsb = 2^-126
+ Frac128 num = m_half - one_neg_hi; // lsb = 2^-126
+ Frac128 den = m_half + one_neg_hi; // lsb = 2^-126
+
+ // Normalize den to d_norm in (1, 2) with lsb = 2^-127 (shift_d is 0 if
+ // m_final + 2^-hi_p >= 2, or 1 if m_final + 2^-hi_p < 2):
+ unsigned shift_d = static_cast<unsigned>((den.val[1] >> 63) == 0);
+ Frac128 d_norm = den << shift_d; // in (1, 2), lsb = 2^-127
+ Frac128 n_norm = num << shift_d; // in (0.12, d_norm), lsb = 2^-127
+
+ // Compute 126-bit reciprocal r2 ~ 1 / d_norm (lsb = 2^-128) via two
+ // Newton-Raphson iterations starting from a 53-bit double reciprocal:
+ // r_{k+1} = r_k * (2 - d_norm * r_k).
+ // Since d_norm * r_k has lsb = 2^-127, ~dr_k represents (2 - d_norm * r_k)
+ // with lsb = 2^-127 to within 2^-127.
+ // Convert (1.0 / d_d) * 2^62 in [2^61, 2^62] via signed int64_t and shift
+ // left by 2 to avoid spurious FE_INVALID on x86-64 double-to-uint64_t
+ // conversions of values >= 2^63:
+ int64_t d_hi_53 = static_cast<int64_t>((d_norm.val[1] >> 11) + 1);
+ double d_d = static_cast<double>(d_hi_53) * 0x1.0p-52;
+ double r0_d = (1.0 / d_d) * 0x1.0p62;
+ uint64_t r0_hi = static_cast<uint64_t>(static_cast<int64_t>(r0_d)) << 2;
+ Frac128 r0({0, r0_hi}); // lsb = 2^-128 (~52-bit accurate)
+
+ Frac128 dr0 = d_norm * r0; // lsb = 2^-127
+ Frac128 r1 = (r0 * (~dr0)) << 1; // lsb = 2^-128 (~104-bit accurate)
+
+ Frac128 dr1 = d_norm * r1; // lsb = 2^-127
+ Frac128 r2 = (r1 * (~dr1)) << 1; // lsb = 2^-128 (~126-bit accurate)
+
+ // Quotient q = n_norm * r2 ~ num / den = tanh(|x|) in (0.12, 1) with
+ // lsb = 2^-127:
+ Frac128 q = n_norm * r2; // lsb = 2^-127
+ DFloat128 r(xbits.sign(), -127, UInt<128>(q));
+ return static_cast<double>(r);
+}
+
+} // namespace double_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_TANH_ACCURATE_128_H
diff --git a/libc/src/__support/math/tanh_double_eval.h b/libc/src/__support/math/tanh_double_eval.h
new file mode 100644
index 00000000000000..358ec633c08953
--- /dev/null
+++ b/libc/src/__support/math/tanh_double_eval.h
@@ -0,0 +1,559 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 tanh.
+///
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_TANH_DOUBLE_EVAL_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_TANH_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/double_double.h"
+#include "src/__support/FPUtil/except_value_utils.h"
+#include "src/__support/FPUtil/fp_bitwise.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/attributes.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+#include "src/__support/macros/properties/cpu_features.h"
+#include "src/__support/math/exp_utils.h"
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+#include "src/__support/math/tanh_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) / 2 ~ hi + lo with 35-bit-precision hi part so that
+// kd * MLOG_2_HALF_HI is exact for |kd| < 2^18:
+// > hi = round(-log(2) / 2, 35, RN);
+// > lo = round(-log(2) / 2 - hi, D, RN);
+LIBC_INLINE_VAR constexpr double MLOG_2_HALF_HI = -0x1.62e42fefcp-2;
+LIBC_INLINE_VAR constexpr double MLOG_2_HALF_LO = 0x1.c610ca86c3899p-38;
+
+// Degree-3 minimax polynomials in v^2 for (cosh(v) - 1) / v^2 and
+// (sinh(v) - v) / v^3 on [-log(2) / 4, log(2) / 4]:
+// cosh(v) - 1 ~ v^2 * P_even(v^2)
+// sinh(v) - v ~ v^3 * P_odd(v^2).
+//
+// Generated by Sollya with:
+// > display = hexadecimal!;
+// > prec = 500;
+// > B = (log(2) / 4)^2;
+// > P_even = fpminimax((cosh(sqrt(x)) - 1) / x, 3, [|D...|], [2^(-100), B]);
+// > for i from 0 to 3 do print(coeff(P_even, i), ",");
+// > dirtyinfnorm((cosh(sqrt(x)) - 1) - x * P_even, [0, B]);
+// 0x1.e49e...p-55 < 2^-54.
+LIBC_INLINE_VAR constexpr double TANH_EVEN_N0_COEFFS[4] = {
+ 0x1.fffffffffffep-2, // ~ 1/2!
+ 0x1.5555555597352p-5, // ~ 1/4!
+ 0x1.6c16bc13d272cp-10, // ~ 1/6!
+ 0x1.a06123b7cd784p-16, // ~ 1/8!
+};
+
+// > P_odd = fpminimax((sinh(sqrt(x)) / sqrt(x) - 1) / x, 3, [|D...|],
+// [2^(-100), B]);
+// > for i from 0 to 3 do print(coeff(P_odd, i), ",");
+// > dirtyinfnorm((sinh(sqrt(x)) - sqrt(x)) - x * sqrt(x) * P_odd, [0, B]);
+// 0x1.ee40...p-61 < 2^-60.
+LIBC_INLINE_VAR constexpr double TANH_ODD_N0_COEFFS[4] = {
+ 0x1.555555555555p-3, // ~ 1/3!
+ 0x1.111111112884cp-7, // ~ 1/5!
+ 0x1.a019fdc55d34fp-13, // ~ 1/7!
+ 0x1.7211ce0221666p-19, // ~ 1/9!
+};
+
+// Tableless (N = 0) evaluation of tanh(x) = sinh(x) / cosh(x) for
+// 2^-27 < |x| < 0x1.30fc1931f09cap+4.
+//
+// Range reduction:
+// Let |x| = k * (log(2) / 2) + v, where:
+// k = round(|x| * 2 * log2(e))
+// |v| <= log(2) / 4
+// then:
+// e^|x| = 2^(k / 2) * e^v
+// e^-|x| = 2^(-k / 2) * e^-v = 2^(k / 2) * 2^-k * e^-v.
+//
+// Factoring out 2^(k / 2) in numerator and denominator and separating
+// e^(+/- v) = cosh(v) +/- sinh(v):
+// tanh(|x|) = (e^v - m_m * e^-v) / (e^v + m_m * e^-v)
+// = (d * cosh(v) + s * sinh(v)) / (s * cosh(v) + d * sinh(v))
+// where:
+// m_m = 2^-k
+// s = 1 + m_m
+// d = 1 - m_m.
+//
+// Key properties across the domain:
+// - For |x| < log(2) / 4, k = 0 gives v = |x| (exact, t2 = 0) and m_m = 1, so
+// d = 0 and s = 2, reducing the ratio to:
+// (2 * sinh(v)) / (2 * cosh(v)) = sinh(x) / cosh(x)
+// with no cancellation down to |x| = 2^-27.
+//
+// - For 1 <= k <= 53, 2^-k is in [2^-53, 2^-1], so both:
+// s = 1.0 + m_m
+// d = 1.0 - m_m
+// are exact in all rounding modes. For k >= 54:
+// 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 tanh_eval(double x) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ double x_a = xbits.abs().get_val();
+
+ // Using ROUND_BIAS, ROUND_SHIFT, and ROUND_MASK (see exp_utils.h), we compute
+ // k_biased = fma(x_a, LOG2_E_2, ROUND_BIAS) and extract in parallel:
+ // kd = (k_biased & ROUND_MASK) - ROUND_SHIFT
+ // k = (k_bits >> 32) & 0x3'ffff.
+ constexpr double LOG2_E_2 = 2.0 * LOG2_E.hi;
+ double k_biased = fputil::multiply_add(x_a, LOG2_E_2, 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 = 2^-k for k < 54
+ // m_m = 0.0 for k >= 54 (where 2^-k <= 2^-54).
+ double m_m = (k < 54) ? scale_exp2(1.0, -k) : 0.0;
+
+ // Reduced argument:
+ // v = (x_a + kd * MLOG_2_HALF_HI) + kd * MLOG_2_HALF_LO.
+ // Since MLOG_2_HALF_HI has 35 significant bits and |kd| <= 55,
+ // t1 = fma(kd, MLOG_2_HALF_HI, x_a)
+ // is exact in all rounding modes with or without FMA.
+ double t1 = fputil::multiply_add(kd, MLOG_2_HALF_HI, x_a); // exact
+ double t2 = kd * MLOG_2_HALF_LO;
+ double v = t1 + t2;
+ double v2 = v * v;
+ double v4 = v2 * v2;
+
+ // Evaluate degree-3 even and odd minimax polynomials P_even(v^2) and
+ // P_odd(v^2) in parallel using Estrin's scheme:
+ // cosh(v) ~ 1 + v^2 * P_even(v^2)
+ // sinh(v) ~ t1 + (t2 + v^3 * P_odd(v^2))
+
+ // First Estrin layer: c_0 + v^2 * c_1 and c_2 + v^2 * c_3
+ double e01 =
+ fputil::multiply_add(v2, TANH_EVEN_N0_COEFFS[1], TANH_EVEN_N0_COEFFS[0]);
+ double o01 =
+ fputil::multiply_add(v2, TANH_ODD_N0_COEFFS[1], TANH_ODD_N0_COEFFS[0]);
+ double e23 =
+ fputil::multiply_add(v2, TANH_EVEN_N0_COEFFS[3], TANH_EVEN_N0_COEFFS[2]);
+ double o23 =
+ fputil::multiply_add(v2, TANH_ODD_N0_COEFFS[3], TANH_ODD_N0_COEFFS[2]);
+
+ // Final polynomial values P_even(v^2) and P_odd(v^2):
+ double p_even = fputil::multiply_add(v4, e23, e01);
+ double p_odd = fputil::multiply_add(v4, o23, o01);
+
+ // Low parts of cosh(v) - 1 and sinh(v) - t1:
+ // c_lo ~ v^2 * P_even(v^2)
+ // s_lo ~ t2 + v^3 * P_odd(v^2)
+ double c_lo = v2 * p_even;
+ double v3 = v * v2;
+ double s_lo = fputil::multiply_add(v3, p_odd, t2);
+
+ // Since m_m is either 0 or 2^-k with 0 <= k <= 53, 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 numerator weights s_num and d_num before accumulating
+ // the polynomial so the numerator rounds in the active rounding direction for
+ // negative x as well:
+ double sign_x = fputil::fp_xor(x, x_a);
+ double s_num = fputil::fp_xor(s, sign_x);
+ double d_num = fputil::fp_xor(d, sign_x);
+
+ // Numerator:
+ // num ~ sign(x) * (d * (1 + c_lo) + s * (t1 + s_lo))
+ // = d_num + (s_num * t1 + (d_num * c_lo + s_num * s_lo))
+ double s_s_lo = s_num * s_lo;
+ double num_lo = fputil::multiply_add(d_num, c_lo, s_s_lo);
+ double num_mid = fputil::multiply_add(s_num, t1, num_lo);
+ double num = d_num + num_mid;
+
+ // Denominator:
+ // den ~ s * (1 + c_lo) + d * (t1 + s_lo)
+ // = s + (d * t1 + (s * c_lo + d * s_lo))
+ double d_s_lo = d * s_lo;
+ double den_lo = fputil::multiply_add(s, c_lo, d_s_lo);
+ double den_mid = fputil::multiply_add(d, t1, den_lo);
+ double den = s + den_mid;
+
+ return num / den;
+}
+
+#else // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// Leading 11-bit split of c_0 ~ -1/3 for exact cubic reconstruction:
+// TANH_C0_HI_11 = round(-1/3, 11, RZ) = -0x1.554p-2
+// -1/3 - TANH_C0_HI_11 = -0x1.5555555555555p-14.
+LIBC_INLINE_VAR constexpr double TANH_C0_HI_11 = -0x1.554p-2;
+
+// Polynomial approximation for ((tanh(x) - x) / x^3 - (-0x1.554p-2)) on
+// [-2^-3, 2^-3] as a degree-5 even polynomial in x^2:
+// tanh(x) - x ~ x^3 * (TANH_C0_HI_11 + P_rem(x^2))
+//
+// Minimax polynomial coefficients generated by Sollya with:
+// > display = hexadecimal!;
+// > prec = 500;
+// > c0 = -0x1.554p-2 + (-0x1.5555555555555p-14);
+// > P = fpminimax((tanh(sqrt(x))/sqrt(x) - 1)/x - c0, [|1, 2, 3, 4, 5|],
+// [|D...|], [2^(-100), 2^(-6)], absolute);
+// > for i from 1 to 5 do print(coeff(P, i), ",");
+// > dirtyinfnorm((tanh(x) - x)/x^3 - (c0 + P(x^2)), [2^(-50), 2^(-3)]);
+// 0x1.9e2e...p-57 < 2^-56.
+// > dirtyinfnorm(tanh(x) - (x + x^3 * (c0 + P(x^2))), [-2^(-3), 2^(-3)]);
+// 0x1.9dd6...p-66 < 2^-65.
+LIBC_INLINE_VAR constexpr double TANH_SMALL_P_COEFFS[6] = {
+ -0x1.5555555555555p-14, // c_0 - (-0x1.554p-2) ~ -1/3 + 0x1.554p-2
+ 0x1.1111111110c74p-3, // c_1 ~ 2/15
+ -0x1.ba1ba1b6047f3p-5, // c_2 ~ -17/315
+ 0x1.664f3f15b86a5p-6, // c_3 ~ 62/2835
+ -0x1.2264b8b90e812p-7, // c_4 ~ -1382/155925
+ 0x1.ce2449b1b72dcp-9, // c_5 ~ 21844/6081075
+};
+
+// Evaluate P_rem(x^2) ~ (tanh(x) - x) / x^3 - (-0x1.554p-2) using Estrin's
+// scheme:
+LIBC_ALWAYS_INLINE double tanh_small_poly(double x2) {
+ double x4 = x2 * x2;
+ double x8 = x4 * x4;
+ double c0 =
+ fputil::multiply_add(x2, TANH_SMALL_P_COEFFS[1], TANH_SMALL_P_COEFFS[0]);
+ double c1 =
+ fputil::multiply_add(x2, TANH_SMALL_P_COEFFS[3], TANH_SMALL_P_COEFFS[2]);
+ double c2 =
+ fputil::multiply_add(x2, TANH_SMALL_P_COEFFS[5], TANH_SMALL_P_COEFFS[4]);
+ double c01 = fputil::multiply_add(x4, c1, c0);
+ return fputil::multiply_add(x8, c2, c01);
+}
+
+// Evaluate tanh(x) in double-double precision for
+// 2^-3 <= |x| < 0x1.30fc1931f09cap+4, falling back to tanh_accurate_128(x) when
+// Ziv's rounding test fails.
+//
+// Range reduction:
+// Let 2 * |x| = k_p * (log(2) / 2^7) + u
+// in which:
+// k_p = round(2 * |x| * 2^7 * log2(e)), k_m = -k_p
+// hi_m = k_m >> 7, idx_m = k_m & 0x7f
+// |u| <= 2^-8 * log(2)
+// then:
+// e^(2 * |x|) = 2^(k_p / 2^7) * e^u = e^u / m
+// where:
+// m = 2^hi_m * 2^(idx_m / 2^7) = 2^hi_m * EXP2_MID[idx_m].
+//
+// Rewriting tanh(|x|) = (e^(2 * |x|) - 1) / (e^(2 * |x|) + 1) in terms of e^u
+// and m avoids multiplying m by e^u:
+// tanh(|x|) = (e^u - m) / (e^u + m) = 1 - (2 * m) / (e^u + m).
+LIBC_MATH_SMALL_TABLES_INLINE double tanh_eval(double x) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ double x_a = xbits.abs().get_val();
+ double x_2 = 2.0 * x_a;
+
+ 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_2, 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_2, 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_m = k_m >> EXP_N_BITS;
+ unsigned idx_m = static_cast<unsigned>(k_m & EXP_MASK);
+
+ // Scale m = EXP2_MID[idx_m] by 2^hi_m:
+ // For 2^-3 <= |x| < 0x1.30fc1931f09cap+4, hi_m is in [-56, 0], so scaling by
+ // 2^hi_m is always in the normal range:
+ // - Since m.hi is in [1, 2), we use scale_exp2(m.hi, hi_m) directly.
+ // - Since EXP2_MID[0].lo == 0.0, we scale m.lo by multiplying by:
+ // s_m = scale_exp2(1.0, hi_m).
+ double s_m = scale_exp2(1.0, hi_m);
+ DoubleDouble m = EXP2_MID[idx_m];
+ m.lo *= s_m;
+ m.hi = scale_exp2(m.hi, hi_m);
+
+ // Reduced argument:
+ // u = (x_2 + kd * MLOG_2_EXP2_HI) + kd * MLOG_2_EXP2_MID.
+ // Since MLOG_2_EXP2_HI has 35 significant bits and |kd| < 2^13,
+ // t1 = fma(kd, MLOG_2_EXP2_HI, x_2)
+ // is exact in all rounding modes with or without FMA.
+ double t1 = fputil::multiply_add(kd, MLOG_2_EXP2_HI, x_2); // exact
+ double t2 = kd * MLOG_2_EXP2_MID;
+ double u = t1 + t2;
+
+ // Evaluate degree-3 minimax polynomial P(u) ~ (e^u - 1 - u) / u^2:
+ // e^u ~ 1 + t1 + exp_lo
+ // where:
+ // exp_lo = t2 + u^2 * P(u).
+ double u2 = u * u;
+ double p = exp_poly(u, u2);
+ double exp_lo = fputil::multiply_add(u2, p, t2);
+
+ // Denominator D = e^u + m ~ (1.0 + m.hi + t1 + exp_lo) + m.lo:
+ // s.hi + s.lo = 1.0 + m.hi (exact FastTwoSum since m.hi <= 1.0)
+ // d1.hi + d1.lo = s.hi + t1 (FastTwoSum, |t1| < 2^-8 < s.hi)
+ // d.hi + d.lo = d1.hi + exp_lo (FastTwoSum, |exp_lo| < 2^-17 < d1.hi).
+ // Accumulating exp_lo into d.hi keeps |d_lo| <= 2^-50 so the linear
+ // Karp-Markstein remainder is accurate to 105 bits:
+ DoubleDouble s = fputil::exact_add(1.0, m.hi);
+ DoubleDouble d1 = fputil::exact_add(s.hi, t1);
+ DoubleDouble d = fputil::exact_add(d1.hi, exp_lo);
+
+ // Combine FastTwoSum residuals and table low part:
+ // d_lo = s.lo + d1.lo + d.lo + m.lo
+ double sd1_lo = s.lo + d1.lo;
+ double sdd_lo = sd1_lo + d.lo;
+ double d_lo = sdd_lo + m.lo;
+
+ // Karp-Markstein division for q_hi + q_lo ~ (2 * m) / (d.hi + d_lo):
+ // y_hi = round(1.0 / d.hi)
+ // q_hi = round(2 * m.hi * y_hi)
+ // rem_hi = 2 * m.hi - q_hi * d.hi
+ // rem_lo = 2 * m.lo - q_hi * d_lo
+ // q_lo = (rem_hi + rem_lo) * y_hi.
+ double two_m_hi = 2.0 * m.hi;
+ double two_m_lo = 2.0 * m.lo;
+ double y_hi = 1.0 / d.hi;
+ double q_hi = two_m_hi * y_hi;
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double rem_hi = fputil::multiply_add(-q_hi, d.hi, two_m_hi);
+ double rem_lo = fputil::multiply_add(-q_hi, d_lo, two_m_lo);
+#else
+ // Split q_hi into 26 + 27 bits and d.hi into 27 + 26 bits via bit-masks so
+ // that q_hi_hi * d_hi_hi (53 bits) and q_hi_hi * d_hi_lo (52 bits) are exact
+ // without FMA:
+ LIBC_BIT_CAST_CONSTEXPR_VAR double HI_MASK_26 =
+ cpp::bit_cast<double>(0xffff'ffff'f800'0000ULL);
+ LIBC_BIT_CAST_CONSTEXPR_VAR double HI_MASK_27 =
+ cpp::bit_cast<double>(0xffff'ffff'fc00'0000ULL);
+ double q_hi_hi = fputil::fp_and(q_hi, HI_MASK_26);
+ double q_hi_lo = q_hi - q_hi_hi;
+ double d_hi_hi = fputil::fp_and(d.hi, HI_MASK_27);
+ double d_hi_lo = d.hi - d_hi_hi;
+
+ // High remainder rem_hi = 2 * m.hi - q_hi * d.hi:
+ double prod_hh = q_hi_hi * d_hi_hi; // exact (53 bits)
+ double prod_hl = q_hi_hi * d_hi_lo; // exact (52 bits)
+ double prod_lh = q_hi_lo * d_hi_hi;
+ double prod_ll = q_hi_lo * d_hi_lo;
+
+ double rem_1 = two_m_hi - prod_hh; // exact
+ double rem_2 = rem_1 - prod_hl;
+ double rem_3 = rem_2 - prod_lh;
+ double rem_hi = rem_3 - prod_ll;
+
+ // Low remainder rem_lo = 2 * m.lo - q_hi * d_lo:
+ double q_d_lo = q_hi * d_lo;
+ double rem_lo = two_m_lo - q_d_lo;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double rem = rem_hi + rem_lo;
+ double q_lo = rem * y_hi;
+
+ // Subtract q_hi + q_lo from 1.0:
+ // r.hi + r.lo = 1.0 - q_hi (FastTwoSum)
+ // r_lo = r.lo - q_lo.
+ // Since 1.0 is an exact power of two and q_hi is in (0, 0.88),
+ // exact_sub(1.0, q_hi) is accurate to within 2^-105 in all rounding modes.
+ DoubleDouble r = fputil::exact_sub(1.0, q_hi);
+ double r_lo = r.lo - q_lo;
+
+ // Apply sign(x) before Ziv's rounding test 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);
+ double r_hi = fputil::fp_xor(r.hi, sign_x);
+ r_lo = fputil::fp_xor(r_lo, sign_x);
+
+ // Error bound for Ziv's rounding test:
+ // Since any error in e^u or q is damped by q_hi ~ 2 * e^(-2 * |x|), we scale
+ // the relative error bound by q_hi:
+ // err = q_hi * TANH_ERR_REL + TANH_ERR_ABS.
+ constexpr double TANH_ERR_REL = 0x1.8p-63;
+ constexpr double TANH_ERR_ABS = 0x1.0p-104;
+ double err = fputil::multiply_add(q_hi, TANH_ERR_REL, TANH_ERR_ABS);
+ double lo_u = r_lo + err;
+ double lo_l = r_lo - err;
+
+ // 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 tanh_accurate_128(x);
+}
+
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
+// Evaluate tanh(x) for |x| < SMALL_BOUND (2^-27 + ulp(2^-27) when
+// LIBC_MATH_HAS_SKIP_ACCURATE_PASS is defined, or 2^-3 otherwise):
+LIBC_INLINE double tanh_small_eval(double x, uint64_t x_abs) {
+ using FPBits = fputil::FPBits<double>;
+
+ // Tiny and subnormal range: |x| <= 2^-27.
+ // For 0 < |x| <= 2^-27:
+ // 0 < (|x| - |tanh(x)|) / |x| ~ x^2 / 3 < 2^-55 < 2^-1 * ulp(x),
+ // so tanh(x) rounds to:
+ // - nextafter(x, 0) when rounding mode is FE_TOWARDZERO, or (x > 0 and
+ // FE_DOWNWARD), or (x < 0 and FE_UPWARD),
+ // - x otherwise, with FE_INEXACT raised for all x != 0 and FE_UNDERFLOW +
+ // errno = ERANGE raised when |tanh(x)| < 2^-1022.
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ if (LIBC_UNLIKELY(x_abs <= 0x3e40'0000'0000'0000ULL))
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ {
+ if (LIBC_UNLIKELY(x_abs == 0))
+ return x;
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ double res = fputil::multiply_add(x, -0x1.0p-55, x);
+#else
+ // Avoid spurious FE_UNDERFLOW on x * -2^-55 when 2^-1022 < |x| < 2^-967,
+ // and preserve directed rounding for subnormal |x| <= 2^-1022:
+ double x_scaled = fputil::multiply_add(x, 0x1.0p55, -x);
+ double res = x_scaled * 0x1.0p-55;
+#endif // LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ if (LIBC_UNLIKELY(x_abs <= 0x0010'0000'0000'0000ULL)) {
+ if (FPBits(res).abs().uintval() < 0x0010'0000'0000'0000ULL) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_underflow_except_if_required<double>();
+ }
+ }
+ return res;
+ }
+
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // 2^-27 < |x| < 2^-3:
+ // 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 * TANH_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 * TANH_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 * (TANH_C0_HI_11 + P_rem(x^2)) = y_hi + y_lo
+ // where:
+ // y_lo ~ x_hi^3 * P_rem(x^2) + x3_lo * (TANH_C0_HI_11 + P_rem(x^2))
+ double p_rem = tanh_small_poly(x2);
+ double p_full = TANH_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^-7 * |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| * 0x1.8p-54, 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.8p-54, 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 tanh_small_accurate_128(x);
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+}
+
+// Handle saturation, inf, and nan for tanh(x) when
+// |x| >= 0x1.30fc1931f09cap+4:
+LIBC_NOINLINE LIBC_INLINE double tanh_exceptional(double x) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+
+ if (xbits.is_nan())
+ return x + 1.0;
+
+ if (xbits.is_inf())
+ return xbits.is_neg() ? -1.0 : 1.0;
+
+ return xbits.is_neg() ? fputil::round_result_slightly_up(-1.0)
+ : fputil::round_result_slightly_down(1.0);
+}
+
+LIBC_ALWAYS_INLINE double tanh(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^-27 + ulp(2^-27) (tanh_eval handles |x| > 2^-27 directly
+ // with k = 0 and zero cancellation when |x| < log(2) / 4)
+ constexpr uint64_t SMALL_BOUND = 0x3e40'0000'0000'0001ULL;
+#else
+ // - Small bound: 2^-3 = 0x1.0p-3
+ constexpr uint64_t SMALL_BOUND = 0x3fc0'0000'0000'0000ULL;
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // - Saturation bound:
+ // > round(atanh(1 - 2^-54), D, RU) = 0x1.30fc1931f09cap+4
+ // > round(atanh(1 - 2^-54), D, RD) = 0x1.30fc1931f09c9p+4
+ // > round(tanh(0x1.30fc1931f09c9p+4), D, RN) = 0x1.fffffffffffffp-1
+ // > round(tanh(0x1.30fc1931f09cap+4), D, RN) = 0x1.0p+0
+ constexpr uint64_t OVERFLOW_BOUND = 0x4033'0fc1'931f'09caULL;
+
+ if (LIBC_UNLIKELY(x_abs < SMALL_BOUND || x_abs >= OVERFLOW_BOUND)) {
+ if (x_abs < SMALL_BOUND)
+ return tanh_small_eval(x, x_abs);
+
+ return tanh_exceptional(x);
+ }
+
+ return tanh_eval(x);
+}
+
+} // namespace double_eval
+} // namespace math
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_TANH_DOUBLE_EVAL_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 9f0325832248ec..f1242e15ed7243 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -3793,6 +3793,16 @@ add_entrypoint_object(
libc.src.__support.math.sinhf16
)
+add_entrypoint_object(
+ tanh
+ SRCS
+ tanh.cpp
+ HDRS
+ ../tanh.h
+ DEPENDS
+ libc.src.__support.math.tanh
+)
+
add_entrypoint_object(
tanhf
SRCS
diff --git a/libc/src/math/generic/tanh.cpp b/libc/src/math/generic/tanh.cpp
new file mode 100644
index 00000000000000..721da730bf059e
--- /dev/null
+++ b/libc/src/math/generic/tanh.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 tanh(x) function.
+///
+//===----------------------------------------------------------------------===//
+
+#include "src/math/tanh.h"
+#include "src/__support/math/tanh.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(double, tanh, (double x)) { return math::tanh(x); }
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index f902fa553e49ef..ca47f4c3f834a1 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -488,6 +488,7 @@ else()
libc.src.__support.math.tanbf16
libc.src.__support.math.tanf
libc.src.__support.math.tanf16
+ libc.src.__support.math.tanh
libc.src.__support.math.tanhf
libc.src.__support.math.tanhf16
libc.src.__support.math.tanpif
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 73a5cb45ed2b9c..5cce92b3bda827 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -357,6 +357,7 @@ TEST(LlvmLibcSharedMathTest, AllDouble) {
EXPECT_FP_EQ(0.0, sin);
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::sqrt(0.0));
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::tan(0.0));
+ EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::tanh(0.0));
EXPECT_EQ(0, LIBC_NAMESPACE::shared::ilogb(1.0));
EXPECT_EQ(0L, LIBC_NAMESPACE::shared::llogb(1.0));
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::logb(1.0));
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 9ec8422b57f081..e20164868eda3f 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2654,6 +2654,22 @@ add_fp_unittest(
libc.src.math.sinhf16
)
+add_fp_unittest(
+ tanh_test
+ NEED_MPFR
+ SUITE
+ libc-math-tests
+ SRCS
+ tanh_test.cpp
+ DEPENDS
+ libc.hdr.fenv_macros
+ libc.src.math.tanh
+ libc.src.__support.math.tanh_double_eval
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.optimization
+)
+
add_fp_unittest(
tanhf_test
NEED_MPFR
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 7485f5409402f6..3c094433ab8bda 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -4814,6 +4814,22 @@ add_fp_unittest(
libc.src.__support.FPUtil.cast
)
+add_fp_unittest(
+ tanh_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ tanh_test.cpp
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.math.tanh
+ libc.src.__support.math.tanh_double_eval
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.macros.attributes
+ libc.src.__support.macros.optimization
+)
+
add_fp_unittest(
tanhf_test
SUITE
diff --git a/libc/test/src/math/smoke/tanh_test.cpp b/libc/test/src/math/smoke/tanh_test.cpp
new file mode 100644
index 00000000000000..2a4fcbf032557e
--- /dev/null
+++ b/libc/test/src/math/smoke/tanh_test.cpp
@@ -0,0 +1,146 @@
+//===----------------------------------------------------------------------===//
+//
+// 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 tanh.
+///
+//===----------------------------------------------------------------------===//
+
+#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/macros/optimization.h"
+#include "src/__support/math/tanh_double_eval.h"
+#include "src/math/tanh.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+LIBC_NOINLINE double tanh_double_eval(double x) {
+ return LIBC_NAMESPACE::math::double_eval::tanh(x);
+}
+
+class TanhTest : 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(1.0, func(inf));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(-1.0, 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_saturation(double (*func)(double)) {
+ // x = 0x1.fffffffffffffp+1023 (max_normal)
+ EXPECT_FP_EQ_ALL_ROUNDING(
+ 1.0, 1.0, 0x1.fffffffffffffp-1, 0x1.fffffffffffffp-1,
+ func(FPBits(uint64_t(0x7fef'ffff'ffff'ffff)).get_val()));
+ EXPECT_MATH_ERRNO(0);
+
+ // x = 0x1.30fc1931f09cap+4 (smallest positive input that saturates to 1.0
+ // in FE_TONEAREST)
+ EXPECT_FP_EQ_ALL_ROUNDING(
+ 1.0, 1.0, 0x1.fffffffffffffp-1, 0x1.fffffffffffffp-1,
+ func(FPBits(uint64_t(0x4033'0fc1'931f'09ca)).get_val()));
+ EXPECT_MATH_ERRNO(0);
+
+ // x = -0x1.30fc1931f09cap+4 (largest negative input that saturates to -1.0
+ // in FE_TONEAREST)
+ EXPECT_FP_EQ_ALL_ROUNDING(
+ -1.0, -0x1.fffffffffffffp-1, -1.0, -0x1.fffffffffffffp-1,
+ func(FPBits(uint64_t(0xc033'0fc1'931f'09ca)).get_val()));
+ EXPECT_MATH_ERRNO(0);
+ }
+
+ 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);
+ }
+
+ void test_accurate_pass(double (*func)(double)) {
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ // Hard-to-round inputs exercising tanh_small_accurate_128 and
+ // tanh_accurate_128 without raising spurious exceptions:
+ EXPECT_FP_EQ(0x1.fde1b6fff6ef9p-4, func(0x1.004524207fad9p-3));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(-0x1.fde1b6fff6ef9p-4, func(-0x1.004524207fad9p-3));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(0x1.fefff82a29e4dp-5, func(0x1.ffaa08fcad74ep-5));
+ EXPECT_MATH_ERRNO(0);
+#else
+ (void)func;
+#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ }
+
+#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_TANH_FTZ_DAZ_TESTS(suffix, func) \
+ TEST_F(LlvmLibcTanhTest##suffix, FTZMode) { \
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::FTZ); \
+ test_denormals(&func); \
+ } \
+ TEST_F(LlvmLibcTanhTest##suffix, DAZMode) { \
+ LIBC_NAMESPACE::testing::ModifyMXCSR mxcsr(LIBC_NAMESPACE::testing::DAZ); \
+ test_denormals(&func); \
+ } \
+ TEST_F(LlvmLibcTanhTest##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_TANH_FTZ_DAZ_TESTS(suffix, func)
+#endif // LIBC_TEST_FTZ_DAZ
+
+#define LIST_TANH_TESTS(suffix, func, check_snan_invalid) \
+ using LlvmLibcTanhTest##suffix = TanhTest; \
+ TEST_F(LlvmLibcTanhTest##suffix, SpecialNumbers) { \
+ test_special_numbers(&func, check_snan_invalid); \
+ } \
+ TEST_F(LlvmLibcTanhTest##suffix, Saturation) { test_saturation(&func); } \
+ TEST_F(LlvmLibcTanhTest##suffix, Underflow) { test_underflow(&func); } \
+ TEST_F(LlvmLibcTanhTest##suffix, AccuratePass) { \
+ test_accurate_pass(&func); \
+ } \
+ LIST_TANH_FTZ_DAZ_TESTS(suffix, func) \
+ static_assert(true, "Require semicolon.")
+
+LIST_TANH_TESTS(Default, LIBC_NAMESPACE::tanh, /*check_snan_invalid=*/true);
+LIST_TANH_TESTS(DoubleEval, tanh_double_eval, /*check_snan_invalid=*/false);
diff --git a/libc/test/src/math/tanh_test.cpp b/libc/test/src/math/tanh_test.cpp
new file mode 100644
index 00000000000000..42e93b2ec3791f
--- /dev/null
+++ b/libc/test/src/math/tanh_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 tanh.
+///
+//===----------------------------------------------------------------------===//
+
+#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/tanh_double_eval.h"
+#include "src/math/tanh.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 tanh_double_eval(double x) {
+ return LIBC_NAMESPACE::math::double_eval::tanh(x);
+}
+
+class TanhTest : public LIBC_NAMESPACE::testing::FPTest<double> {
+public:
+ void test_special_numbers(double (*func)(double)) {
+ EXPECT_FP_EQ(aNaN, func(aNaN));
+ EXPECT_FP_EQ_ALL_ROUNDING(1.0, func(inf));
+ EXPECT_FP_EQ_ALL_ROUNDING(-1.0, 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_tricky_inputs(double (*func)(double), double ulp_tolerance = 0.5,
+ bool all_rounding = true) {
+ constexpr int N = 22;
+ constexpr uint64_t INPUTS[N] = {
+ 0x3e40'0000'0000'0000ULL, // x = 0x1.0p-27
+ 0x3e40'0000'0000'0001ULL, // x = 0x1.0000000000001p-27
+ 0x3e4d'12ed'0af1'a27fULL, // x = 0x1.d12ed0af1a27fp-27
+ 0x3fbf'faa0'8fca'd74eULL, // x = 0x1.ffaa08fcad74ep-5
+ 0x3fbf'ffff'ffff'ffffULL, // x = 0x1.fffffffffffffp-4
+ 0x3fc0'0000'0000'0000ULL, // x = 0x1.0p-3
+ 0x3fc0'0452'4207'fad9ULL, // x = 0x1.004524207fad9p-3
+ 0x3fca'c343'b179'fec4ULL, // x = 0x1.ac343b179fec4p-3
+ 0x3fd0'0764'a988'bf73ULL, // x = 0x1.00764a988bf73p-2
+ 0x3fd1'7d1e'8a63'711fULL, // x = 0x1.17d1e8a63711fp-2
+ 0x3fd2'91c6'01a0'5276ULL, // x = 0x1.291c601a05276p-2
+ 0x3fd3'6f33'd51c'264dULL, // x = 0x1.36f33d51c264dp-2
+ 0x3fd4'3eae'a236'49c3ULL, // x = 0x1.43eaea23649c3p-2
+ 0x3fdd'88d7'550b'2826ULL, // x = 0x1.d88d7550b2826p-2
+ 0x3fde'611a'a58a'b608ULL, // x = 0x1.e611aa58ab608p-2
+ 0x3fe0'1efe'7ac8'c15dULL, // x = 0x1.01efe7ac8c15dp-1
+ 0x3fe1'005e'c0bc'cabbULL, // x = 0x1.1005ec0bccabbp-1
+ 0x3fe3'3dfe'b0fa'4bfeULL, // x = 0x1.33dfeb0fa4bfep-1
+ 0x3fe4'9f24'ac5c'ac35ULL, // x = 0x1.49f24ac5cac35p-1
+ 0x4010'9cc2'de69'e78cULL, // x = 0x1.09cc2de69e78cp+2
+ 0x4033'0fc1'931f'09c9ULL, // x = 0x1.30fc1931f09c9p+4
+ 0x4033'0fc1'931f'09caULL, // x = 0x1.30fc1931f09cap+4
+ };
+ for (int i = 0; i < N; ++i) {
+ double x = FPBits(INPUTS[i]).get_val();
+ if (all_rounding) {
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tanh, x, func(x),
+ ulp_tolerance);
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tanh, -x, func(-x),
+ ulp_tolerance);
+ } else {
+ EXPECT_MPFR_MATCH(mpfr::Operation::Tanh, x, func(x), ulp_tolerance);
+ EXPECT_MPFR_MATCH(mpfr::Operation::Tanh, -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(20.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::Tanh, x, result,
+ ulp_tolerance, rounding_mode)) {
+ ++fails;
+ while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(
+ mpfr::Operation::Tanh, x, result, tol, rounding_mode)) {
+ mx = x;
+ mr = result;
+
+ if (tol > 1000.0)
+ break;
+
+ tol *= 2.0;
+ }
+ }
+ }
+ if (fails) {
+ tlog << " Tanh failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol)
+ << ".\n";
+ EXPECT_MPFR_MATCH(mpfr::Operation::Tanh, 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_TANH_TESTS(suffix, func, ulp_tolerance, all_rounding) \
+ using LlvmLibcTanhTest##suffix = TanhTest; \
+ TEST_F(LlvmLibcTanhTest##suffix, SpecialNumbers) { \
+ test_special_numbers(&func); \
+ } \
+ TEST_F(LlvmLibcTanhTest##suffix, TrickyInputs) { \
+ test_tricky_inputs(&func, ulp_tolerance, all_rounding); \
+ } \
+ TEST_F(LlvmLibcTanhTest##suffix, InDoubleRange) { \
+ test_in_range(&func, ulp_tolerance, all_rounding); \
+ } \
+ static_assert(true, "Require semicolon.")
+
+LIST_TANH_TESTS(Default, LIBC_NAMESPACE::tanh,
+ /*ulp_tolerance=*/TOLERANCE + 0.5,
+ /*all_rounding=*/ALL_ROUNDING);
+LIST_TANH_TESTS(DoubleEval, tanh_double_eval,
+ /*ulp_tolerance=*/TOLERANCE + 0.5,
+ /*all_rounding=*/ALL_ROUNDING);
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 83cd98d6422cbc..a2713132f5c45f 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -12073,6 +12073,57 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_tanh_accurate_128",
+ hdrs = ["src/__support/math/tanh_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_tanh_double_eval",
+ hdrs = ["src/__support/math/tanh_double_eval.h"],
+ deps = [
+ ":__support_common",
+ ":__support_cpp_bit",
+ ":__support_fputil_double_double",
+ ":__support_fputil_except_value_utils",
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_fp_bitwise",
+ ":__support_fputil_multiply_add",
+ ":__support_macros_attributes",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_macros_properties_cpu_features",
+ ":__support_math_exp_utils",
+ ":__support_math_tanh_accurate_128",
+ ":hdr_errno_macros",
+ ":hdr_fenv_macros",
+ ],
+)
+
+libc_support_library(
+ name = "__support_math_tanh",
+ hdrs = ["src/__support/math/tanh.h"],
+ deps = [
+ ":__support_common",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_math_tanh_double_eval",
+ ],
+)
+
libc_support_library(
name = "__support_math_tanhf",
hdrs = ["src/__support/math/tanhf.h"],
@@ -15612,6 +15663,11 @@ libc_math_function(
],
)
+libc_math_function(
+ name = "tanh",
+ additional_deps = [":__support_math_tanh"],
+)
+
libc_math_function(
name = "tanhf",
additional_deps = [
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 46cf02bf9e7271..9f2369d6f11aa2 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
@@ -1044,6 +1044,14 @@ math_mpfr_test(
],
)
+math_mpfr_test(
+ name = "tanh",
+ deps = [
+ "//libc:__support_macros_attributes",
+ "//libc:__support_math_tanh_double_eval",
+ ],
+)
+
math_mpfr_test(name = "tanhf")
math_mpfr_test(
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 b47989fcb6b5fb..1c2c5629789c73 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
@@ -2177,6 +2177,15 @@ math_test(name = "tan")
math_test(name = "tanf")
+math_test(
+ name = "tanh",
+ deps = [
+ "//libc:__support_macros_attributes",
+ "//libc:__support_macros_optimization",
+ "//libc:__support_math_tanh_double_eval",
+ ],
+)
+
math_test(name = "tanhf")
math_test(name = "tanpif")
More information about the llvm-branch-commits
mailing list