[libc-commits] [libc] [libc][math][c23] Add atanhf16 C23 math function. (PR #132612)
Harrison Hao via libc-commits
libc-commits at lists.llvm.org
Sat Apr 12 03:52:30 PDT 2025
https://github.com/harrisonGPU updated https://github.com/llvm/llvm-project/pull/132612
>From 5f74e2fac4fc2e32fd1f01db40469253c6f353a1 Mon Sep 17 00:00:00 2001
From: Harrison Hao <tsworld1314 at gmail.com>
Date: Mon, 17 Mar 2025 15:34:24 +0000
Subject: [PATCH 1/8] [libc][math][c23] Add atanhf16 C23 math function.
---
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/docs/headers/math/index.rst | 2 +-
libc/include/math.yaml | 7 ++
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/atanhf16.h | 21 ++++++
libc/src/math/generic/CMakeLists.txt | 18 +++++
libc/src/math/generic/atanhf16.cpp | 86 ++++++++++++++++++++++
libc/test/src/math/CMakeLists.txt | 11 +++
libc/test/src/math/atanhf16_test.cpp | 39 ++++++++++
libc/test/src/math/smoke/CMakeLists.txt | 12 +++
libc/test/src/math/smoke/atanhf16_test.cpp | 58 +++++++++++++++
11 files changed, 255 insertions(+), 1 deletion(-)
create mode 100644 libc/src/math/atanhf16.h
create mode 100644 libc/src/math/generic/atanhf16.cpp
create mode 100644 libc/test/src/math/atanhf16_test.cpp
create mode 100644 libc/test/src/math/smoke/atanhf16_test.cpp
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 1ac3a781d5279..1215258cb6bc9 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -660,6 +660,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.acoshf16
libc.src.math.asinf16
libc.src.math.asinhf16
+ libc.src.math.atanhf16
libc.src.math.canonicalizef16
libc.src.math.ceilf16
libc.src.math.copysignf16
diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst
index 947bd4b60b391..dad10886adc09 100644
--- a/libc/docs/headers/math/index.rst
+++ b/libc/docs/headers/math/index.rst
@@ -267,7 +267,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| atan2pi | | | | | | 7.12.4.11 | F.10.1.11 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| atanh | |check| | | | | | 7.12.5.3 | F.10.2.3 |
+| atanh | |check| | | | |check| | | 7.12.5.3 | F.10.2.3 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| atanpi | | | | | | 7.12.4.10 | F.10.1.10 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/include/math.yaml b/libc/include/math.yaml
index fef829422244d..3f513ef40870b 100644
--- a/libc/include/math.yaml
+++ b/libc/include/math.yaml
@@ -113,6 +113,13 @@ functions:
return_type: float
arguments:
- type: float
+ - name: atanhf16
+ standards:
+ - stdc
+ return_type: _Float16
+ arguments:
+ - type: _Float16
+ guard: LIBC_TYPES_HAS_FLOAT16
- name: canonicalize
standards:
- stdc
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index d177ff79141c0..42ec04190d02c 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -66,6 +66,7 @@ add_math_entrypoint_object(atan2f128)
add_math_entrypoint_object(atanh)
add_math_entrypoint_object(atanhf)
+add_math_entrypoint_object(atanhf16)
add_math_entrypoint_object(canonicalize)
add_math_entrypoint_object(canonicalizef)
diff --git a/libc/src/math/atanhf16.h b/libc/src/math/atanhf16.h
new file mode 100644
index 0000000000000..9fbb262c16514
--- /dev/null
+++ b/libc/src/math/atanhf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for atanhf16 ----------------------*- 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_MATH_ATANHF16_H
+#define LLVM_LIBC_SRC_MATH_ATANHF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 atanhf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_ATANHF16_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index adbed5b2de48c..4440e7e28e2ef 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4029,6 +4029,24 @@ add_entrypoint_object(
libc.src.__support.macros.optimization
)
+add_entrypoint_object(
+ atanhf16
+ SRCS
+ atanhf16.cpp
+ HDRS
+ ../atanhf16.h
+ DEPENDS
+ .explogxf
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.polyeval
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
+)
+
add_object_library(
inv_trigf_utils
HDRS
diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp
new file mode 100644
index 0000000000000..94d6aa149cf00
--- /dev/null
+++ b/libc/src/math/generic/atanhf16.cpp
@@ -0,0 +1,86 @@
+//===-- Implementation of atanh(x) function -------------------------------===//
+//
+// 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/atanhf16.h"
+#include "explogxf.h"
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FPBits.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/common.h"
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/optimization.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
+ using FPBits = typename fputil::FPBits<float16>;
+
+ FPBits xbits(x);
+ Sign sign = xbits.sign();
+ uint16_t x_abs = xbits.abs().uintval();
+
+ if (LIBC_UNLIKELY(x_abs >= 0x3c00U)) {
+ if (xbits.is_nan()) {
+ return x;
+ }
+ // |x| == 1.0
+ if (x_abs == 0x3c00U) {
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_DIVBYZERO);
+ return FPBits::inf(sign).get_val();
+ } else {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
+ }
+ }
+
+ // For |x| less than approximately 0.10
+ if (LIBC_UNLIKELY(x_abs <= 0x2e66U)) {
+ // The Taylor expansion of atanh(x) is:
+ // atanh(x) = x + x^3/3 + x^5/5 + x^7/7 + x^9/9 + x^11/11
+ // = x * [1 + x^2/3 + x^4/5 + x^6/7 + x^8/9 + x^10/11]
+ // When |x| < 0x0100U, this can be approximated by:
+ // atanh(x) ≈ x + (1/3)*x^3
+ if (LIBC_UNLIKELY(x_abs < 0x0100U)) {
+ return static_cast<float16>(
+ LIBC_UNLIKELY(x_abs == 0) ? x : (x + 0x1.555556p-2 * x * x * x));
+ }
+
+ // For 0x0100U <= |x| <= 0x2e66U:
+ // Let t = x^2.
+ // Define P(t) ≈ (1/3)*t + (1/5)*t^2 + (1/7)*t^3 + (1/9)*t^4 + (1/11)*t^5.
+ // The coefficients below were derived using Sollya:
+ // > display = hexadecimal;
+ // > round(1/3, SG, RN);
+ // > round(1/5, SG, RN);
+ // > round(1/7, SG, RN);
+ // > round(1/9, SG, RN);
+ // > round(1/11, SG, RN);
+ // This yields:
+ // 0x1.555556p-2
+ // 0x1.99999ap-3
+ // 0x1.24924ap-3
+ // 0x1.c71c72p-4
+ // 0x1.745d18p-4f
+ // Thus, atanh(x) ≈ x * (1 + P(x^2)).
+ float xf = x;
+ float x2 = xf * xf;
+ float pe = fputil::polyeval(x2, 0.0f, 0x1.555556p-2f, 0x1.99999ap-3f,
+ 0x1.24924ap-3f, 0x1.c71c72p-4f, 0x1.745d18p-4f);
+ return static_cast<float16>(fputil::multiply_add(xf, pe, xf));
+ }
+
+ float xf = x;
+ return static_cast<float16>(0.5 * log_eval((xf + 1.0) / (xf - 1.0)));
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 7ee8b86135557..c80760f728335 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -2143,6 +2143,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ atanhf16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ atanhf16_test.cpp
+ DEPENDS
+ libc.src.math.atanhf16
+)
+
add_fp_unittest(
fmul_test
NEED_MPFR
diff --git a/libc/test/src/math/atanhf16_test.cpp b/libc/test/src/math/atanhf16_test.cpp
new file mode 100644
index 0000000000000..ce0179a1962df
--- /dev/null
+++ b/libc/test/src/math/atanhf16_test.cpp
@@ -0,0 +1,39 @@
+//===-- Unittests for atanhf16 --------------------------------------------===//
+//
+// 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/errno/libc_errno.h"
+#include "src/math/atanhf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+#include <stdint.h>
+
+using LlvmLibcAtanhf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+
+static constexpr uint16_t POS_START = 0x0000U;
+static constexpr uint16_t POS_STOP = 0x3BFFU;
+static constexpr uint16_t NEG_START = 0xBBFFU;
+static constexpr uint16_t NEG_STOP = 0x8000U;
+
+TEST_F(LlvmLibcAtanhf16Test, PositiveRange) {
+ for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atanh, x,
+ LIBC_NAMESPACE::atanhf16(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcAtanhf16Test, NegativeRange) {
+ for (uint16_t v = NEG_START; v <= NEG_STOP; --v) {
+ float16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atanh, x,
+ LIBC_NAMESPACE::atanhf16(x), 0.5);
+ }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 223d1933bca38..ec867339152e3 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3947,6 +3947,18 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ atanhf16_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ atanhf16_test.cpp
+ DEPENDS
+ libc.src.errno.errno
+ libc.src.math.atanhf16
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
asinhf_test
SUITE
diff --git a/libc/test/src/math/smoke/atanhf16_test.cpp b/libc/test/src/math/smoke/atanhf16_test.cpp
new file mode 100644
index 0000000000000..1ac483d2ae758
--- /dev/null
+++ b/libc/test/src/math/smoke/atanhf16_test.cpp
@@ -0,0 +1,58 @@
+//===-- Unittests for atanhf16 --------------------------------------------===//
+//
+// 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/cast.h"
+#include "src/errno/libc_errno.h"
+#include "src/math/atanhf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcAtanhf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcAtanhf16Test, SpecialNumbers) {
+ LIBC_NAMESPACE::libc_errno = 0;
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::atanhf16(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::atanhf16(zero));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(
+ -0.0f,
+ LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast<float16>(-0.0f)));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf,
+ LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast<float16>(1.0f)),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ neg_inf,
+ LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast<float16>(-1.0f)),
+ FE_DIVBYZERO);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_IS_NAN_WITH_EXCEPTION(
+ LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast<float16>(2.0f)),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_IS_NAN_WITH_EXCEPTION(
+ LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast<float16>(-2.0f)),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_IS_NAN_WITH_EXCEPTION(LIBC_NAMESPACE::atanhf16(inf), FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_IS_NAN_WITH_EXCEPTION(LIBC_NAMESPACE::atanhf16(neg_inf),
+ FE_INVALID);
+ EXPECT_MATH_ERRNO(EDOM);
+}
>From 6aa1640dc7f2530e6b9abd6ea600166c6b8fc8e5 Mon Sep 17 00:00:00 2001
From: Harrison Hao <tsworld1314 at gmail.com>
Date: Wed, 26 Mar 2025 10:18:23 +0800
Subject: [PATCH 2/8] [libc][math] Update range.
---
libc/test/src/math/atanhf16_test.cpp | 9 ++++++---
1 file changed, 6 insertions(+), 3 deletions(-)
diff --git a/libc/test/src/math/atanhf16_test.cpp b/libc/test/src/math/atanhf16_test.cpp
index ce0179a1962df..6637c7eead8f5 100644
--- a/libc/test/src/math/atanhf16_test.cpp
+++ b/libc/test/src/math/atanhf16_test.cpp
@@ -17,13 +17,16 @@
using LlvmLibcAtanhf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
+// Range for positive numbers: [0, 1)
static constexpr uint16_t POS_START = 0x0000U;
-static constexpr uint16_t POS_STOP = 0x3BFFU;
+static constexpr uint16_t POS_STOP = 0x3C00;
+
+// Range for negative numbers: (-1, 0]
static constexpr uint16_t NEG_START = 0xBBFFU;
static constexpr uint16_t NEG_STOP = 0x8000U;
TEST_F(LlvmLibcAtanhf16Test, PositiveRange) {
- for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+ for (uint16_t v = POS_START; v < POS_STOP; ++v) {
float16 x = FPBits(v).get_val();
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atanh, x,
LIBC_NAMESPACE::atanhf16(x), 0.5);
@@ -31,7 +34,7 @@ TEST_F(LlvmLibcAtanhf16Test, PositiveRange) {
}
TEST_F(LlvmLibcAtanhf16Test, NegativeRange) {
- for (uint16_t v = NEG_START; v <= NEG_STOP; --v) {
+ for (uint16_t v = NEG_START; v >= NEG_STOP; --v) {
float16 x = FPBits(v).get_val();
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atanh, x,
LIBC_NAMESPACE::atanhf16(x), 0.5);
>From afb32e244b1854166f956ce958724620d67a009d Mon Sep 17 00:00:00 2001
From: Harrison Hao <tsworld1314 at gmail.com>
Date: Sat, 29 Mar 2025 06:10:45 +0000
Subject: [PATCH 3/8] [libc][math][c23] Update for comments.
---
libc/src/math/generic/CMakeLists.txt | 1 +
libc/src/math/generic/atanhf16.cpp | 40 +++++++++-------------
libc/test/src/math/atanhf16_test.cpp | 13 ++++---
libc/test/src/math/smoke/atanhf16_test.cpp | 4 +--
4 files changed, 24 insertions(+), 34 deletions(-)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 4440e7e28e2ef..a83107f248562 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4040,6 +4040,7 @@ add_entrypoint_object(
libc.hdr.errno_macros
libc.hdr.fenv_macros
libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.multiply_add
libc.src.__support.FPUtil.polyeval
diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp
index 94d6aa149cf00..bb16a78a9862d 100644
--- a/libc/src/math/generic/atanhf16.cpp
+++ b/libc/src/math/generic/atanhf16.cpp
@@ -1,4 +1,4 @@
-//===-- Implementation of atanh(x) function -------------------------------===//
+//===-- Half-precision atanh(x) function ----------------------------------===//
//
// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
// See https://llvm.org/LICENSE.txt for license information.
@@ -10,6 +10,7 @@
#include "explogxf.h"
#include "hdr/errno_macros.h"
#include "hdr/fenv_macros.h"
+#include "src/__support/FPUtil/FEnvImpl.h"
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/cast.h"
@@ -21,26 +22,27 @@
namespace LIBC_NAMESPACE_DECL {
LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
- using FPBits = typename fputil::FPBits<float16>;
+ using FPBits = fputil::FPBits<float16>;
FPBits xbits(x);
Sign sign = xbits.sign();
uint16_t x_abs = xbits.abs().uintval();
+ // |x| >= 1
if (LIBC_UNLIKELY(x_abs >= 0x3c00U)) {
- if (xbits.is_nan()) {
+ if (xbits.is_nan())
return x;
- }
+
// |x| == 1.0
if (x_abs == 0x3c00U) {
fputil::set_errno_if_required(ERANGE);
fputil::raise_except_if_required(FE_DIVBYZERO);
return FPBits::inf(sign).get_val();
- } else {
- fputil::set_errno_if_required(EDOM);
- fputil::raise_except_if_required(FE_INVALID);
- return FPBits::quiet_nan().get_val();
}
+ // |x| > 1.0
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ return FPBits::quiet_nan().get_val();
}
// For |x| less than approximately 0.10
@@ -52,35 +54,25 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
// atanh(x) ≈ x + (1/3)*x^3
if (LIBC_UNLIKELY(x_abs < 0x0100U)) {
return static_cast<float16>(
- LIBC_UNLIKELY(x_abs == 0) ? x : (x + 0x1.555556p-2 * x * x * x));
+ LIBC_UNLIKELY(x_abs == 0) ? x : (x + 0x1.555556p-2f * x * x * x));
}
// For 0x0100U <= |x| <= 0x2e66U:
// Let t = x^2.
// Define P(t) ≈ (1/3)*t + (1/5)*t^2 + (1/7)*t^3 + (1/9)*t^4 + (1/11)*t^5.
- // The coefficients below were derived using Sollya:
- // > display = hexadecimal;
- // > round(1/3, SG, RN);
- // > round(1/5, SG, RN);
- // > round(1/7, SG, RN);
- // > round(1/9, SG, RN);
- // > round(1/11, SG, RN);
- // This yields:
- // 0x1.555556p-2
- // 0x1.99999ap-3
- // 0x1.24924ap-3
- // 0x1.c71c72p-4
- // 0x1.745d18p-4f
+ // Coefficients (from Sollya, RN, hexadecimal):
+ // 1/3 = 0x1.555556p-2, 1/5 = 0x1.99999ap-3, 1/7 = 0x1.24924ap-3,
+ // 1/9 = 0x1.c71c72p-4, 1/11 = 0x1.745d18p-4
// Thus, atanh(x) ≈ x * (1 + P(x^2)).
float xf = x;
float x2 = xf * xf;
float pe = fputil::polyeval(x2, 0.0f, 0x1.555556p-2f, 0x1.99999ap-3f,
0x1.24924ap-3f, 0x1.c71c72p-4f, 0x1.745d18p-4f);
- return static_cast<float16>(fputil::multiply_add(xf, pe, xf));
+ return fputil::cast<float16>(fputil::multiply_add(xf, pe, xf));
}
float xf = x;
- return static_cast<float16>(0.5 * log_eval((xf + 1.0) / (xf - 1.0)));
+ return fputil::cast<float16>(0.5 * log_eval((xf + 1.0) / (xf - 1.0)));
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/atanhf16_test.cpp b/libc/test/src/math/atanhf16_test.cpp
index 6637c7eead8f5..a534e05230e9d 100644
--- a/libc/test/src/math/atanhf16_test.cpp
+++ b/libc/test/src/math/atanhf16_test.cpp
@@ -7,7 +7,6 @@
//===----------------------------------------------------------------------===//
#include "src/__support/FPUtil/FPBits.h"
-#include "src/errno/libc_errno.h"
#include "src/math/atanhf16.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
@@ -17,13 +16,13 @@
using LlvmLibcAtanhf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
-// Range for positive numbers: [0, 1)
+// Range for positive numbers: [0, +Inf]
static constexpr uint16_t POS_START = 0x0000U;
-static constexpr uint16_t POS_STOP = 0x3C00;
+static constexpr uint16_t POS_STOP = 0x7C00U;
-// Range for negative numbers: (-1, 0]
-static constexpr uint16_t NEG_START = 0xBBFFU;
-static constexpr uint16_t NEG_STOP = 0x8000U;
+// Range for negative numbers: [-Inf, 0]
+static constexpr uint16_t NEG_START = 0x8000U;
+static constexpr uint16_t NEG_STOP = 0xFC00U;
TEST_F(LlvmLibcAtanhf16Test, PositiveRange) {
for (uint16_t v = POS_START; v < POS_STOP; ++v) {
@@ -34,7 +33,7 @@ TEST_F(LlvmLibcAtanhf16Test, PositiveRange) {
}
TEST_F(LlvmLibcAtanhf16Test, NegativeRange) {
- for (uint16_t v = NEG_START; v >= NEG_STOP; --v) {
+ for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
float16 x = FPBits(v).get_val();
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atanh, x,
LIBC_NAMESPACE::atanhf16(x), 0.5);
diff --git a/libc/test/src/math/smoke/atanhf16_test.cpp b/libc/test/src/math/smoke/atanhf16_test.cpp
index 1ac483d2ae758..1e81d5f18c0b2 100644
--- a/libc/test/src/math/smoke/atanhf16_test.cpp
+++ b/libc/test/src/math/smoke/atanhf16_test.cpp
@@ -22,9 +22,7 @@ TEST_F(LlvmLibcAtanhf16Test, SpecialNumbers) {
EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::atanhf16(zero));
EXPECT_MATH_ERRNO(0);
- EXPECT_FP_EQ_ALL_ROUNDING(
- -0.0f,
- LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast<float16>(-0.0f)));
+ EXPECT_FP_EQ_ALL_ROUNDING(neg_zero, LIBC_NAMESPACE::atanhf16(neg_zero));
EXPECT_MATH_ERRNO(0);
EXPECT_FP_EQ_WITH_EXCEPTION(
>From 9276a1515a0d84508de1fc13ba108a4eb56c8458 Mon Sep 17 00:00:00 2001
From: Harrison Hao <tsworld1314 at gmail.com>
Date: Sat, 29 Mar 2025 08:31:20 +0000
Subject: [PATCH 4/8] [libc][math] Update.
---
libc/src/math/generic/atanhf16.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp
index bb16a78a9862d..6b4d6abd15ba0 100644
--- a/libc/src/math/generic/atanhf16.cpp
+++ b/libc/src/math/generic/atanhf16.cpp
@@ -72,7 +72,7 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
}
float xf = x;
- return fputil::cast<float16>(0.5 * log_eval((xf + 1.0) / (xf - 1.0)));
+ return fputil::cast<float16>(0.5f * log_eval((xf + 1.0f) / (xf - 1.0f)));
}
} // namespace LIBC_NAMESPACE_DECL
>From 2a2c9ff6e3a442379fc5ae0c00b69a33549ff7d8 Mon Sep 17 00:00:00 2001
From: Harrison Hao <tsworld1314 at gmail.com>
Date: Sat, 29 Mar 2025 09:49:10 +0000
Subject: [PATCH 5/8] [libc][math] Add except results.
---
libc/src/math/generic/atanhf16.cpp | 14 ++++++++++++++
1 file changed, 14 insertions(+)
diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp
index 6b4d6abd15ba0..4569384a83349 100644
--- a/libc/src/math/generic/atanhf16.cpp
+++ b/libc/src/math/generic/atanhf16.cpp
@@ -14,6 +14,7 @@
#include "src/__support/FPUtil/FPBits.h"
#include "src/__support/FPUtil/PolyEval.h"
#include "src/__support/FPUtil/cast.h"
+#include "src/__support/FPUtil/except_value_utils.h"
#include "src/__support/FPUtil/multiply_add.h"
#include "src/__support/common.h"
#include "src/__support/macros/config.h"
@@ -21,6 +22,15 @@
namespace LIBC_NAMESPACE_DECL {
+static constexpr size_t N_EXCEPTS = 2;
+static constexpr fputil::ExceptValues<float16, N_EXCEPTS> ATANHF16_EXCEPTS{{
+ // (input, RZ output, RU offset, RD offset, RN offset)
+ // x = 0x1.a5cp-4, atanhf16(x) = 0x1.a74p-4 (RZ)
+ {0x2E97, 0x2E9D, 1, 0, 0},
+ // x = -0x1.a5cp-4, atanhf16(x) = -0x1.a74p-4 (RZ)
+ {0xAE97, 0xAE9D, 0, 1, 0},
+}};
+
LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
using FPBits = fputil::FPBits<float16>;
@@ -45,6 +55,10 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
return FPBits::quiet_nan().get_val();
}
+ if (auto r = ATANHF16_EXCEPTS.lookup(xbits.uintval());
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+
// For |x| less than approximately 0.10
if (LIBC_UNLIKELY(x_abs <= 0x2e66U)) {
// The Taylor expansion of atanh(x) is:
>From c092cad2d31082dc04ce020f1a5ec11a88393e6d Mon Sep 17 00:00:00 2001
From: Harrison Hao <tsworld1314 at gmail.com>
Date: Mon, 31 Mar 2025 23:00:57 +0800
Subject: [PATCH 6/8] [libc][math] Update for comments.
---
libc/src/math/generic/CMakeLists.txt | 1 +
libc/src/math/generic/atanhf16.cpp | 13 ++++++++-----
libc/test/src/math/atanhf16_test.cpp | 1 -
libc/test/src/math/smoke/CMakeLists.txt | 2 +-
4 files changed, 10 insertions(+), 7 deletions(-)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index a83107f248562..b5672e9d4b60f 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -4040,6 +4040,7 @@ add_entrypoint_object(
libc.hdr.errno_macros
libc.hdr.fenv_macros
libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.except_value_utils
libc.src.__support.FPUtil.fenv_impl
libc.src.__support.FPUtil.fp_bits
libc.src.__support.FPUtil.multiply_add
diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp
index 4569384a83349..a951112727eec 100644
--- a/libc/src/math/generic/atanhf16.cpp
+++ b/libc/src/math/generic/atanhf16.cpp
@@ -61,17 +61,20 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
// For |x| less than approximately 0.10
if (LIBC_UNLIKELY(x_abs <= 0x2e66U)) {
+ // atanh(+/-0) = +/-0
+ if (LIBC_UNLIKELY(x_abs == 0U))
+ return x;
// The Taylor expansion of atanh(x) is:
// atanh(x) = x + x^3/3 + x^5/5 + x^7/7 + x^9/9 + x^11/11
// = x * [1 + x^2/3 + x^4/5 + x^6/7 + x^8/9 + x^10/11]
- // When |x| < 0x0100U, this can be approximated by:
+ // When |x| < 2^-16, this can be approximated by:
// atanh(x) ≈ x + (1/3)*x^3
if (LIBC_UNLIKELY(x_abs < 0x0100U)) {
- return static_cast<float16>(
- LIBC_UNLIKELY(x_abs == 0) ? x : (x + 0x1.555556p-2f * x * x * x));
+ float xf = x;
+ return fputil::cast<float16>(xf + 0x1.555556p-2f * xf * xf * xf);
}
- // For 0x0100U <= |x| <= 0x2e66U:
+ // For 2^-16 <= |x| <= 0x1.998p-4 (~0.10):
// Let t = x^2.
// Define P(t) ≈ (1/3)*t + (1/5)*t^2 + (1/7)*t^3 + (1/9)*t^4 + (1/11)*t^5.
// Coefficients (from Sollya, RN, hexadecimal):
@@ -86,7 +89,7 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
}
float xf = x;
- return fputil::cast<float16>(0.5f * log_eval((xf + 1.0f) / (xf - 1.0f)));
+ return fputil::cast<float16>(0.5 * log_eval((xf + 1.0f) / (xf - 1.0f)));
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/test/src/math/atanhf16_test.cpp b/libc/test/src/math/atanhf16_test.cpp
index a534e05230e9d..9f0a493f1bb0e 100644
--- a/libc/test/src/math/atanhf16_test.cpp
+++ b/libc/test/src/math/atanhf16_test.cpp
@@ -6,7 +6,6 @@
//
//===----------------------------------------------------------------------===//
-#include "src/__support/FPUtil/FPBits.h"
#include "src/math/atanhf16.h"
#include "test/UnitTest/FPMatcher.h"
#include "test/UnitTest/Test.h"
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index ec867339152e3..c989539b30ed9 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3956,7 +3956,7 @@ add_fp_unittest(
DEPENDS
libc.src.errno.errno
libc.src.math.atanhf16
- libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.cast
)
add_fp_unittest(
>From 81c5cf7543fbf0d3789a48ad34ff33e93116ede9 Mon Sep 17 00:00:00 2001
From: Harrison Hao <tsworld1314 at gmail.com>
Date: Mon, 31 Mar 2025 23:48:05 +0800
Subject: [PATCH 7/8] [libc][math] Update for comments again.
---
libc/test/src/math/atanhf16_test.cpp | 2 +-
1 file changed, 1 insertion(+), 1 deletion(-)
diff --git a/libc/test/src/math/atanhf16_test.cpp b/libc/test/src/math/atanhf16_test.cpp
index 9f0a493f1bb0e..e35cc775b0609 100644
--- a/libc/test/src/math/atanhf16_test.cpp
+++ b/libc/test/src/math/atanhf16_test.cpp
@@ -24,7 +24,7 @@ static constexpr uint16_t NEG_START = 0x8000U;
static constexpr uint16_t NEG_STOP = 0xFC00U;
TEST_F(LlvmLibcAtanhf16Test, PositiveRange) {
- for (uint16_t v = POS_START; v < POS_STOP; ++v) {
+ for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
float16 x = FPBits(v).get_val();
EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Atanh, x,
LIBC_NAMESPACE::atanhf16(x), 0.5);
>From 94ff9a7ceeed3bb5d050a8dadb25fab5779ec8fc Mon Sep 17 00:00:00 2001
From: Harrison Hao <tsworld1314 at gmail.com>
Date: Sat, 12 Apr 2025 18:51:33 +0800
Subject: [PATCH 8/8] [libc] Add float log_eval.
---
libc/src/math/generic/atanhf16.cpp | 8 ++-
libc/src/math/generic/common_constants.cpp | 72 ++++++++++++++++++++++
libc/src/math/generic/common_constants.h | 8 +++
libc/src/math/generic/explogxf.h | 34 ++++++++++
4 files changed, 120 insertions(+), 2 deletions(-)
diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp
index a951112727eec..c46fa961b432a 100644
--- a/libc/src/math/generic/atanhf16.cpp
+++ b/libc/src/math/generic/atanhf16.cpp
@@ -22,13 +22,17 @@
namespace LIBC_NAMESPACE_DECL {
-static constexpr size_t N_EXCEPTS = 2;
+static constexpr size_t N_EXCEPTS = 4;
static constexpr fputil::ExceptValues<float16, N_EXCEPTS> ATANHF16_EXCEPTS{{
// (input, RZ output, RU offset, RD offset, RN offset)
// x = 0x1.a5cp-4, atanhf16(x) = 0x1.a74p-4 (RZ)
{0x2E97, 0x2E9D, 1, 0, 0},
// x = -0x1.a5cp-4, atanhf16(x) = -0x1.a74p-4 (RZ)
{0xAE97, 0xAE9D, 0, 1, 0},
+ // x = -0x1.99cp-4, atanhf16(x) = -0x1.9bp-4 (RZ)
+ {0xAE67, 0xAE6C, 0, 1, 1},
+ // x = -0x1.b8cp-3, atanhf16(x) = -0x1.bfcp-3 (RZ)
+ {0xB2E3, 0xB2FF, 0, 1, 0},
}};
LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
@@ -89,7 +93,7 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) {
}
float xf = x;
- return fputil::cast<float16>(0.5 * log_eval((xf + 1.0f) / (xf - 1.0f)));
+ return fputil::cast<float16>(0.5 * log_eval_f((xf + 1.0f) / (xf - 1.0f)));
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/common_constants.cpp b/libc/src/math/generic/common_constants.cpp
index 3088ef96e3b93..7f8a1244488c7 100644
--- a/libc/src/math/generic/common_constants.cpp
+++ b/libc/src/math/generic/common_constants.cpp
@@ -13,6 +13,42 @@
namespace LIBC_NAMESPACE_DECL {
+// Lookup table for logf(f) = logf(1 + n*2^(-7)) where n = 0..127,
+// computed and stored as float precision constants.
+const float ONE_OVER_F_FLOAT[128] = {
+ 0x1p0f, 0x1.fc07fp-1f, 0x1.f81f82p-1f, 0x1.f4465ap-1f,
+ 0x1.f07c2p-1f, 0x1.ecc07cp-1f, 0x1.e9131ap-1f, 0x1.e573acp-1f,
+ 0x1.e1e1e2p-1f, 0x1.de5d6ep-1f, 0x1.dae608p-1f, 0x1.d77b66p-1f,
+ 0x1.d41d42p-1f, 0x1.d0cb58p-1f, 0x1.cd8568p-1f, 0x1.ca4b3p-1f,
+ 0x1.c71c72p-1f, 0x1.c3f8fp-1f, 0x1.c0e07p-1f, 0x1.bdd2b8p-1f,
+ 0x1.bacf92p-1f, 0x1.b7d6c4p-1f, 0x1.b4e81cp-1f, 0x1.b20364p-1f,
+ 0x1.af286cp-1f, 0x1.ac5702p-1f, 0x1.a98ef6p-1f, 0x1.a6d01ap-1f,
+ 0x1.a41a42p-1f, 0x1.a16d4p-1f, 0x1.9ec8eap-1f, 0x1.9c2d14p-1f,
+ 0x1.99999ap-1f, 0x1.970e5p-1f, 0x1.948b1p-1f, 0x1.920fb4p-1f,
+ 0x1.8f9c18p-1f, 0x1.8d3018p-1f, 0x1.8acb9p-1f, 0x1.886e6p-1f,
+ 0x1.861862p-1f, 0x1.83c978p-1f, 0x1.818182p-1f, 0x1.7f406p-1f,
+ 0x1.7d05f4p-1f, 0x1.7ad22p-1f, 0x1.78a4c8p-1f, 0x1.767dcep-1f,
+ 0x1.745d18p-1f, 0x1.724288p-1f, 0x1.702e06p-1f, 0x1.6e1f76p-1f,
+ 0x1.6c16c2p-1f, 0x1.6a13cep-1f, 0x1.681682p-1f, 0x1.661ec6p-1f,
+ 0x1.642c86p-1f, 0x1.623fa8p-1f, 0x1.605816p-1f, 0x1.5e75bcp-1f,
+ 0x1.5c9882p-1f, 0x1.5ac056p-1f, 0x1.58ed24p-1f, 0x1.571ed4p-1f,
+ 0x1.555556p-1f, 0x1.539094p-1f, 0x1.51d07ep-1f, 0x1.501502p-1f,
+ 0x1.4e5e0ap-1f, 0x1.4cab88p-1f, 0x1.4afd6ap-1f, 0x1.49539ep-1f,
+ 0x1.47ae14p-1f, 0x1.460cbcp-1f, 0x1.446f86p-1f, 0x1.42d662p-1f,
+ 0x1.414142p-1f, 0x1.3fb014p-1f, 0x1.3e22ccp-1f, 0x1.3c995ap-1f,
+ 0x1.3b13b2p-1f, 0x1.3991c2p-1f, 0x1.381382p-1f, 0x1.3698ep-1f,
+ 0x1.3521dp-1f, 0x1.33ae46p-1f, 0x1.323e34p-1f, 0x1.30d19p-1f,
+ 0x1.2f684cp-1f, 0x1.2e025cp-1f, 0x1.2c9fb4p-1f, 0x1.2b404ap-1f,
+ 0x1.29e412p-1f, 0x1.288b02p-1f, 0x1.27350cp-1f, 0x1.25e228p-1f,
+ 0x1.24924ap-1f, 0x1.234568p-1f, 0x1.21fb78p-1f, 0x1.20b47p-1f,
+ 0x1.1f7048p-1f, 0x1.1e2ef4p-1f, 0x1.1cf06ap-1f, 0x1.1bb4a4p-1f,
+ 0x1.1a7b96p-1f, 0x1.194538p-1f, 0x1.181182p-1f, 0x1.16e068p-1f,
+ 0x1.15b1e6p-1f, 0x1.1485fp-1f, 0x1.135c82p-1f, 0x1.12358ep-1f,
+ 0x1.111112p-1f, 0x1.0fef02p-1f, 0x1.0ecf56p-1f, 0x1.0db20ap-1f,
+ 0x1.0c9714p-1f, 0x1.0b7e6ep-1f, 0x1.0a681p-1f, 0x1.0953f4p-1f,
+ 0x1.08421p-1f, 0x1.07326p-1f, 0x1.0624dep-1f, 0x1.05198p-1f,
+ 0x1.041042p-1f, 0x1.03091cp-1f, 0x1.020408p-1f, 0x1.010102p-1f};
+
// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
const double ONE_OVER_F[128] = {
0x1.0000000000000p+0, 0x1.fc07f01fc07f0p-1, 0x1.f81f81f81f820p-1,
@@ -59,6 +95,42 @@ const double ONE_OVER_F[128] = {
0x1.05197f7d73404p-1, 0x1.0410410410410p-1, 0x1.03091b51f5e1ap-1,
0x1.0204081020408p-1, 0x1.0101010101010p-1};
+// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127,
+// computed and stored as float precision constants.
+const float LOG_F_FLOAT[128] = {
+ 0.0f, 0x1.fe02a6p-8f, 0x1.fc0a8cp-7f, 0x1.7b91bp-6f,
+ 0x1.f829bp-6f, 0x1.39e87cp-5f, 0x1.77459p-5f, 0x1.b42dd8p-5f,
+ 0x1.f0a30cp-5f, 0x1.16536ep-4f, 0x1.341d7ap-4f, 0x1.51b074p-4f,
+ 0x1.6f0d28p-4f, 0x1.8c345ep-4f, 0x1.a926d4p-4f, 0x1.c5e548p-4f,
+ 0x1.e27076p-4f, 0x1.fec914p-4f, 0x1.0d77e8p-3f, 0x1.1b72aep-3f,
+ 0x1.29553p-3f, 0x1.371fc2p-3f, 0x1.44d2b6p-3f, 0x1.526e5ep-3f,
+ 0x1.5ff308p-3f, 0x1.6d60fep-3f, 0x1.7ab89p-3f, 0x1.87fa06p-3f,
+ 0x1.9525aap-3f, 0x1.a23bc2p-3f, 0x1.af3c94p-3f, 0x1.bc2868p-3f,
+ 0x1.c8ff7cp-3f, 0x1.d5c216p-3f, 0x1.e27076p-3f, 0x1.ef0adcp-3f,
+ 0x1.fb9186p-3f, 0x1.04025ap-2f, 0x1.0a324ep-2f, 0x1.1058cp-2f,
+ 0x1.1675cap-2f, 0x1.1c898cp-2f, 0x1.22942p-2f, 0x1.2895a2p-2f,
+ 0x1.2e8e2cp-2f, 0x1.347ddap-2f, 0x1.3a64c6p-2f, 0x1.404308p-2f,
+ 0x1.4618bcp-2f, 0x1.4be5fap-2f, 0x1.51aad8p-2f, 0x1.576772p-2f,
+ 0x1.5d1bdcp-2f, 0x1.62c83p-2f, 0x1.686c82p-2f, 0x1.6e08eap-2f,
+ 0x1.739d8p-2f, 0x1.792a56p-2f, 0x1.7eaf84p-2f, 0x1.842d1ep-2f,
+ 0x1.89a338p-2f, 0x1.8f11e8p-2f, 0x1.947942p-2f, 0x1.99d958p-2f,
+ 0x1.9f323ep-2f, 0x1.a4840ap-2f, 0x1.a9cecap-2f, 0x1.af1294p-2f,
+ 0x1.b44f78p-2f, 0x1.b9858ap-2f, 0x1.beb4dap-2f, 0x1.c3dd7ap-2f,
+ 0x1.c8ff7cp-2f, 0x1.ce1afp-2f, 0x1.d32fe8p-2f, 0x1.d83e72p-2f,
+ 0x1.dd46ap-2f, 0x1.e24882p-2f, 0x1.e74426p-2f, 0x1.ec399ep-2f,
+ 0x1.f128f6p-2f, 0x1.f6124p-2f, 0x1.faf588p-2f, 0x1.ffd2ep-2f,
+ 0x1.02552ap-1f, 0x1.04bdfap-1f, 0x1.0723e6p-1f, 0x1.0986f4p-1f,
+ 0x1.0be72ep-1f, 0x1.0e4498p-1f, 0x1.109f3ap-1f, 0x1.12f71ap-1f,
+ 0x1.154c3ep-1f, 0x1.179eacp-1f, 0x1.19ee6cp-1f, 0x1.1c3b82p-1f,
+ 0x1.1e85f6p-1f, 0x1.20cdcep-1f, 0x1.23130ep-1f, 0x1.2555bcp-1f,
+ 0x1.2795e2p-1f, 0x1.29d38p-1f, 0x1.2c0e9ep-1f, 0x1.2e4744p-1f,
+ 0x1.307d74p-1f, 0x1.32b134p-1f, 0x1.34e28ap-1f, 0x1.37117cp-1f,
+ 0x1.393e0ep-1f, 0x1.3b6844p-1f, 0x1.3d9026p-1f, 0x1.3fb5b8p-1f,
+ 0x1.41d8fep-1f, 0x1.43f9fep-1f, 0x1.4618bcp-1f, 0x1.48353ep-1f,
+ 0x1.4a4f86p-1f, 0x1.4c679ap-1f, 0x1.4e7d82p-1f, 0x1.50913cp-1f,
+ 0x1.52a2d2p-1f, 0x1.54b246p-1f, 0x1.56bf9ep-1f, 0x1.58cadcp-1f,
+ 0x1.5ad404p-1f, 0x1.5cdb1ep-1f, 0x1.5ee02ap-1f, 0x1.60e33p-1f};
+
// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127.
const double LOG_F[128] = {
0x0.0000000000000p+0, 0x1.fe02a6b106788p-8, 0x1.fc0a8b0fc03e3p-7,
diff --git a/libc/src/math/generic/common_constants.h b/libc/src/math/generic/common_constants.h
index dc1a90c1cb24a..e65f002845953 100644
--- a/libc/src/math/generic/common_constants.h
+++ b/libc/src/math/generic/common_constants.h
@@ -15,9 +15,17 @@
namespace LIBC_NAMESPACE_DECL {
+// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127,
+// computed and stored as float precision constants.
+extern const float ONE_OVER_F_FLOAT[128];
+
// Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127.
extern const double ONE_OVER_F[128];
+// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127,
+// computed and stored as float precision constants.
+extern const float LOG_F_FLOAT[128];
+
// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127.
extern const double LOG_F[128];
diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h
index e79aa13eb57f7..91285c1631111 100644
--- a/libc/src/math/generic/explogxf.h
+++ b/libc/src/math/generic/explogxf.h
@@ -297,6 +297,40 @@ LIBC_INLINE static double log2_eval(double x) {
return result;
}
+// x should be positive, normal finite value
+LIBC_INLINE static float log_eval_f(float x) {
+ // For x = 2^ex * (1 + mx), logf(x) = ex * logf(2) + logf(1 + mx).
+ using FPB = fputil::FPBits<float>;
+ FPB bs(x);
+
+ float ex = static_cast<float>(bs.get_exponent());
+ // p1 is the leading 7 bits of mx, i.e.
+ // p1 * 2^(-7) <= m_x < (p1 + 1) * 2^(-7).
+ int p1 = static_cast<int>(bs.get_mantissa() >> (FPB::FRACTION_LEN - 7));
+
+ // Set bs to (1 + (mx - p1*2^(-7))
+ bs.set_uintval(bs.uintval() & (FPB::FRACTION_MASK >> 7));
+ bs.set_biased_exponent(FPB::EXP_BIAS);
+ // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)).
+ float dx = (bs.get_val() - 1.0f) * ONE_OVER_F_FLOAT[p1];
+
+ // Minimax polynomial of log(1 + dx) generated by Sollya with:
+ // > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]);
+ const float COEFFS[6] = {-0x1.fffffep-2f, 0x1.555556p-2f, -0x1.fffefep-3f,
+ 0x1.99999ap-3f, -0x1.554318p-3f, 0x1.1dc5c4p-3f};
+
+ float dx2 = dx * dx;
+
+ float c1 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]);
+ float c2 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]);
+ float c3 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]);
+
+ float p = fputil::polyeval(dx2, dx, c1, c2, c3);
+
+ float result = fputil::multiply_add(ex, 0x1.62e42ep-1f, LOG_F_FLOAT[p1] + p);
+ return result;
+}
+
// x should be positive, normal finite value
LIBC_INLINE static double log_eval(double x) {
// For x = 2^ex * (1 + mx)
More information about the libc-commits
mailing list