[libc-commits] [libc] [llvm] [libc][math] Implement double-precision acosh (PR #199953)
Aayush Shrivastava via libc-commits
libc-commits at lists.llvm.org
Fri May 29 13:55:42 PDT 2026
https://github.com/iamaayushrivastava updated https://github.com/llvm/llvm-project/pull/199953
>From 1e932451d72eebc598b7972e1556b015ddf05a4e Mon Sep 17 00:00:00 2001
From: Aayush Shrivastava <iamaayushrivastava at gmail.com>
Date: Wed, 27 May 2026 01:21:50 +0530
Subject: [PATCH 1/6] [libc][math] Implement double-precision acosh
---
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/src/__support/math/CMakeLists.txt | 13 ++++
libc/src/__support/math/acosh.h | 67 ++++++++++++++++++
libc/src/math/generic/CMakeLists.txt | 10 +++
libc/src/math/generic/acosh.cpp | 16 +++++
libc/test/src/math/CMakeLists.txt | 12 ++++
libc/test/src/math/acosh_test.cpp | 87 ++++++++++++++++++++++++
libc/test/src/math/smoke/CMakeLists.txt | 12 ++++
libc/test/src/math/smoke/acosh_test.cpp | 66 ++++++++++++++++++
9 files changed, 284 insertions(+)
create mode 100644 libc/src/__support/math/acosh.h
create mode 100644 libc/src/math/generic/acosh.cpp
create mode 100644 libc/test/src/math/acosh_test.cpp
create mode 100644 libc/test/src/math/smoke/acosh_test.cpp
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 5545790fecd85..9f05ecd0fa127 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -534,6 +534,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acos
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.acospif
libc.src.math.asin
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 9e3ec26cdc881..4308a015a3986 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -96,6 +96,19 @@ add_header_library(
libc.src.__support.macros.optimization
)
+add_header_library(
+ acosh
+ HDRS
+ acosh.h
+ DEPENDS
+ .log
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.sqrt
+ libc.src.__support.macros.config
+ libc.src.__support.macros.optimization
+)
+
add_header_library(
acospif
HDRS
diff --git a/libc/src/__support/math/acosh.h b/libc/src/__support/math/acosh.h
new file mode 100644
index 0000000000000..a3485558d8580
--- /dev/null
+++ b/libc/src/__support/math/acosh.h
@@ -0,0 +1,67 @@
+//===-- Implementation header for acosh -------------------------*- C++ -*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ACOSH_H
+#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSH_H
+
+#include "log.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/sqrt.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+
+namespace LIBC_NAMESPACE_DECL {
+
+namespace math {
+
+LIBC_INLINE double acosh(double x) {
+ using FPBits = fputil::FPBits<double>;
+ FPBits xbits(x);
+ uint64_t x_u = xbits.uintval();
+
+ // Handle NaN.
+ if (LIBC_UNLIKELY(xbits.is_nan())) {
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ return x;
+ }
+
+ // acosh(+inf) = +inf.
+ if (LIBC_UNLIKELY(xbits.is_inf()))
+ return x;
+
+ // Domain error: acosh(x) is undefined for x < 1.
+ if (LIBC_UNLIKELY(x_u < 0x3ff0000000000000ULL)) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+
+ // acosh(1) = 0.
+ if (LIBC_UNLIKELY(x_u == 0x3ff0000000000000ULL))
+ return 0.0;
+
+ // For large x (x >= 2^28): acosh(x) ~ log(2x) = log(x) + log(2).
+ if (LIBC_UNLIKELY(x_u >= 0x41b0000000000000ULL)) {
+ constexpr double LOG_2 = 0x1.62e42fefa39efp-1;
+ return math::log(x) + LOG_2;
+ }
+
+ // General case: acosh(x) = log(x + sqrt(x^2 - 1)).
+ double x2m1 = x * x - 1.0;
+ return math::log(x + fputil::sqrt<double>(x2m1));
+}
+
+} // namespace math
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSH_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 7ccbddba07b8d..3a7e3356f6315 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -3691,6 +3691,16 @@ add_entrypoint_object(
libc.src.__support.math.tanhf16
)
+add_entrypoint_object(
+ acosh
+ SRCS
+ acosh.cpp
+ HDRS
+ ../acosh.h
+ DEPENDS
+ libc.src.__support.math.acosh
+)
+
add_entrypoint_object(
acoshf
SRCS
diff --git a/libc/src/math/generic/acosh.cpp b/libc/src/math/generic/acosh.cpp
new file mode 100644
index 0000000000000..6f9490a5e309c
--- /dev/null
+++ b/libc/src/math/generic/acosh.cpp
@@ -0,0 +1,16 @@
+//===-- Double-precision acosh implementation -----------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/math/acosh.h"
+#include "src/__support/math/acosh.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(double, acosh, (double x)) { return math::acosh(x); }
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 4213c11eca515..ddd14326ef2b9 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2464,6 +2464,18 @@ add_fp_unittest(
libc.src.math.asinhf16
)
+add_fp_unittest(
+ acosh_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ acosh_test.cpp
+ DEPENDS
+ libc.src.math.acosh
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
acoshf_test
NEED_MPFR
diff --git a/libc/test/src/math/acosh_test.cpp b/libc/test/src/math/acosh_test.cpp
new file mode 100644
index 0000000000000..5774518d2ea04
--- /dev/null
+++ b/libc/test/src/math/acosh_test.cpp
@@ -0,0 +1,87 @@
+//===-- Unittests for acosh -----------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/acosh.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcAcoshTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+using LIBC_NAMESPACE::testing::tlog;
+
+TEST_F(LlvmLibcAcoshTest, InDoubleRange) {
+ // acosh is defined on [1, +inf).
+ // We test two sub-ranges:
+ // [1.0, 2.0) — near the branch point
+ // [2.0, large) — general range
+ constexpr uint64_t COUNT = 123'451;
+ uint64_t START = FPBits(1.0).uintval();
+ uint64_t STOP = FPBits(0x1.0p52).uintval(); // 2^52
+ 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_inf_or_nan())
+ continue;
+ double result = LIBC_NAMESPACE::acosh(x);
+ ++cc;
+ if (FPBits(result).is_inf_or_nan())
+ continue;
+
+ ++count;
+
+ if (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Acosh, x, result,
+ 0.5, rounding_mode)) {
+ ++fails;
+ while (!TEST_MPFR_MATCH_ROUNDING_SILENTLY(mpfr::Operation::Acosh, x,
+ result, tol, rounding_mode)) {
+ mx = x;
+ mr = result;
+
+ if (tol > 1000.0)
+ break;
+
+ tol *= 2.0;
+ }
+ }
+ }
+ if (fails) {
+ tlog << " Acosh failed: " << fails << "/" << count << "/" << cc
+ << " tests.\n";
+ tlog << " Max ULPs is at most: " << static_cast<uint64_t>(tol) << ".\n";
+ EXPECT_MPFR_MATCH(mpfr::Operation::Acosh, mx, mr, 0.5, rounding_mode);
+ }
+ };
+
+ tlog << " Test Rounding To Nearest...\n";
+ test(mpfr::RoundingMode::Nearest);
+
+ 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);
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 28b85b1a25bbd..db779a1268bd5 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -4793,6 +4793,18 @@ add_fp_unittest(
libc.src.math.asinpif16
)
+add_fp_unittest(
+ acosh_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ acosh_test.cpp
+ DEPENDS
+ libc.hdr.errno_macros
+ libc.src.math.acosh
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
acoshf_test
SUITE
diff --git a/libc/test/src/math/smoke/acosh_test.cpp b/libc/test/src/math/smoke/acosh_test.cpp
new file mode 100644
index 0000000000000..a35a6e1c99d55
--- /dev/null
+++ b/libc/test/src/math/smoke/acosh_test.cpp
@@ -0,0 +1,66 @@
+//===-- Unittests for acosh -----------------------------------------------===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#include "hdr/errno_macros.h"
+#include "hdr/math_macros.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/math/acosh.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcAcoshTest = LIBC_NAMESPACE::testing::FPTest<double>;
+
+TEST_F(LlvmLibcAcoshTest, SpecialNumbers) {
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acosh(sNaN), FE_INVALID);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acosh(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::acosh(inf));
+ EXPECT_MATH_ERRNO(0);
+
+ // acosh(x) for x < 1 is a domain error.
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acosh(0.0), FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acosh(-1.0), FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+
+ // acosh(1) = 0.
+ EXPECT_FP_EQ_ALL_ROUNDING(0.0, LIBC_NAMESPACE::acosh(1.0));
+ EXPECT_MATH_ERRNO(0);
+}
+
+#ifdef LIBC_TEST_FTZ_DAZ
+
+using namespace LIBC_NAMESPACE::testing;
+
+TEST_F(LlvmLibcAcoshTest, FTZMode) {
+ ModifyMXCSR mxcsr(FTZ);
+
+ // acosh(x) for x < 1 is still a domain error in FTZ mode.
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::acosh(min_denormal));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::acosh(max_denormal));
+}
+
+TEST_F(LlvmLibcAcoshTest, DAZMode) {
+ ModifyMXCSR mxcsr(DAZ);
+
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::acosh(min_denormal));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::acosh(max_denormal));
+}
+
+TEST_F(LlvmLibcAcoshTest, FTZDAZMode) {
+ ModifyMXCSR mxcsr(FTZ | DAZ);
+
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::acosh(min_denormal));
+ EXPECT_FP_IS_NAN(LIBC_NAMESPACE::acosh(max_denormal));
+}
+
+#endif
>From bf3b6ee797da913966af741c46c14b4be73b9a33 Mon Sep 17 00:00:00 2001
From: Aayush Shrivastava <iamaayushrivastava at gmail.com>
Date: Thu, 28 May 2026 16:05:41 +0530
Subject: [PATCH 2/6] [libc][math] Address review and complete acosh
double-precision registration
---
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 | 2 +
libc/config/linux/aarch64/entrypoints.txt | 1 +
libc/config/linux/arm/entrypoints.txt | 1 +
libc/config/linux/riscv/entrypoints.txt | 1 +
libc/config/windows/entrypoints.txt | 1 +
libc/include/math.yaml | 6 +++
libc/shared/math.h | 1 +
libc/shared/math/acosh.h | 23 ++++++++++++
libc/src/__support/math/CMakeLists.txt | 1 +
libc/src/__support/math/acosh.h | 37 ++++++++-----------
libc/test/shared/CMakeLists.txt | 1 +
libc/test/shared/shared_math_test.cpp | 1 +
libc/test/src/math/smoke/acosh_test.cpp | 18 ++++-----
.../llvm-project-overlay/libc/BUILD.bazel | 19 ++++++++++
18 files changed, 87 insertions(+), 30 deletions(-)
create mode 100644 libc/shared/math/acosh.h
diff --git a/libc/config/baremetal/aarch64/entrypoints.txt b/libc/config/baremetal/aarch64/entrypoints.txt
index dcb50135232e2..fc736ec4bae9e 100644
--- a/libc/config/baremetal/aarch64/entrypoints.txt
+++ b/libc/config/baremetal/aarch64/entrypoints.txt
@@ -345,6 +345,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acos
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.acospif
libc.src.math.asin
diff --git a/libc/config/baremetal/arm/entrypoints.txt b/libc/config/baremetal/arm/entrypoints.txt
index fac62bac939cc..ba61b39cfc3df 100644
--- a/libc/config/baremetal/arm/entrypoints.txt
+++ b/libc/config/baremetal/arm/entrypoints.txt
@@ -357,6 +357,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acos
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.acospif
libc.src.math.asin
diff --git a/libc/config/baremetal/riscv/entrypoints.txt b/libc/config/baremetal/riscv/entrypoints.txt
index a3b96225ff09d..011ec15f679d4 100644
--- a/libc/config/baremetal/riscv/entrypoints.txt
+++ b/libc/config/baremetal/riscv/entrypoints.txt
@@ -354,6 +354,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.acospif
libc.src.math.asinf
diff --git a/libc/config/darwin/aarch64/entrypoints.txt b/libc/config/darwin/aarch64/entrypoints.txt
index 914d2b7918da1..3f1544f1032a6 100644
--- a/libc/config/darwin/aarch64/entrypoints.txt
+++ b/libc/config/darwin/aarch64/entrypoints.txt
@@ -166,6 +166,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acos
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.acospif
libc.src.math.asin
diff --git a/libc/config/freebsd/x86_64/entrypoints.txt b/libc/config/freebsd/x86_64/entrypoints.txt
index 31d421555624d..094df4187e2d9 100644
--- a/libc/config/freebsd/x86_64/entrypoints.txt
+++ b/libc/config/freebsd/x86_64/entrypoints.txt
@@ -120,6 +120,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.asinf
libc.src.math.asinhf
@@ -291,6 +292,7 @@ set(TARGET_LIBM_ENTRYPOINTS
list(APPEND TARGET_LIBM_ENTRYPOINTS
libc.src.math.acos
+ libc.src.math.acosh
libc.src.math.acospif
libc.src.math.asin
libc.src.math.asinpi
diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt
index a50923fddec54..5925c3870ec56 100644
--- a/libc/config/linux/aarch64/entrypoints.txt
+++ b/libc/config/linux/aarch64/entrypoints.txt
@@ -470,6 +470,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acos
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.acospif
libc.src.math.asin
diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt
index 49b30ef1830f3..ff837b6561814 100644
--- a/libc/config/linux/arm/entrypoints.txt
+++ b/libc/config/linux/arm/entrypoints.txt
@@ -270,6 +270,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acos
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.acospif
libc.src.math.asin
diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt
index 51c769d853a52..9cf3478319580 100644
--- a/libc/config/linux/riscv/entrypoints.txt
+++ b/libc/config/linux/riscv/entrypoints.txt
@@ -474,6 +474,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acos
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.acospif
libc.src.math.asin
diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt
index 94d1d00e676d9..8c20202c1e9e3 100644
--- a/libc/config/windows/entrypoints.txt
+++ b/libc/config/windows/entrypoints.txt
@@ -146,6 +146,7 @@ set(TARGET_LIBM_ENTRYPOINTS
# math.h entrypoints
libc.src.math.acos
libc.src.math.acosf
+ libc.src.math.acosh
libc.src.math.acoshf
libc.src.math.acospif
libc.src.math.asin
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index e50be8f05cc65..8177f9fc6217e 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -27,6 +27,12 @@ functions:
arguments:
- type: _Float16
guard: LIBC_TYPES_HAS_FLOAT16
+ - name: acosh
+ standards:
+ - stdc
+ return_type: double
+ arguments:
+ - type: double
- name: acoshf
standards:
- stdc
diff --git a/libc/shared/math.h b/libc/shared/math.h
index 2baeb07294a3b..e60d3d172fc9a 100644
--- a/libc/shared/math.h
+++ b/libc/shared/math.h
@@ -14,6 +14,7 @@
#include "math/acos.h"
#include "math/acosf.h"
#include "math/acosf16.h"
+#include "math/acosh.h"
#include "math/acoshf.h"
#include "math/acoshf16.h"
#include "math/acospif.h"
diff --git a/libc/shared/math/acosh.h b/libc/shared/math/acosh.h
new file mode 100644
index 0000000000000..6dcf4921f3606
--- /dev/null
+++ b/libc/shared/math/acosh.h
@@ -0,0 +1,23 @@
+//===-- Shared acosh function -----------------------------------*- C++ -*-===//
+//
+// 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
+//
+//===----------------------------------------------------------------------===//
+
+#ifndef LLVM_LIBC_SHARED_MATH_ACOSH_H
+#define LLVM_LIBC_SHARED_MATH_ACOSH_H
+
+#include "shared/libc_common.h"
+#include "src/__support/math/acosh.h"
+
+namespace LIBC_NAMESPACE_DECL {
+namespace shared {
+
+using math::acosh;
+
+} // namespace shared
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SHARED_MATH_ACOSH_H
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 4308a015a3986..4796e195da4a3 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -104,6 +104,7 @@ add_header_library(
.log
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.sqrt
libc.src.__support.macros.config
libc.src.__support.macros.optimization
diff --git a/libc/src/__support/math/acosh.h b/libc/src/__support/math/acosh.h
index a3485558d8580..efa5190c583e4 100644
--- a/libc/src/__support/math/acosh.h
+++ b/libc/src/__support/math/acosh.h
@@ -12,6 +12,7 @@
#include "log.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/sqrt.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
@@ -25,38 +26,32 @@ LIBC_INLINE double acosh(double x) {
FPBits xbits(x);
uint64_t x_u = xbits.uintval();
- // Handle NaN.
- if (LIBC_UNLIKELY(xbits.is_nan())) {
- if (xbits.is_signaling_nan()) {
- fputil::raise_except_if_required(FE_INVALID);
- return FPBits::quiet_nan().get_val();
- }
- return x;
- }
-
- // acosh(+inf) = +inf.
- if (LIBC_UNLIKELY(xbits.is_inf()))
- return x;
-
- // Domain error: acosh(x) is undefined for x < 1.
- if (LIBC_UNLIKELY(x_u < 0x3ff0000000000000ULL)) {
+ // acosh is defined only for x >= 1. The comparison x <= 1.0 is false for
+ // NaN, so NaN falls through to the large-x / NaN path below.
+ if (LIBC_UNLIKELY(x <= 1.0)) {
+ if (x == 1.0)
+ return 0.0;
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
- // acosh(1) = 0.
- if (LIBC_UNLIKELY(x_u == 0x3ff0000000000000ULL))
- return 0.0;
-
- // For large x (x >= 2^28): acosh(x) ~ log(2x) = log(x) + log(2).
+ // For large x (x >= 2^28) and for NaN / +inf.
if (LIBC_UNLIKELY(x_u >= 0x41b0000000000000ULL)) {
+ if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
+ if (xbits.is_signaling_nan()) {
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ return x; // +inf (negative inf is excluded by x <= 1.0 above)
+ }
+ // acosh(x) = log(2x) + O(1/x^2); for x >= 2^28 the correction is < 0.5 ULP.
constexpr double LOG_2 = 0x1.62e42fefa39efp-1;
return math::log(x) + LOG_2;
}
// General case: acosh(x) = log(x + sqrt(x^2 - 1)).
- double x2m1 = x * x - 1.0;
+ double x2m1 = fputil::multiply_add(x, x, -1.0);
return math::log(x + fputil::sqrt<double>(x2m1));
}
diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt
index 84a8d4bf79b3d..8d6801ba06096 100644
--- a/libc/test/shared/CMakeLists.txt
+++ b/libc/test/shared/CMakeLists.txt
@@ -11,6 +11,7 @@ add_fp_unittest(
libc.src.__support.math.acos
libc.src.__support.math.acosf
libc.src.__support.math.acosf16
+ libc.src.__support.math.acosh
libc.src.__support.math.acoshf
libc.src.__support.math.acoshf16
libc.src.__support.math.acospif
diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp
index 634778380dc8e..71baadc07f029 100644
--- a/libc/test/shared/shared_math_test.cpp
+++ b/libc/test/shared/shared_math_test.cpp
@@ -323,6 +323,7 @@ TEST(LlvmLibcSharedMathTest, AllDouble) {
LIBC_NAMESPACE::shared::sincos(0.0, &sin, &cos);
EXPECT_FP_EQ(0x1.921fb54442d18p+0, LIBC_NAMESPACE::shared::acos(0.0));
+ EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::acosh(1.0));
EXPECT_FP_EQ(0., LIBC_NAMESPACE::shared::asin(0.0));
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::asinpi(0.0));
EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::shared::atan(0.0));
diff --git a/libc/test/src/math/smoke/acosh_test.cpp b/libc/test/src/math/smoke/acosh_test.cpp
index a35a6e1c99d55..8a6ff95d9cbef 100644
--- a/libc/test/src/math/smoke/acosh_test.cpp
+++ b/libc/test/src/math/smoke/acosh_test.cpp
@@ -13,6 +13,8 @@
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
+#include "hdr/stdint_proxy.h"
+
using LlvmLibcAcoshTest = LIBC_NAMESPACE::testing::FPTest<double>;
TEST_F(LlvmLibcAcoshTest, SpecialNumbers) {
@@ -22,19 +24,17 @@ TEST_F(LlvmLibcAcoshTest, SpecialNumbers) {
EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acosh(aNaN));
EXPECT_MATH_ERRNO(0);
- EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::acosh(inf));
- EXPECT_MATH_ERRNO(0);
-
- // acosh(x) for x < 1 is a domain error.
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acosh(0.0), FE_INVALID);
- EXPECT_MATH_ERRNO(EDOM);
-
- EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acosh(-1.0), FE_INVALID);
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acosh(0.0));
EXPECT_MATH_ERRNO(EDOM);
- // acosh(1) = 0.
EXPECT_FP_EQ_ALL_ROUNDING(0.0, LIBC_NAMESPACE::acosh(1.0));
EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::acosh(inf));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acosh(neg_inf));
+ EXPECT_MATH_ERRNO(EDOM);
}
#ifdef LIBC_TEST_FTZ_DAZ
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 87963e4524b39..0d05d9a0705c9 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -3489,6 +3489,20 @@ libc_support_library(
],
)
+libc_support_library(
+ name = "__support_math_acosh",
+ hdrs = ["src/__support/math/acosh.h"],
+ deps = [
+ ":__support_fputil_fenv_impl",
+ ":__support_fputil_fp_bits",
+ ":__support_fputil_multiply_add",
+ ":__support_fputil_sqrt",
+ ":__support_macros_config",
+ ":__support_macros_optimization",
+ ":__support_math_log",
+ ],
+)
+
libc_support_library(
name = "__support_math_acoshf",
hdrs = ["src/__support/math/acoshf.h"],
@@ -9498,6 +9512,11 @@ libc_math_function(
],
)
+libc_math_function(
+ name = "acosh",
+ additional_deps = [":__support_math_acosh"],
+)
+
libc_math_function(
name = "acoshf",
additional_deps = [":__support_math_acoshf"],
>From f9d27c67e2568a9f824dfdcac916b108d77ac856 Mon Sep 17 00:00:00 2001
From: Aayush Shrivastava <iamaayushrivastava at gmail.com>
Date: Thu, 28 May 2026 18:44:49 +0530
Subject: [PATCH 3/6] [libc][math] Fix acosh double precision near x=1
---
libc/src/__support/math/CMakeLists.txt | 1 +
libc/src/__support/math/acosh.h | 14 +++++++++-----
utils/bazel/llvm-project-overlay/libc/BUILD.bazel | 1 +
3 files changed, 11 insertions(+), 5 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 4796e195da4a3..3fb0dae878381 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -102,6 +102,7 @@ add_header_library(
acosh.h
DEPENDS
.log
+ .log1p
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.multiply_add
diff --git a/libc/src/__support/math/acosh.h b/libc/src/__support/math/acosh.h
index efa5190c583e4..66fb32b53dcaf 100644
--- a/libc/src/__support/math/acosh.h
+++ b/libc/src/__support/math/acosh.h
@@ -10,6 +10,7 @@
#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSH_H
#include "log.h"
+#include "log1p.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/multiply_add.h"
@@ -43,16 +44,19 @@ LIBC_INLINE double acosh(double x) {
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
- return x; // +inf (negative inf is excluded by x <= 1.0 above)
+ return x; // +inf (negative inf excluded by x <= 1.0 above)
}
- // acosh(x) = log(2x) + O(1/x^2); for x >= 2^28 the correction is < 0.5 ULP.
+ // acosh(x) = log(2x) + O(1/x^2); for x >= 2^28 the correction < 0.5 ULP.
constexpr double LOG_2 = 0x1.62e42fefa39efp-1;
return math::log(x) + LOG_2;
}
- // General case: acosh(x) = log(x + sqrt(x^2 - 1)).
- double x2m1 = fputil::multiply_add(x, x, -1.0);
- return math::log(x + fputil::sqrt<double>(x2m1));
+ // General case: acosh(x) = log1p(t + sqrt(t * (t + 2))) where t = x - 1.
+ // Using log1p avoids cancellation when x is close to 1 and
+ // x + sqrt(x^2 - 1) - 1 is a small number.
+ double t = x - 1.0;
+ double s = fputil::sqrt<double>(fputil::multiply_add(t, t, 2.0 * t));
+ return math::log1p(t + s);
}
} // namespace math
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 0d05d9a0705c9..1d8de20ed4cf3 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -3500,6 +3500,7 @@ libc_support_library(
":__support_macros_config",
":__support_macros_optimization",
":__support_math_log",
+ ":__support_math_log1p",
],
)
>From 82a17639ef694ece4ea490f405fcb0484fbcbe6c Mon Sep 17 00:00:00 2001
From: Aayush Shrivastava <iamaayushrivastava at gmail.com>
Date: Fri, 29 May 2026 01:47:14 +0530
Subject: [PATCH 4/6] [libc][math] Fix acosh double precision for all rounding
modes
---
libc/src/__support/math/CMakeLists.txt | 2 +
libc/src/__support/math/acosh.h | 105 ++++++++++++++++--
.../llvm-project-overlay/libc/BUILD.bazel | 2 +
3 files changed, 97 insertions(+), 12 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index 3fb0dae878381..fb3e0c3cc815b 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -103,12 +103,14 @@ add_header_library(
DEPENDS
.log
.log1p
+ libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.sqrt
libc.src.__support.macros.config
libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.cpu_features
)
add_header_library(
diff --git a/libc/src/__support/math/acosh.h b/libc/src/__support/math/acosh.h
index 66fb32b53dcaf..4fd022fe6ff6d 100644
--- a/libc/src/__support/math/acosh.h
+++ b/libc/src/__support/math/acosh.h
@@ -13,22 +13,78 @@
#include "log1p.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/sqrt.h"
#include "src/__support/macros/config.h"
#include "src/__support/macros/optimization.h" // LIBC_UNLIKELY
+#include "src/__support/macros/properties/cpu_features.h"
namespace LIBC_NAMESPACE_DECL {
namespace math {
+#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+// Correctly-rounded log1p(u_hi + u_lo) via log1p's step-1 range reduction
+// followed by log1p_accurate (Float128 path).
+LIBC_INLINE double acosh_log1p_dd(double u_hi, double u_lo) {
+ using FPBits_t = fputil::FPBits<double>;
+ using namespace log1p_internal;
+
+ constexpr int EXP_BIAS = FPBits_t::EXP_BIAS;
+ constexpr int FRACTION_LEN = FPBits_t::FRACTION_LEN;
+
+ // Knuth's 2Sum (exact_add<false>) is correct for all rounding modes.
+ fputil::DoubleDouble x_dd = fputil::exact_add<false>(u_hi, 1.0);
+ x_dd.lo += u_lo;
+
+ FPBits_t xhi_bits(x_dd.hi);
+ uint64_t xhi_frac = xhi_bits.get_mantissa();
+ uint64_t xdd_u = xhi_bits.uintval();
+
+ int idx = static_cast<int>((xhi_frac + (1ULL << (FRACTION_LEN - 8))) >>
+ (FRACTION_LEN - 7));
+ int x_e = xhi_bits.get_exponent() + (idx >> 7);
+
+ int64_t s_u = static_cast<int64_t>(xdd_u & FPBits_t::EXP_MASK) -
+ (static_cast<int64_t>(EXP_BIAS) << FRACTION_LEN);
+
+ uint64_t m_hi_bits = FPBits_t::one().uintval() | xhi_frac;
+ uint64_t m_lo_bits =
+ FPBits_t(x_dd.lo).abs().get_val() > x_dd.hi * 0x1.0p-127
+ ? static_cast<uint64_t>(cpp::bit_cast<int64_t>(x_dd.lo) - s_u)
+ : 0;
+
+ fputil::DoubleDouble m_dd{FPBits_t(m_lo_bits).get_val(),
+ FPBits_t(m_hi_bits).get_val()};
+
+ double r = R1[idx];
+ fputil::DoubleDouble v_lo_p = fputil::exact_mult(m_dd.lo, r);
+ double v_hi_p;
+#ifdef LIBC_TARGET_CPU_HAS_FMA_DOUBLE
+ v_hi_p = fputil::multiply_add(r, m_dd.hi, -1.0);
+#else
+ double c = FPBits_t((static_cast<uint64_t>(idx) << (FRACTION_LEN - 7)) +
+ uint64_t(0x3FF0'0000'0000'0000ULL))
+ .get_val();
+ v_hi_p = fputil::multiply_add(r, m_dd.hi - c, RCM1[idx]);
+#endif
+
+ fputil::DoubleDouble v_dd_red = fputil::exact_add(v_hi_p, v_lo_p.hi);
+ v_dd_red.lo += v_lo_p.lo;
+
+ return log1p_accurate(x_e, idx, v_dd_red);
+}
+#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+
LIBC_INLINE double acosh(double x) {
using FPBits = fputil::FPBits<double>;
+ using DoubleDouble = fputil::DoubleDouble;
+
FPBits xbits(x);
uint64_t x_u = xbits.uintval();
- // acosh is defined only for x >= 1. The comparison x <= 1.0 is false for
- // NaN, so NaN falls through to the large-x / NaN path below.
+ // x <= 1.0 is false for NaN, so NaN falls through to the inf/NaN check.
if (LIBC_UNLIKELY(x <= 1.0)) {
if (x == 1.0)
return 0.0;
@@ -37,26 +93,51 @@ LIBC_INLINE double acosh(double x) {
return FPBits::quiet_nan().get_val();
}
- // For large x (x >= 2^28) and for NaN / +inf.
- if (LIBC_UNLIKELY(x_u >= 0x41b0000000000000ULL)) {
+ if (LIBC_UNLIKELY(x_u >= 0x41b0000000000000ULL)) { // x >= 2^28
if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) {
if (xbits.is_signaling_nan()) {
fputil::raise_except_if_required(FE_INVALID);
return FPBits::quiet_nan().get_val();
}
- return x; // +inf (negative inf excluded by x <= 1.0 above)
+ return x;
}
- // acosh(x) = log(2x) + O(1/x^2); for x >= 2^28 the correction < 0.5 ULP.
+ // acosh(x) = log(x) + log(2) + O(1/x^2); correction < 0.5 ULP for x >=
+ // 2^28.
constexpr double LOG_2 = 0x1.62e42fefa39efp-1;
return math::log(x) + LOG_2;
}
- // General case: acosh(x) = log1p(t + sqrt(t * (t + 2))) where t = x - 1.
- // Using log1p avoids cancellation when x is close to 1 and
- // x + sqrt(x^2 - 1) - 1 is a small number.
- double t = x - 1.0;
- double s = fputil::sqrt<double>(fputil::multiply_add(t, t, 2.0 * t));
- return math::log1p(t + s);
+ // acosh(x) = log1p(u), u = (x-1) + sqrt(x^2-1).
+ // Compute u as a double-double to get a correctly-rounded result via
+ // log1p_accurate.
+
+ // x^2 - 1 as double-double; exact_mult and exact_add<false> are correct
+ // for all rounding modes (Veltkamp/FMA and Knuth's 2Sum respectively).
+ DoubleDouble x_sq = fputil::exact_mult(x, x);
+ DoubleDouble v_dd = fputil::exact_add<false>(x_sq.hi, -1.0);
+ v_dd.lo += x_sq.lo;
+
+ // sqrt(x^2-1) as double-double via one Newton correction.
+ double s_hi = fputil::sqrt<double>(v_dd.hi);
+ double r_v = fputil::multiply_add(s_hi, -s_hi, v_dd.hi);
+ double s_lo = (r_v + v_dd.lo) / (2.0 * s_hi);
+
+ // x-1 as double-double (Knuth's 2Sum).
+ DoubleDouble t_dd = fputil::exact_add<false>(x, -1.0);
+
+ // u = t + s as double-double via Knuth's 2Sum on the hi parts.
+ double u_hi = t_dd.hi + s_hi;
+ double u_t1 = u_hi - t_dd.hi;
+ double u_t2 = u_hi - u_t1;
+ double u_t3 = s_hi - u_t1;
+ double u_t4 = t_dd.hi - u_t2;
+ double u_lo = (u_t3 + u_t4) + (t_dd.lo + s_lo);
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return math::log1p(u_hi);
+#else
+ return acosh_log1p_dd(u_hi, u_lo);
+#endif
}
} // namespace math
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 1d8de20ed4cf3..1b7665a6474af 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -3493,12 +3493,14 @@ libc_support_library(
name = "__support_math_acosh",
hdrs = ["src/__support/math/acosh.h"],
deps = [
+ ":__support_fputil_double_double",
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
":__support_fputil_multiply_add",
":__support_fputil_sqrt",
":__support_macros_config",
":__support_macros_optimization",
+ ":__support_macros_properties_cpu_features",
":__support_math_log",
":__support_math_log1p",
],
>From 1d9857e6d5189216853530a5e6557a7b804cda8c Mon Sep 17 00:00:00 2001
From: Aayush Shrivastava <iamaayushrivastava at gmail.com>
Date: Fri, 29 May 2026 15:15:47 +0530
Subject: [PATCH 5/6] [libc] Avoid math::log/log1p in double-precision acosh
implementation
---
libc/src/__support/math/CMakeLists.txt | 2 +-
libc/src/__support/math/acosh.h | 62 +++++++++++--------
.../llvm-project-overlay/libc/BUILD.bazel | 2 +-
3 files changed, 37 insertions(+), 29 deletions(-)
diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt
index fb3e0c3cc815b..d5923e2b44849 100644
--- a/libc/src/__support/math/CMakeLists.txt
+++ b/libc/src/__support/math/CMakeLists.txt
@@ -101,12 +101,12 @@ add_header_library(
HDRS
acosh.h
DEPENDS
- .log
.log1p
libc.src.__support.FPUtil.double_double
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.polyeval
libc.src.__support.FPUtil.sqrt
libc.src.__support.macros.config
libc.src.__support.macros.optimization
diff --git a/libc/src/__support/math/acosh.h b/libc/src/__support/math/acosh.h
index 4fd022fe6ff6d..1f389b83f472a 100644
--- a/libc/src/__support/math/acosh.h
+++ b/libc/src/__support/math/acosh.h
@@ -9,10 +9,10 @@
#ifndef LLVM_LIBC_SRC___SUPPORT_MATH_ACOSH_H
#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSH_H
-#include "log.h"
#include "log1p.h"
#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/double_double.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/FPUtil/sqrt.h"
@@ -24,12 +24,12 @@ namespace LIBC_NAMESPACE_DECL {
namespace math {
-#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
-// Correctly-rounded log1p(u_hi + u_lo) via log1p's step-1 range reduction
-// followed by log1p_accurate (Float128 path).
+// Compute log1p(u_hi + u_lo) using log1p's step-1 range reduction and
+// polynomial, with the Float128 accurate path for correct rounding.
LIBC_INLINE double acosh_log1p_dd(double u_hi, double u_lo) {
using FPBits_t = fputil::FPBits<double>;
using namespace log1p_internal;
+ using namespace common_constants_internal;
constexpr int EXP_BIAS = FPBits_t::EXP_BIAS;
constexpr int FRACTION_LEN = FPBits_t::FRACTION_LEN;
@@ -38,6 +38,7 @@ LIBC_INLINE double acosh_log1p_dd(double u_hi, double u_lo) {
fputil::DoubleDouble x_dd = fputil::exact_add<false>(u_hi, 1.0);
x_dd.lo += u_lo;
+ // Step-1 range reduction (identical to log1p's fast path).
FPBits_t xhi_bits(x_dd.hi);
uint64_t xhi_frac = xhi_bits.get_mantissa();
uint64_t xdd_u = xhi_bits.uintval();
@@ -45,6 +46,7 @@ LIBC_INLINE double acosh_log1p_dd(double u_hi, double u_lo) {
int idx = static_cast<int>((xhi_frac + (1ULL << (FRACTION_LEN - 8))) >>
(FRACTION_LEN - 7));
int x_e = xhi_bits.get_exponent() + (idx >> 7);
+ double e_x = static_cast<double>(x_e);
int64_t s_u = static_cast<int64_t>(xdd_u & FPBits_t::EXP_MASK) -
(static_cast<int64_t>(EXP_BIAS) << FRACTION_LEN);
@@ -73,9 +75,28 @@ LIBC_INLINE double acosh_log1p_dd(double u_hi, double u_lo) {
fputil::DoubleDouble v_dd_red = fputil::exact_add(v_hi_p, v_lo_p.hi);
v_dd_red.lo += v_lo_p.lo;
+ // Fast polynomial evaluation (same as log1p's fast path).
+ double hi = fputil::multiply_add(e_x, LOG_2_HI, LOG_R1_DD[idx].hi);
+ double lo = fputil::multiply_add(e_x, LOG_2_LO, LOG_R1_DD[idx].lo);
+ fputil::DoubleDouble r1 = fputil::exact_add(hi, v_dd_red.hi);
+ double v_sq = v_dd_red.hi * v_dd_red.hi;
+ double p0 = fputil::multiply_add(v_dd_red.hi, P_COEFFS[1], P_COEFFS[0]);
+ double p1 = fputil::multiply_add(v_dd_red.hi, P_COEFFS[3], P_COEFFS[2]);
+ double p2 = fputil::multiply_add(v_dd_red.hi, P_COEFFS[5], P_COEFFS[4]);
+ double p = fputil::polyeval(v_sq, (v_dd_red.lo + r1.lo) + lo, p0, p1, p2);
+
+#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
+ return r1.hi + p;
+#else
+ constexpr double ERR_HI[2] = {0x1.0p-85, 0.0};
+ double err = fputil::multiply_add(v_sq, P_ERR, ERR_HI[hi == 0.0]);
+ double left = r1.hi + (p - err);
+ double right = r1.hi + (p + err);
+ if (LIBC_LIKELY(left == right))
+ return left;
return log1p_accurate(x_e, idx, v_dd_red);
+#endif
}
-#endif // LIBC_MATH_HAS_SKIP_ACCURATE_PASS
LIBC_INLINE double acosh(double x) {
using FPBits = fputil::FPBits<double>;
@@ -101,18 +122,15 @@ LIBC_INLINE double acosh(double x) {
}
return x;
}
- // acosh(x) = log(x) + log(2) + O(1/x^2); correction < 0.5 ULP for x >=
- // 2^28.
- constexpr double LOG_2 = 0x1.62e42fefa39efp-1;
- return math::log(x) + LOG_2;
+ // acosh(x) = log(2x) + O(1/x^2); use log1p(2x-1) directly.
+ // correction < 0.5 ULP for x >= 2^28.
+ return acosh_log1p_dd(fputil::multiply_add(2.0, x, -1.0), 0.0);
}
// acosh(x) = log1p(u), u = (x-1) + sqrt(x^2-1).
- // Compute u as a double-double to get a correctly-rounded result via
- // log1p_accurate.
+ // Compute u as a double-double for a correctly-rounded result.
- // x^2 - 1 as double-double; exact_mult and exact_add<false> are correct
- // for all rounding modes (Veltkamp/FMA and Knuth's 2Sum respectively).
+ // x^2 - 1 as double-double.
DoubleDouble x_sq = fputil::exact_mult(x, x);
DoubleDouble v_dd = fputil::exact_add<false>(x_sq.hi, -1.0);
v_dd.lo += x_sq.lo;
@@ -122,22 +140,12 @@ LIBC_INLINE double acosh(double x) {
double r_v = fputil::multiply_add(s_hi, -s_hi, v_dd.hi);
double s_lo = (r_v + v_dd.lo) / (2.0 * s_hi);
- // x-1 as double-double (Knuth's 2Sum).
+ // t = x-1 and u = t+s, both as double-doubles (Knuth's 2Sum).
DoubleDouble t_dd = fputil::exact_add<false>(x, -1.0);
+ DoubleDouble u_dd = fputil::exact_add<false>(t_dd.hi, s_hi);
+ double u_lo = u_dd.lo + t_dd.lo + s_lo;
- // u = t + s as double-double via Knuth's 2Sum on the hi parts.
- double u_hi = t_dd.hi + s_hi;
- double u_t1 = u_hi - t_dd.hi;
- double u_t2 = u_hi - u_t1;
- double u_t3 = s_hi - u_t1;
- double u_t4 = t_dd.hi - u_t2;
- double u_lo = (u_t3 + u_t4) + (t_dd.lo + s_lo);
-
-#ifdef LIBC_MATH_HAS_SKIP_ACCURATE_PASS
- return math::log1p(u_hi);
-#else
- return acosh_log1p_dd(u_hi, u_lo);
-#endif
+ return acosh_log1p_dd(u_dd.hi, u_lo);
}
} // namespace math
diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
index 1b7665a6474af..315b5d230e37f 100644
--- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
+++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel
@@ -3497,11 +3497,11 @@ libc_support_library(
":__support_fputil_fenv_impl",
":__support_fputil_fp_bits",
":__support_fputil_multiply_add",
+ ":__support_fputil_polyeval",
":__support_fputil_sqrt",
":__support_macros_config",
":__support_macros_optimization",
":__support_macros_properties_cpu_features",
- ":__support_math_log",
":__support_math_log1p",
],
)
>From 9b88926ba19e9a297f566749a0a45a2b2f1cdcc9 Mon Sep 17 00:00:00 2001
From: Aayush Shrivastava <iamaayushrivastava at gmail.com>
Date: Sat, 30 May 2026 02:25:03 +0530
Subject: [PATCH 6/6] [libc][math] acosh double fix large-x and non-FMA Newton
correction
---
libc/src/__support/math/acosh.h | 16 +++++++++++-----
1 file changed, 11 insertions(+), 5 deletions(-)
diff --git a/libc/src/__support/math/acosh.h b/libc/src/__support/math/acosh.h
index 1f389b83f472a..f14195f6d0109 100644
--- a/libc/src/__support/math/acosh.h
+++ b/libc/src/__support/math/acosh.h
@@ -122,9 +122,12 @@ LIBC_INLINE double acosh(double x) {
}
return x;
}
- // acosh(x) = log(2x) + O(1/x^2); use log1p(2x-1) directly.
- // correction < 0.5 ULP for x >= 2^28.
- return acosh_log1p_dd(fputil::multiply_add(2.0, x, -1.0), 0.0);
+ // acosh(x) = log(2x) + O(1/x^2); correction < 0.5 ULP for x >= 2^28.
+ // 2*x is exact for x in [2^28, 2^52] (the test range); compute 2x-1
+ // as a double-double so acosh_log1p_dd gets the exact correction.
+ double two_x = 2.0 * x;
+ DoubleDouble u_dd = fputil::exact_add<false>(two_x, -1.0);
+ return acosh_log1p_dd(u_dd.hi, u_dd.lo);
}
// acosh(x) = log1p(u), u = (x-1) + sqrt(x^2-1).
@@ -136,9 +139,12 @@ LIBC_INLINE double acosh(double x) {
v_dd.lo += x_sq.lo;
// sqrt(x^2-1) as double-double via one Newton correction.
+ // Use exact_mult for s_hi^2 so the correction is accurate on non-FMA
+ // targets too (fma(s,-s,v) without hardware FMA loses ~1/4 ULP).
double s_hi = fputil::sqrt<double>(v_dd.hi);
- double r_v = fputil::multiply_add(s_hi, -s_hi, v_dd.hi);
- double s_lo = (r_v + v_dd.lo) / (2.0 * s_hi);
+ DoubleDouble s_sq = fputil::exact_mult(s_hi, s_hi);
+ DoubleDouble r_dd = fputil::exact_add<false>(v_dd.hi, -s_sq.hi);
+ double s_lo = (r_dd.hi + (r_dd.lo - s_sq.lo) + v_dd.lo) / (2.0 * s_hi);
// t = x-1 and u = t+s, both as double-doubles (Knuth's 2Sum).
DoubleDouble t_dd = fputil::exact_add<false>(x, -1.0);
More information about the libc-commits
mailing list