[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:54:06 PDT 2026
llvmorg-github-actions[bot] wrote:
<!--LLVM PR SUMMARY COMMENT-->
@llvm/pr-subscribers-libc
@llvm/pr-subscribers-backend-amdgpu
Author: lntue
<details>
<summary>Changes</summary>
- 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.
---
Patch is 64.93 KiB, truncated to 20.00 KiB below, full version: https://github.com/llvm/llvm-project/pull/235618.diff
31 Files Affected:
- (modified) libc/config/baremetal/aarch64/entrypoints.txt (+1)
- (modified) libc/config/baremetal/arm/entrypoints.txt (+1)
- (modified) libc/config/baremetal/riscv/entrypoints.txt (+1)
- (modified) libc/config/darwin/aarch64/entrypoints.txt (+1)
- (modified) libc/config/freebsd/x86_64/entrypoints.txt (+1)
- (modified) libc/config/gpu/amdgpu/entrypoints.txt (+1)
- (modified) libc/config/gpu/nvptx/entrypoints.txt (+1)
- (modified) libc/config/linux/aarch64/entrypoints.txt (+1)
- (modified) libc/config/linux/arm/entrypoints.txt (+1)
- (modified) libc/config/linux/riscv/entrypoints.txt (+1)
- (modified) libc/config/linux/x86_64/entrypoints.txt (+1)
- (modified) libc/config/windows/entrypoints.txt (+1)
- (modified) libc/docs/headers/math/index.md (+1-1)
- (modified) libc/include/math.yaml (+6)
- (modified) libc/shared/math.h (+1)
- (added) libc/shared/math/tanh.h (+28)
- (modified) libc/src/__support/math/CMakeLists.txt (+51)
- (added) libc/src/__support/math/tanh.h (+34)
- (added) libc/src/__support/math/tanh_accurate_128.h (+221)
- (added) libc/src/__support/math/tanh_double_eval.h (+559)
- (modified) libc/src/math/generic/CMakeLists.txt (+10)
- (added) libc/src/math/generic/tanh.cpp (+21)
- (modified) libc/test/shared/CMakeLists.txt (+1)
- (modified) libc/test/shared/shared_math_test.cpp (+1)
- (modified) libc/test/src/math/CMakeLists.txt (+16)
- (modified) libc/test/src/math/smoke/CMakeLists.txt (+16)
- (added) libc/test/src/math/smoke/tanh_test.cpp (+146)
- (added) libc/test/src/math/tanh_test.cpp (+187)
- (modified) utils/bazel/llvm-project-overlay/libc/BUILD.bazel (+56)
- (modified) utils/bazel/llvm-project-overlay/libc/test/src/math/BUILD.bazel (+8)
- (modified) utils/bazel/llvm-project-overlay/libc/test/src/math/smoke/BUILD.bazel (+9)
``````````diff
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; ...
[truncated]
``````````
</details>
https://github.com/llvm/llvm-project/pull/235618
More information about the llvm-branch-commits
mailing list