From 41e7e797ce6bdefd9b561abeaa716c69ae27f252 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Mon, 17 Mar 2025 15:34:24 +0000 Subject: [PATCH 01/17] [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 73dfeae1a2c94..33d16b6c01ac2 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -664,6 +664,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; + + 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( + 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(fputil::multiply_add(xf, pe, xf)); + } + + float xf = x; + return static_cast(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 + +using LlvmLibcAtanhf16Test = LIBC_NAMESPACE::testing::FPTest; +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; + +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(-0.0f))); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ_WITH_EXCEPTION( + inf, + LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast(1.0f)), + FE_DIVBYZERO); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ_WITH_EXCEPTION( + neg_inf, + LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast(-1.0f)), + FE_DIVBYZERO); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_IS_NAN_WITH_EXCEPTION( + LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast(2.0f)), + FE_INVALID); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_IS_NAN_WITH_EXCEPTION( + LIBC_NAMESPACE::atanhf16(LIBC_NAMESPACE::fputil::cast(-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 89c8cc08d93565c42ec8860014051e955500350f Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Wed, 26 Mar 2025 10:18:23 +0800 Subject: [PATCH 02/17] [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; 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 988f79b83f12ba0c4448a86985a8fd692254d86f Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Sat, 29 Mar 2025 06:10:45 +0000 Subject: [PATCH 03/17] [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; + using FPBits = fputil::FPBits; 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( - 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(fputil::multiply_add(xf, pe, xf)); + return fputil::cast(fputil::multiply_add(xf, pe, xf)); } float xf = x; - return static_cast(0.5 * log_eval((xf + 1.0) / (xf - 1.0))); + return fputil::cast(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; 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(-0.0f))); + EXPECT_FP_EQ_ALL_ROUNDING(neg_zero, LIBC_NAMESPACE::atanhf16(neg_zero)); EXPECT_MATH_ERRNO(0); EXPECT_FP_EQ_WITH_EXCEPTION( From ac3bf294c9f7bb0a57db411082f54ce831307a6f Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Sat, 29 Mar 2025 08:31:20 +0000 Subject: [PATCH 04/17] [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(0.5 * log_eval((xf + 1.0) / (xf - 1.0))); + return fputil::cast(0.5f * log_eval((xf + 1.0f) / (xf - 1.0f))); } } // namespace LIBC_NAMESPACE_DECL From 2a5dd72fc6bf91e1f62fa79bb5450ca0d3d4fe9c Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Sat, 29 Mar 2025 09:49:10 +0000 Subject: [PATCH 05/17] [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 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; @@ -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 5f995ee0605efb268995949c6756507760d3fc50 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Mon, 31 Mar 2025 23:00:57 +0800 Subject: [PATCH 06/17] [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( - LIBC_UNLIKELY(x_abs == 0) ? x : (x + 0x1.555556p-2f * x * x * x)); + float xf = x; + return fputil::cast(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(0.5f * log_eval((xf + 1.0f) / (xf - 1.0f))); + return fputil::cast(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 7057a6e41102a966a4ec4d762118b9fac72af080 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Mon, 31 Mar 2025 23:48:05 +0800 Subject: [PATCH 07/17] [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 7d3a74e2d645d5ee6e1be554f9b712f3b16e061b Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Sat, 12 Apr 2025 18:51:33 +0800 Subject: [PATCH 08/17] [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 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(0.5 * log_eval((xf + 1.0f) / (xf - 1.0f))); + return fputil::cast(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; + FPB bs(x); + + float ex = static_cast(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(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) From da773e0525d5117feafd3a187b6145cc1f9be1ec Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 15 Apr 2025 15:53:18 +0000 Subject: [PATCH 09/17] [libc] Add a new version log_eval. --- libc/src/math/generic/atanhf16.cpp | 9 +- libc/src/math/generic/common_constants.cpp | 6 ++ libc/src/math/generic/explogxf.h | 106 ++++++++++----------- libc/test/src/math/smoke/atanhf16_test.cpp | 4 + 4 files changed, 68 insertions(+), 57 deletions(-) diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp index c46fa961b432a..406a8144bc62c 100644 --- a/libc/src/math/generic/atanhf16.cpp +++ b/libc/src/math/generic/atanhf16.cpp @@ -44,8 +44,13 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) { // |x| >= 1 if (LIBC_UNLIKELY(x_abs >= 0x3c00U)) { - if (xbits.is_nan()) + if (xbits.is_nan()) { + if (xbits.is_signaling_nan()) { + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } return x; + } // |x| == 1.0 if (x_abs == 0x3c00U) { @@ -93,7 +98,7 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) { } float xf = x; - return fputil::cast(0.5 * log_eval_f((xf + 1.0f) / (xf - 1.0f))); + return fputil::cast(0.5 * log_eval((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 7f8a1244488c7..a53641cf37bdb 100644 --- a/libc/src/math/generic/common_constants.cpp +++ b/libc/src/math/generic/common_constants.cpp @@ -15,6 +15,9 @@ 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. +// Generated by Sollya with the following commands: +// display = hexadecimal; +// for n from 0 to 127 do { print(single(1 / (1 + n / 128.0))); }; 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, @@ -97,6 +100,9 @@ const double ONE_OVER_F[128] = { // Lookup table for (1/f) where f = 1 + n*2^(-7), n = 0..127, // computed and stored as float precision constants. +// Generated by Sollya with the following commands: +// display = hexadecimal; +// for n from 0 to 127 do { print(single(log(1 + n / 128.0))); }; 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, diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h index 91285c1631111..7e68994a1bd82 100644 --- a/libc/src/math/generic/explogxf.h +++ b/libc/src/math/generic/explogxf.h @@ -133,6 +133,29 @@ struct exp_b_reduc_t { double lo; }; +// Coefficients for double (6th-degree minimax polynomial on [0, 2^-7]). +// Minimax polynomial of log(1 + dx) generated by Sollya with: +// > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]); +static constexpr double LOG_COEFFS_DOUBLE[6] = { + -0x1.ffffffffffffcp-2, 0x1.5555555552ddep-2, -0x1.ffffffefe562dp-3, + 0x1.9999817d3a50fp-3, -0x1.554317b3f67a5p-3, 0x1.1dc5c45e09c18p-3}; + +// Coefficients for float (6th-degree minimax polynomial on [0, 2^-7]). +// Minimax polynomial of log(1 + dx) generated by Sollya with: +// > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]); +static constexpr float LOG_COEFFS_FLOAT[6] = {-0x1.fffffep-2f, 0x1.555556p-2f, + -0x1.fffefep-3f, 0x1.99999ap-3f, + -0x1.554318p-3f, 0x1.1dc5c4p-3f}; + +// log(2) in double precision. +static constexpr double LOG2_DOUBLE = 0x1.62e42fefa39efp-1; + +// log(2) in float precision. +// Generated by Sollya with the following commands: +// > display = hexadecimal; +// > round(log(2), SG, RN); +static constexpr float LOG2_FLOAT = 0x1.62e43p-1f; + // The function correctly calculates b^x value with at least float precision // in a limited range. // Range reduction: @@ -297,73 +320,46 @@ 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) { +template LIBC_INLINE static T log_eval(T x) { // For x = 2^ex * (1 + mx), logf(x) = ex * logf(2) + logf(1 + mx). - using FPB = fputil::FPBits; - FPB bs(x); + using FPBits = fputil::FPBits; + FPBits xbits(x); - float ex = static_cast(bs.get_exponent()); + T ex = static_cast(xbits.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(bs.get_mantissa() >> (FPB::FRACTION_LEN - 7)); + int p1 = static_cast(xbits.get_mantissa() >> (FPBits::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); + xbits.set_uintval(xbits.uintval() & (FPBits::FRACTION_MASK >> 7)); + xbits.set_biased_exponent(FPBits::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]); + T dx = static_cast(xbits.get_val() - 1.0) * + (std::is_same::value ? static_cast(ONE_OVER_F[p1]) + : ONE_OVER_F_FLOAT[p1]); - float p = fputil::polyeval(dx2, dx, c1, c2, c3); + T dx2 = dx * dx; - 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) - // log(x) = ex * log(2) + log(1 + mx) - using FPB = fputil::FPBits; - FPB bs(x); - - double ex = static_cast(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(bs.get_mantissa() >> (FPB::FRACTION_LEN - 7)); + const T *coeffs = nullptr; + T log2_val; + if constexpr (std::is_same::value) { + coeffs = LOG_COEFFS_DOUBLE; + log2_val = LOG2_DOUBLE; + } else { + coeffs = LOG_COEFFS_FLOAT; + log2_val = LOG2_FLOAT; + } - // 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)). - double dx = (bs.get_val() - 1.0) * ONE_OVER_F[p1]; + T c1 = fputil::multiply_add(dx, coeffs[1], coeffs[0]); + T c2 = fputil::multiply_add(dx, coeffs[3], coeffs[2]); + T c3 = fputil::multiply_add(dx, coeffs[5], coeffs[4]); - // Minimax polynomial of log(1 + dx) generated by Sollya with: - // > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]); - const double COEFFS[6] = {-0x1.ffffffffffffcp-2, 0x1.5555555552ddep-2, - -0x1.ffffffefe562dp-3, 0x1.9999817d3a50fp-3, - -0x1.554317b3f67a5p-3, 0x1.1dc5c45e09c18p-3}; - double dx2 = dx * dx; - double c1 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]); - double c2 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]); - double c3 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]); + T p = fputil::polyeval(dx2, dx, c1, c2, c3); - double p = fputil::polyeval(dx2, dx, c1, c2, c3); - double result = - fputil::multiply_add(ex, /*log(2)*/ 0x1.62e42fefa39efp-1, LOG_F[p1] + p); - return result; + if constexpr (std::is_same::value) + return fputil::multiply_add(ex, log2_val, LOG_F[p1] + p); + else + return fputil::multiply_add(ex, log2_val, LOG_F_FLOAT[p1] + p); } // Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We diff --git a/libc/test/src/math/smoke/atanhf16_test.cpp b/libc/test/src/math/smoke/atanhf16_test.cpp index 1e81d5f18c0b2..81df6da8cee26 100644 --- a/libc/test/src/math/smoke/atanhf16_test.cpp +++ b/libc/test/src/math/smoke/atanhf16_test.cpp @@ -16,6 +16,10 @@ using LlvmLibcAtanhf16Test = LIBC_NAMESPACE::testing::FPTest; TEST_F(LlvmLibcAtanhf16Test, SpecialNumbers) { LIBC_NAMESPACE::libc_errno = 0; + EXPECT_FP_EQ_WITH_EXCEPTION_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::atanhf16(sNaN), + FE_INVALID); + EXPECT_MATH_ERRNO(0); + EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::atanhf16(aNaN)); EXPECT_MATH_ERRNO(0); From 7dd2bd24cac3bece725f70674d19760be7131a20 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 15 Apr 2025 15:54:42 +0000 Subject: [PATCH 10/17] [libc] Update. --- libc/src/math/generic/explogxf.h | 1 + 1 file changed, 1 insertion(+) diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h index 7e68994a1bd82..fbd0e4119c6bb 100644 --- a/libc/src/math/generic/explogxf.h +++ b/libc/src/math/generic/explogxf.h @@ -320,6 +320,7 @@ LIBC_INLINE static double log2_eval(double x) { return result; } +// x should be positive, normal finite value template LIBC_INLINE static T log_eval(T x) { // For x = 2^ex * (1 + mx), logf(x) = ex * logf(2) + logf(1 + mx). using FPBits = fputil::FPBits; From 735a6672108a3ada96e3e6a4c70afbfc6b94b786 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 15 Apr 2025 16:08:00 +0000 Subject: [PATCH 11/17] [libc] Update code. --- libc/src/math/generic/atanhf16.cpp | 2 +- libc/src/math/generic/explogxf.h | 102 ++++++++++++++++------------- 2 files changed, 56 insertions(+), 48 deletions(-) diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp index 406a8144bc62c..7e73b69d3a66b 100644 --- a/libc/src/math/generic/atanhf16.cpp +++ b/libc/src/math/generic/atanhf16.cpp @@ -98,7 +98,7 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) { } float xf = x; - return fputil::cast(0.5 * log_eval((xf + 1.0f) / (xf - 1.0f))); + return fputil::cast(0.5 * log_eval_f((xf + 1.0f) / (xf - 1.0f))); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h index fbd0e4119c6bb..f02806fd2b3f8 100644 --- a/libc/src/math/generic/explogxf.h +++ b/libc/src/math/generic/explogxf.h @@ -133,29 +133,6 @@ struct exp_b_reduc_t { double lo; }; -// Coefficients for double (6th-degree minimax polynomial on [0, 2^-7]). -// Minimax polynomial of log(1 + dx) generated by Sollya with: -// > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]); -static constexpr double LOG_COEFFS_DOUBLE[6] = { - -0x1.ffffffffffffcp-2, 0x1.5555555552ddep-2, -0x1.ffffffefe562dp-3, - 0x1.9999817d3a50fp-3, -0x1.554317b3f67a5p-3, 0x1.1dc5c45e09c18p-3}; - -// Coefficients for float (6th-degree minimax polynomial on [0, 2^-7]). -// Minimax polynomial of log(1 + dx) generated by Sollya with: -// > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]); -static constexpr float LOG_COEFFS_FLOAT[6] = {-0x1.fffffep-2f, 0x1.555556p-2f, - -0x1.fffefep-3f, 0x1.99999ap-3f, - -0x1.554318p-3f, 0x1.1dc5c4p-3f}; - -// log(2) in double precision. -static constexpr double LOG2_DOUBLE = 0x1.62e42fefa39efp-1; - -// log(2) in float precision. -// Generated by Sollya with the following commands: -// > display = hexadecimal; -// > round(log(2), SG, RN); -static constexpr float LOG2_FLOAT = 0x1.62e43p-1f; - // The function correctly calculates b^x value with at least float precision // in a limited range. // Range reduction: @@ -321,12 +298,12 @@ LIBC_INLINE static double log2_eval(double x) { } // x should be positive, normal finite value -template LIBC_INLINE static T log_eval(T x) { +LIBC_INLINE static float log_eval_f(float x) { // For x = 2^ex * (1 + mx), logf(x) = ex * logf(2) + logf(1 + mx). - using FPBits = fputil::FPBits; + using FPBits = fputil::FPBits; FPBits xbits(x); - T ex = static_cast(xbits.get_exponent()); + float ex = static_cast(xbits.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(xbits.get_mantissa() >> (FPBits::FRACTION_LEN - 7)); @@ -335,32 +312,63 @@ template LIBC_INLINE static T log_eval(T x) { xbits.set_uintval(xbits.uintval() & (FPBits::FRACTION_MASK >> 7)); xbits.set_biased_exponent(FPBits::EXP_BIAS); // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)). - T dx = static_cast(xbits.get_val() - 1.0) * - (std::is_same::value ? static_cast(ONE_OVER_F[p1]) - : ONE_OVER_F_FLOAT[p1]); + float dx = (xbits.get_val() - 1.0f) * ONE_OVER_F_FLOAT[p1]; - T dx2 = dx * dx; + // 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}; - const T *coeffs = nullptr; - T log2_val; - if constexpr (std::is_same::value) { - coeffs = LOG_COEFFS_DOUBLE; - log2_val = LOG2_DOUBLE; - } else { - coeffs = LOG_COEFFS_FLOAT; - log2_val = LOG2_FLOAT; - } + float dx2 = dx * dx; - T c1 = fputil::multiply_add(dx, coeffs[1], coeffs[0]); - T c2 = fputil::multiply_add(dx, coeffs[3], coeffs[2]); - T c3 = fputil::multiply_add(dx, coeffs[5], coeffs[4]); + 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]); - T p = fputil::polyeval(dx2, dx, c1, c2, c3); + float p = fputil::polyeval(dx2, dx, c1, c2, c3); - if constexpr (std::is_same::value) - return fputil::multiply_add(ex, log2_val, LOG_F[p1] + p); - else - return fputil::multiply_add(ex, log2_val, LOG_F_FLOAT[p1] + p); + // Generated by Sollya with the following commands: + // > display = hexadecimal; + // > round(log(2), SG, RN); + static constexpr float LOGF_2 = 0x1.62e43p-1f; + + float result = fputil::multiply_add(ex, LOGF_2, 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) + // log(x) = ex * log(2) + log(1 + mx) + using FPB = fputil::FPBits; + FPB bs(x); + + double ex = static_cast(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(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)). + double dx = (bs.get_val() - 1.0) * ONE_OVER_F[p1]; + + // Minimax polynomial of log(1 + dx) generated by Sollya with: + // > P = fpminimax(log(1 + x)/x, 6, [|D...|], [0, 2^-7]); + const double COEFFS[6] = {-0x1.ffffffffffffcp-2, 0x1.5555555552ddep-2, + -0x1.ffffffefe562dp-3, 0x1.9999817d3a50fp-3, + -0x1.554317b3f67a5p-3, 0x1.1dc5c45e09c18p-3}; + double dx2 = dx * dx; + double c1 = fputil::multiply_add(dx, COEFFS[1], COEFFS[0]); + double c2 = fputil::multiply_add(dx, COEFFS[3], COEFFS[2]); + double c3 = fputil::multiply_add(dx, COEFFS[5], COEFFS[4]); + + double p = fputil::polyeval(dx2, dx, c1, c2, c3); + double result = + fputil::multiply_add(ex, /*log(2)*/ 0x1.62e42fefa39efp-1, LOG_F[p1] + p); + return result; } // Rounding tests for 2^hi * (mid + lo) when the output might be denormal. We From a4061e14c501876eaa2adfa3c9d6c36392434275 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Sat, 19 Apr 2025 16:30:44 +0800 Subject: [PATCH 12/17] [libc] Update COEFF. --- libc/src/math/generic/atanhf16.cpp | 12 +++--------- libc/src/math/generic/explogxf.h | 13 ++++++++----- 2 files changed, 11 insertions(+), 14 deletions(-) diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp index 7e73b69d3a66b..29743347144f5 100644 --- a/libc/src/math/generic/atanhf16.cpp +++ b/libc/src/math/generic/atanhf16.cpp @@ -22,17 +22,11 @@ namespace LIBC_NAMESPACE_DECL { -static constexpr size_t N_EXCEPTS = 4; +static constexpr size_t N_EXCEPTS = 1; static constexpr fputil::ExceptValues 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)) { @@ -76,9 +70,9 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 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| < 2^-16, this can be approximated by: + // When |x| < 2^-6 (0x2400U), this can be approximated by: // atanh(x) ≈ x + (1/3)*x^3 - if (LIBC_UNLIKELY(x_abs < 0x0100U)) { + if (LIBC_UNLIKELY(x_abs < 0x2400U)) { float xf = x; return fputil::cast(xf + 0x1.555556p-2f * xf * xf * xf); } diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h index f02806fd2b3f8..2310d7c3d2731 100644 --- a/libc/src/math/generic/explogxf.h +++ b/libc/src/math/generic/explogxf.h @@ -314,10 +314,13 @@ LIBC_INLINE static float log_eval_f(float x) { // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)). float dx = (xbits.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}; + // Minimax polynomial for log(1 + dx), generated using Sollya: + // > P = fpminimax(log(1 + x)/x, 6, [|SG...|], [0, 2^-7]); + // > Q = (P - 1) / x; + // > for i from 0 to degree(Q) do print(coeff(Q, i)); + static constexpr float COEFFS[6] = {-0x1p-1f, 0x1.555556p-2f, + -0x1.00022ep-2f, 0x1.9ea056p-3f, + -0x1.e50324p-2f, 0x1.c018fp3f}; float dx2 = dx * dx; @@ -330,7 +333,7 @@ LIBC_INLINE static float log_eval_f(float x) { // Generated by Sollya with the following commands: // > display = hexadecimal; // > round(log(2), SG, RN); - static constexpr float LOGF_2 = 0x1.62e43p-1f; + constexpr float LOGF_2 = 0x1.62e43p-1f; float result = fputil::multiply_add(ex, LOGF_2, LOG_F_FLOAT[p1] + p); return result; From fd7ff4c87651a0cb8eda3de7c5d412ee78abd10a Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Thu, 24 Apr 2025 15:00:12 +0000 Subject: [PATCH 13/17] [libc] Update range. --- libc/src/math/generic/atanhf16.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp index 29743347144f5..8d601d92f0910 100644 --- a/libc/src/math/generic/atanhf16.cpp +++ b/libc/src/math/generic/atanhf16.cpp @@ -62,22 +62,22 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 x)) { LIBC_UNLIKELY(r.has_value())) return r.value(); - // For |x| less than approximately 0.10 - if (LIBC_UNLIKELY(x_abs <= 0x2e66U)) { + // For |x| less than approximately 0.24 + if (LIBC_UNLIKELY(x_abs <= 0x33f3U)) { // 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| < 2^-6 (0x2400U), this can be approximated by: + // When |x| < 2^-7 (0x0200U), this can be approximated by: // atanh(x) ≈ x + (1/3)*x^3 - if (LIBC_UNLIKELY(x_abs < 0x2400U)) { + if (LIBC_UNLIKELY(x_abs < 0x0200U)) { float xf = x; return fputil::cast(xf + 0x1.555556p-2f * xf * xf * xf); } - // For 2^-16 <= |x| <= 0x1.998p-4 (~0.10): + // For 2^-7 <= |x| <= 0x1.fccp-3 (~0.24): // 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): From ba79fa0602ed998fea3fc28ada988e3ecf4aa5b2 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Thu, 24 Apr 2025 15:13:57 +0000 Subject: [PATCH 14/17] [libc] Update range again. --- libc/src/math/generic/atanhf16.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/libc/src/math/generic/atanhf16.cpp b/libc/src/math/generic/atanhf16.cpp index 8d601d92f0910..57885ac0b767f 100644 --- a/libc/src/math/generic/atanhf16.cpp +++ b/libc/src/math/generic/atanhf16.cpp @@ -70,14 +70,14 @@ LLVM_LIBC_FUNCTION(float16, atanhf16, (float16 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| < 2^-7 (0x0200U), this can be approximated by: + // When |x| < 2^-5 (0x0800U), this can be approximated by: // atanh(x) ≈ x + (1/3)*x^3 - if (LIBC_UNLIKELY(x_abs < 0x0200U)) { + if (LIBC_UNLIKELY(x_abs < 0x0800U)) { float xf = x; return fputil::cast(xf + 0x1.555556p-2f * xf * xf * xf); } - // For 2^-7 <= |x| <= 0x1.fccp-3 (~0.24): + // For 2^-5 <= |x| <= 0x1.fccp-3 (~0.24): // 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): From 822f67269824025fcb518a7d1857f4da1e9b8132 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Thu, 24 Apr 2025 15:43:34 +0000 Subject: [PATCH 15/17] [libc] Add to do comments. --- libc/src/math/generic/explogxf.h | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h index 2310d7c3d2731..8acdb7e75dfe2 100644 --- a/libc/src/math/generic/explogxf.h +++ b/libc/src/math/generic/explogxf.h @@ -298,6 +298,8 @@ LIBC_INLINE static double log2_eval(double x) { } // x should be positive, normal finite value +// TODO: Simplify range reduction and polynomial degree for float16. +// See issue #137190. LIBC_INLINE static float log_eval_f(float x) { // For x = 2^ex * (1 + mx), logf(x) = ex * logf(2) + logf(1 + mx). using FPBits = fputil::FPBits; @@ -318,9 +320,8 @@ LIBC_INLINE static float log_eval_f(float x) { // > P = fpminimax(log(1 + x)/x, 6, [|SG...|], [0, 2^-7]); // > Q = (P - 1) / x; // > for i from 0 to degree(Q) do print(coeff(Q, i)); - static constexpr float COEFFS[6] = {-0x1p-1f, 0x1.555556p-2f, - -0x1.00022ep-2f, 0x1.9ea056p-3f, - -0x1.e50324p-2f, 0x1.c018fp3f}; + constexpr float COEFFS[6] = {-0x1p-1f, 0x1.555556p-2f, -0x1.00022ep-2f, + 0x1.9ea056p-3f, -0x1.e50324p-2f, 0x1.c018fp3f}; float dx2 = dx * dx; From 03fc016c9b4cca77bd2ae8c7edaaba44f2294c70 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Fri, 25 Apr 2025 11:17:26 +0000 Subject: [PATCH 16/17] [libc] Update comments. --- libc/src/math/generic/common_constants.cpp | 2 +- libc/src/math/generic/explogxf.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/libc/src/math/generic/common_constants.cpp b/libc/src/math/generic/common_constants.cpp index a53641cf37bdb..b2c1293c6326d 100644 --- a/libc/src/math/generic/common_constants.cpp +++ b/libc/src/math/generic/common_constants.cpp @@ -98,7 +98,7 @@ 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, +// Lookup table for log(f) = log(1 + n*2^(-7)) where n = 0..127, // computed and stored as float precision constants. // Generated by Sollya with the following commands: // display = hexadecimal; diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h index 8acdb7e75dfe2..002b7659d4848 100644 --- a/libc/src/math/generic/explogxf.h +++ b/libc/src/math/generic/explogxf.h @@ -310,7 +310,7 @@ LIBC_INLINE static float log_eval_f(float x) { // p1 * 2^(-7) <= m_x < (p1 + 1) * 2^(-7). int p1 = static_cast(xbits.get_mantissa() >> (FPBits::FRACTION_LEN - 7)); - // Set bs to (1 + (mx - p1*2^(-7)) + // Set bits to (1 + (mx - p1*2^(-7)) xbits.set_uintval(xbits.uintval() & (FPBits::FRACTION_MASK >> 7)); xbits.set_biased_exponent(FPBits::EXP_BIAS); // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)). From 5d67867c674b1a430bdaa84f380a2790e40129f6 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Fri, 25 Apr 2025 11:19:16 +0000 Subject: [PATCH 17/17] [libc] Update comments again. --- libc/src/math/generic/explogxf.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/src/math/generic/explogxf.h b/libc/src/math/generic/explogxf.h index 002b7659d4848..212ede4758549 100644 --- a/libc/src/math/generic/explogxf.h +++ b/libc/src/math/generic/explogxf.h @@ -310,7 +310,7 @@ LIBC_INLINE static float log_eval_f(float x) { // p1 * 2^(-7) <= m_x < (p1 + 1) * 2^(-7). int p1 = static_cast(xbits.get_mantissa() >> (FPBits::FRACTION_LEN - 7)); - // Set bits to (1 + (mx - p1*2^(-7)) + // Set bits to (1 + (mx - p1*2^(-7))) xbits.set_uintval(xbits.uintval() & (FPBits::FRACTION_MASK >> 7)); xbits.set_biased_exponent(FPBits::EXP_BIAS); // dx = (mx - p1*2^(-7)) / (1 + p1*2^(-7)).