From 3922caef603af2708796f5641bbceef9fa3b3842 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Mon, 10 Mar 2025 20:04:11 +0800 Subject: [PATCH 01/17] [libc][math][c23] Add acoshf16 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/acoshf16.h | 21 +++++++++ libc/src/math/generic/CMakeLists.txt | 20 ++++++++ libc/src/math/generic/acoshf16.cpp | 55 ++++++++++++++++++++++ libc/test/src/math/CMakeLists.txt | 11 +++++ libc/test/src/math/acoshf16_test.cpp | 52 ++++++++++++++++++++ libc/test/src/math/smoke/CMakeLists.txt | 11 +++++ libc/test/src/math/smoke/acoshf16_test.cpp | 41 ++++++++++++++++ 11 files changed, 221 insertions(+), 1 deletion(-) create mode 100644 libc/src/math/acoshf16.h create mode 100644 libc/src/math/generic/acoshf16.cpp create mode 100644 libc/test/src/math/acoshf16_test.cpp create mode 100644 libc/test/src/math/smoke/acoshf16_test.cpp diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index 3cb9ee82752b3..cfb42344e7953 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -653,6 +653,7 @@ if(LIBC_TYPES_HAS_FLOAT16) # math.h C23 _Float16 entrypoints libc.src.math.asinf16 libc.src.math.acosf16 + libc.src.math.acoshf16 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 5b855ce4881c3..e137626361f48 100644 --- a/libc/docs/headers/math/index.rst +++ b/libc/docs/headers/math/index.rst @@ -251,7 +251,7 @@ Higher Math Functions +===========+==================+=================+========================+======================+========================+========================+============================+ | acos | |check| | | | |check| | | 7.12.4.1 | F.10.1.1 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| acosh | |check| | | | | | 7.12.5.1 | F.10.2.1 | +| acosh | |check| | | | |check| | | 7.12.5.1 | F.10.2.1 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | acospi | | | | | | 7.12.4.8 | F.10.1.8 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/include/math.yaml b/libc/include/math.yaml index a66f981030864..7fd73445cc74e 100644 --- a/libc/include/math.yaml +++ b/libc/include/math.yaml @@ -27,6 +27,13 @@ functions: return_type: float arguments: - type: float + name: acoshf16 + standards: + - stdc + return_type: _Float16 + arguments: + - type: _Float16 + guard: LIBC_TYPES_HAS_FLOAT16 - name: asin standards: - stdc diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index f18a73d46f9aa..bdf8b923416d8 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -46,6 +46,7 @@ add_math_entrypoint_object(acosf16) add_math_entrypoint_object(acosh) add_math_entrypoint_object(acoshf) +add_math_entrypoint_object(acoshf16) add_math_entrypoint_object(asin) add_math_entrypoint_object(asinf) diff --git a/libc/src/math/acoshf16.h b/libc/src/math/acoshf16.h new file mode 100644 index 0000000000000..a23426ffe589d --- /dev/null +++ b/libc/src/math/acoshf16.h @@ -0,0 +1,21 @@ +//===-- Implementation header for acoshf16 -----------------------*- 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_ACOSHF16_H +#define LLVM_LIBC_SRC_MATH_ACOSHF16_H + +#include "src/__support/macros/config.h" +#include "src/__support/macros/properties/types.h" + +namespace LIBC_NAMESPACE_DECL { + +float16 acoshf16(float16 x); + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC_MATH_ACOSHF16_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 3114289bad486..be3634dfebdc7 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3944,6 +3944,26 @@ add_entrypoint_object( libc.src.__support.macros.optimization ) +add_entrypoint_object( + acoshf16 + SRCS + acoshf16.cpp + HDRS + ../acoshf16.h + DEPENDS + 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 + libc.src.__support.FPUtil.polyeval + libc.src.__support.FPUtil.sqrt + libc.src.__support.macros.optimization + libc.src.__support.macros.properties.types +) + add_entrypoint_object( asinhf SRCS diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp new file mode 100644 index 0000000000000..01f97f536c256 --- /dev/null +++ b/libc/src/math/generic/acoshf16.cpp @@ -0,0 +1,55 @@ +//===-- Half-precision acoshf16 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/acoshf16.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/sqrt.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/cast.h" +#include "src/__support/macros/optimization.h" +#include "src/math/generic/explogxf.h" + +namespace LIBC_NAMESPACE_DECL { + +LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { + using FPBits = fputil::FPBits; + FPBits xbits(x); + + uint16_t x_u = xbits.uintval(); + uint16_t x_abs = x_u & 0x7fff; + + // acoshf16(x) domain: x >= 1.0 + if (LIBC_UNLIKELY(x_abs < 0x3c00)) { + fputil::set_errno_if_required(EDOM); + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + + if (LIBC_UNLIKELY(x == float16(1.0f))) + return float16(0.0f); + + // acosh(inf) = inf, acosh(NaN) = NaN + if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) + return x; + + // Compute in float32 for accuracy + float xf32 = static_cast(x); + // sqrt(x^2 - 1) + float sqrt_term = fputil::sqrt( + fputil::multiply_add(xf32, xf32, -1.0f)); + + // log(x + sqrt(x^2 - 1)) + float result = static_cast(log_eval(xf32 + sqrt_term)); + + return fputil::cast(result); +} + +} // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index 032050bb06ec3..5dd53de0234a7 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -2173,6 +2173,17 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + acoshf16_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + acoshf16_test.cpp + DEPENDS + libc.src.math.acoshf16 +) + add_fp_unittest( asinf_test NEED_MPFR diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp new file mode 100644 index 0000000000000..d4ebef87f3b61 --- /dev/null +++ b/libc/test/src/math/acoshf16_test.cpp @@ -0,0 +1,52 @@ +//===-- Unittests for acoshf16 ---------------------------------------------===// +// +// 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/acoshf16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +using LlvmLibcAcoshf16Test = LIBC_NAMESPACE::testing::FPTest; +namespace mpfr = LIBC_NAMESPACE::testing::mpfr; + +TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { + LIBC_NAMESPACE::libc_errno = 0; + + // NaN input + EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(aNaN)); + EXPECT_MATH_ERRNO(0); + + // acoshf16(1.0) = 0 + EXPECT_FP_EQ_ALL_ROUNDING(float16(0.0f), LIBC_NAMESPACE::acoshf16(float16(1.0f))); + EXPECT_MATH_ERRNO(0); + + // Domain error (x < 1) + EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(float16(0.5f))); + EXPECT_MATH_ERRNO(EDOM); + + // acoshf16(+inf) = inf + EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::acoshf16(inf)); + EXPECT_MATH_ERRNO(0); + + // acoshf16(x) domain error (negative infinity) + EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(neg_inf)); + EXPECT_MATH_ERRNO(EDOM); +} + +TEST_F(LlvmLibcAcoshf16Test, InFloat16Range) { + constexpr uint16_t START = 0x3c00U; // 1.0 + constexpr uint16_t STOP = 0x7bffU; // Largest finite float16 value + + for (uint16_t bits = START; bits <= STOP; ++bits) { + float16 x = FPBits(bits).get_val(); + ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, float(x), + float(LIBC_NAMESPACE::acoshf16(x)), 0.5); + } +} diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index aef7c83ba0215..2c152e39e1c39 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -3945,6 +3945,17 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + acoshf16_test + SUITE + libc-math-smoke-tests + SRCS + acoshf16_test.cpp + DEPENDS + libc.src.errno.errno + libc.src.math.acoshf16 +) + add_fp_unittest( asinf_test SUITE diff --git a/libc/test/src/math/smoke/acoshf16_test.cpp b/libc/test/src/math/smoke/acoshf16_test.cpp new file mode 100644 index 0000000000000..e197f3f305834 --- /dev/null +++ b/libc/test/src/math/smoke/acoshf16_test.cpp @@ -0,0 +1,41 @@ +//===-- Unittests for acoshf16 --------------------------------------------===// +// +// 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/errno/libc_errno.h" +#include "src/math/acoshf16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" + +using LlvmLibcAcoshf16Test = LIBC_NAMESPACE::testing::FPTest; + +TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { + LIBC_NAMESPACE::libc_errno = 0; + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acoshf16(sNaN), FE_INVALID); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(zero)); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(neg_zero)); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ(inf, LIBC_NAMESPACE::acoshf16(inf)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(neg_inf)); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ(float16(0.0f), LIBC_NAMESPACE::acoshf16(float16(1.0f))); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(float16(0.5f))); + EXPECT_MATH_ERRNO(EDOM); +} \ No newline at end of file From 8906a102bf60f64b8438a6cb6009e582ddcf5ec3 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Mon, 10 Mar 2025 21:57:22 +0800 Subject: [PATCH 02/17] [libc][math][c23] Fix some build issues. --- libc/src/math/generic/CMakeLists.txt | 2 + libc/src/math/generic/acoshf16.cpp | 58 +++++++++++++++++++++------- libc/test/src/math/acoshf16_test.cpp | 51 +++++++++++++++++------- 3 files changed, 83 insertions(+), 28 deletions(-) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index be3634dfebdc7..db5a23c199eea 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3951,6 +3951,8 @@ add_entrypoint_object( HDRS ../acoshf16.h DEPENDS + .common_constants + .explogxf libc.hdr.errno_macros libc.hdr.fenv_macros libc.src.__support.FPUtil.cast diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index 01f97f536c256..dc8bdebc90abb 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -11,10 +11,11 @@ #include "hdr/fenv_macros.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/sqrt.h" -#include "src/__support/FPUtil/multiply_add.h" #include "src/__support/FPUtil/cast.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/sqrt.h" #include "src/__support/macros/optimization.h" +#include "src/math/generic/common_constants.h" #include "src/math/generic/explogxf.h" namespace LIBC_NAMESPACE_DECL { @@ -26,30 +27,57 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { uint16_t x_u = xbits.uintval(); uint16_t x_abs = x_u & 0x7fff; - // acoshf16(x) domain: x >= 1.0 + // Check for NaN input first. + if (LIBC_UNLIKELY(xbits.is_nan())) { + if (xbits.is_signaling_nan()) { + fputil::raise_except(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + return x; + } + + // Check for infinite inputs. + if (LIBC_UNLIKELY(xbits.is_inf())) { + if (xbits.is_neg()) { + fputil::set_errno_if_required(EDOM); + fputil::raise_except(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + return x; + } + + // Domain error for inputs less than 1.0. if (LIBC_UNLIKELY(x_abs < 0x3c00)) { fputil::set_errno_if_required(EDOM); - fputil::raise_except_if_required(FE_INVALID); + fputil::raise_except(FE_INVALID); return FPBits::quiet_nan().get_val(); } - if (LIBC_UNLIKELY(x == float16(1.0f))) + // acosh(1.0) exactly equals 0.0 + if (LIBC_UNLIKELY(xbits.uintval() == 0x3c00U)) return float16(0.0f); - // acosh(inf) = inf, acosh(NaN) = NaN - if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) - return x; - - // Compute in float32 for accuracy float xf32 = static_cast(x); - // sqrt(x^2 - 1) - float sqrt_term = fputil::sqrt( - fputil::multiply_add(xf32, xf32, -1.0f)); - // log(x + sqrt(x^2 - 1)) + // High precision for inputs very close to 1.0 + if (LIBC_UNLIKELY(xf32 < 1.25f)) { + float delta = xf32 - 1.0f; + float sqrt_2_delta = fputil::sqrt(2.0f * delta); + double correction = (double)delta / 12.0 - (3.0 * (double)delta * delta) / 160.0; + double precise_result = (double)sqrt_2_delta * (1.0 - correction); + return fputil::cast(static_cast(precise_result)); + } + // Special optimization for large input values. + if (LIBC_UNLIKELY(xf32 >= 32.0f)) { + float result = static_cast(log_eval(2.0f * xf32)); + return fputil::cast(result); + } + + // Standard computation for general case. + float sqrt_term = fputil::sqrt(fputil::multiply_add(xf32, xf32, -1.0f)); float result = static_cast(log_eval(xf32 + sqrt_term)); return fputil::cast(result); } -} // namespace LIBC_NAMESPACE_DECL +} // namespace LIBC_NAMESPACE_DECL \ No newline at end of file diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index d4ebef87f3b61..099486ba338c0 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -1,17 +1,17 @@ -//===-- Unittests for acoshf16 ---------------------------------------------===// +//===-- Unittests for acoshf16 ----------------------------------------------===// // -// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// 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/acoshf16.h" +#include "src/__support/FPUtil/FPBits.h" #include "test/UnitTest/FPMatcher.h" -#include "test/UnitTest/Test.h" #include "utils/MPFRWrapper/MPFRUtils.h" +#include using LlvmLibcAcoshf16Test = LIBC_NAMESPACE::testing::FPTest; namespace mpfr = LIBC_NAMESPACE::testing::mpfr; @@ -19,34 +19,59 @@ namespace mpfr = LIBC_NAMESPACE::testing::mpfr; TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { LIBC_NAMESPACE::libc_errno = 0; - // NaN input EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(aNaN)); EXPECT_MATH_ERRNO(0); - // acoshf16(1.0) = 0 - EXPECT_FP_EQ_ALL_ROUNDING(float16(0.0f), LIBC_NAMESPACE::acoshf16(float16(1.0f))); + EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acoshf16(sNaN), FE_INVALID); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ_ALL_ROUNDING(float16(0.0f), + LIBC_NAMESPACE::acoshf16(float16(1.0f))); EXPECT_MATH_ERRNO(0); - // Domain error (x < 1) EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(float16(0.5f))); EXPECT_MATH_ERRNO(EDOM); - // acoshf16(+inf) = inf EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::acoshf16(inf)); EXPECT_MATH_ERRNO(0); - // acoshf16(x) domain error (negative infinity) EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(neg_inf)); EXPECT_MATH_ERRNO(EDOM); } TEST_F(LlvmLibcAcoshf16Test, InFloat16Range) { - constexpr uint16_t START = 0x3c00U; // 1.0 - constexpr uint16_t STOP = 0x7bffU; // Largest finite float16 value + constexpr uint16_t START = 0x3C00U; // 1.0 + constexpr uint16_t STOP = 0x7BFFU; // Largest finite float16 value for (uint16_t bits = START; bits <= STOP; ++bits) { float16 x = FPBits(bits).get_val(); + if (FPBits(bits).is_nan() || FPBits(bits).is_inf()) + continue; ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, float(x), - float(LIBC_NAMESPACE::acoshf16(x)), 0.5); + float(LIBC_NAMESPACE::acoshf16(x)), + 5.0); + } +} + +TEST_F(LlvmLibcAcoshf16Test, SpecificBitPatterns) { + constexpr int N = 12; + constexpr uint16_t INPUTS[N] = { + 0x3C00, // 1.0 + 0x3C01, // just above 1.0 (minimally larger than 1) + 0x3E00, // 1.5 + 0x4200, // 3.0 + 0x4500, // 5.0 + 0x4900, // 10.0 + 0x51FF, // ~47.94 (random mid-range value) + 0x5CB0, // ~300.0 (random mid-range value) + 0x643F, // ~1087.6 (random large value) + 0x77FF, // just below next exponent interval (max for exponent 0x1D) + 0x7801, // just above previous value (min for exponent 0x1E) + 0x7BFF // 65504.0 (max finite half) + }; + for (int i = 0; i < N; ++i) { + float16 x = FPBits(INPUTS[i]).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, x, + LIBC_NAMESPACE::acoshf16(x), 0.5); } } From e9eca1b845039bdab11bbacc3b7850053fdd8dc7 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 11 Mar 2025 11:13:10 +0800 Subject: [PATCH 03/17] [libc][math][c23] Fix some test failures. --- libc/src/math/generic/acoshf16.cpp | 25 ++++++++++++---- libc/test/src/math/acoshf16_test.cpp | 45 +++++++++++++++------------- 2 files changed, 45 insertions(+), 25 deletions(-) diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index dc8bdebc90abb..8d629bad55c92 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -23,10 +23,18 @@ namespace LIBC_NAMESPACE_DECL { LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { using FPBits = fputil::FPBits; FPBits xbits(x); - uint16_t x_u = xbits.uintval(); uint16_t x_abs = x_u & 0x7fff; +// if (LIBC_UNLIKELY(x <= 1.0f)) { +// if (x == 1.0f) +// return 0.0f; +// // x < 1. +// fputil::set_errno_if_required(EDOM); +// fputil::raise_except_if_required(FE_INVALID); +// return FPBits::quiet_nan().get_val(); +// } + // Check for NaN input first. if (LIBC_UNLIKELY(xbits.is_nan())) { if (xbits.is_signaling_nan()) { @@ -62,11 +70,18 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { // High precision for inputs very close to 1.0 if (LIBC_UNLIKELY(xf32 < 1.25f)) { float delta = xf32 - 1.0f; - float sqrt_2_delta = fputil::sqrt(2.0f * delta); - double correction = (double)delta / 12.0 - (3.0 * (double)delta * delta) / 160.0; - double precise_result = (double)sqrt_2_delta * (1.0 - correction); - return fputil::cast(static_cast(precise_result)); + float sqrt_2 = fputil::sqrt(2.0f * delta); + float sqrt_2d = fputil::sqrt(2.0f * delta); + float d32 = delta * fputil::sqrt(delta); + float term2 = d32 / (6.0f * fputil::sqrt(2.0f)); + float d52 = d32 * delta; + float term3 = 3.0f * d52 / (80.0f * sqrt_2); + float d72 = d52 * delta; + float term4 = 5.0f * d72 / (1792.0f * sqrt_2); + float result = sqrt_2d - term2 + term3 - term4; + return fputil::cast(result); } + // Special optimization for large input values. if (LIBC_UNLIKELY(xf32 >= 32.0f)) { float result = static_cast(log_eval(2.0f * xf32)); diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index 099486ba338c0..d13e8f693f80b 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -6,6 +6,7 @@ // //===----------------------------------------------------------------------===// +#include "src/__support/FPUtil/cast.h" #include "src/errno/libc_errno.h" #include "src/math/acoshf16.h" #include "src/__support/FPUtil/FPBits.h" @@ -40,34 +41,38 @@ TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { } TEST_F(LlvmLibcAcoshf16Test, InFloat16Range) { - constexpr uint16_t START = 0x3C00U; // 1.0 - constexpr uint16_t STOP = 0x7BFFU; // Largest finite float16 value + constexpr uint32_t COUNT = 100'000; + constexpr uint32_t STEP = UINT32_MAX / COUNT; - for (uint16_t bits = START; bits <= STOP; ++bits) { - float16 x = FPBits(bits).get_val(); - if (FPBits(bits).is_nan() || FPBits(bits).is_inf()) + for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) { + LIBC_NAMESPACE::fputil::FPBits bits(v); + float xf32 = bits.get_val(); + if (bits.is_nan() || bits.is_inf()) continue; - ASSERT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, float(x), - float(LIBC_NAMESPACE::acoshf16(x)), - 5.0); + if (xf32 < 1.0f) + continue; + float16 xh = LIBC_NAMESPACE::fputil::cast(xf32); + + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, xh, + LIBC_NAMESPACE::acoshf16(xh), 3.0); } } TEST_F(LlvmLibcAcoshf16Test, SpecificBitPatterns) { constexpr int N = 12; constexpr uint16_t INPUTS[N] = { - 0x3C00, // 1.0 - 0x3C01, // just above 1.0 (minimally larger than 1) - 0x3E00, // 1.5 - 0x4200, // 3.0 - 0x4500, // 5.0 - 0x4900, // 10.0 - 0x51FF, // ~47.94 (random mid-range value) - 0x5CB0, // ~300.0 (random mid-range value) - 0x643F, // ~1087.6 (random large value) - 0x77FF, // just below next exponent interval (max for exponent 0x1D) - 0x7801, // just above previous value (min for exponent 0x1E) - 0x7BFF // 65504.0 (max finite half) + 0x3C00, // x = 1.0 + 0x3C01, // x = just above 1.0 (minimally larger than 1) + 0x3E00, // x = 1.5 + 0x4200, // x = 3.0 + 0x4500, // x = 5.0 + 0x4900, // x = 10.0 + 0x51FF, // x = ~47.94 + 0x5CB0, // x = ~300.0 + 0x643F, // x = ~1087.6 + 0x77FF, // x = just below next exponent interval (max for exponent 0x1D) + 0x7801, // x = just above previous value (min for exponent 0x1E) + 0x7BFF // x = 65504.0 (max finite half) }; for (int i = 0; i < N; ++i) { float16 x = FPBits(INPUTS[i]).get_val(); From c6d781de0c7adf6e1fe93a7d315209b497498946 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 11 Mar 2025 20:19:01 +0800 Subject: [PATCH 04/17] [libc][math][c23] Do some change for comments. --- libc/docs/headers/math/index.rst | 2 +- libc/src/math/acoshf16.h | 2 +- libc/src/math/generic/acoshf16.cpp | 46 ++++++++++++++-------------- libc/test/src/math/acoshf16_test.cpp | 22 +++++-------- 4 files changed, 33 insertions(+), 39 deletions(-) diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst index e137626361f48..8dc7c8c1a0966 100644 --- a/libc/docs/headers/math/index.rst +++ b/libc/docs/headers/math/index.rst @@ -251,7 +251,7 @@ Higher Math Functions +===========+==================+=================+========================+======================+========================+========================+============================+ | acos | |check| | | | |check| | | 7.12.4.1 | F.10.1.1 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| acosh | |check| | | | |check| | | 7.12.5.1 | F.10.2.1 | +| acosh | |check| | | | |check| | | 7.12.5.1 | F.10.2.1 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | acospi | | | | | | 7.12.4.8 | F.10.1.8 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/src/math/acoshf16.h b/libc/src/math/acoshf16.h index a23426ffe589d..f471ecfeb6154 100644 --- a/libc/src/math/acoshf16.h +++ b/libc/src/math/acoshf16.h @@ -1,4 +1,4 @@ -//===-- Implementation header for acoshf16 -----------------------*- C++ -*-===// +//===-- Implementation header for acoshf16 ----------------------*- C++ -*-===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index 8d629bad55c92..8ebaabdf670ad 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -1,4 +1,4 @@ -//===-- Half-precision acoshf16 function -----------------------------------===// +//===-- Half-precision acoshf16 function ----------------------------------===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. @@ -38,7 +38,7 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { // Check for NaN input first. if (LIBC_UNLIKELY(xbits.is_nan())) { if (xbits.is_signaling_nan()) { - fputil::raise_except(FE_INVALID); + fputil::raise_except_if_required(FE_INVALID); return FPBits::quiet_nan().get_val(); } return x; @@ -48,7 +48,7 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { if (LIBC_UNLIKELY(xbits.is_inf())) { if (xbits.is_neg()) { fputil::set_errno_if_required(EDOM); - fputil::raise_except(FE_INVALID); + fputil::raise_except_if_required(FE_INVALID); return FPBits::quiet_nan().get_val(); } return x; @@ -57,36 +57,36 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { // Domain error for inputs less than 1.0. if (LIBC_UNLIKELY(x_abs < 0x3c00)) { fputil::set_errno_if_required(EDOM); - fputil::raise_except(FE_INVALID); + fputil::raise_except_if_required(FE_INVALID); return FPBits::quiet_nan().get_val(); } // acosh(1.0) exactly equals 0.0 - if (LIBC_UNLIKELY(xbits.uintval() == 0x3c00U)) + if (LIBC_UNLIKELY(x_u == 0x3c00U)) return float16(0.0f); float xf32 = static_cast(x); // High precision for inputs very close to 1.0 - if (LIBC_UNLIKELY(xf32 < 1.25f)) { - float delta = xf32 - 1.0f; - float sqrt_2 = fputil::sqrt(2.0f * delta); - float sqrt_2d = fputil::sqrt(2.0f * delta); - float d32 = delta * fputil::sqrt(delta); - float term2 = d32 / (6.0f * fputil::sqrt(2.0f)); - float d52 = d32 * delta; - float term3 = 3.0f * d52 / (80.0f * sqrt_2); - float d72 = d52 * delta; - float term4 = 5.0f * d72 / (1792.0f * sqrt_2); - float result = sqrt_2d - term2 + term3 - term4; - return fputil::cast(result); - } + // if (LIBC_UNLIKELY(xf32 < 1.25f)) { + // float delta = xf32 - 1.0f; + // float sqrt_2 = fputil::sqrt(2.0f * delta); + // float sqrt_2d = fputil::sqrt(2.0f * delta); + // float d32 = delta * fputil::sqrt(delta); + // float term2 = d32 / (6.0f * fputil::sqrt(2.0f)); + // float d52 = d32 * delta; + // float term3 = 3.0f * d52 / (80.0f * sqrt_2); + // float d72 = d52 * delta; + // float term4 = 5.0f * d72 / (1792.0f * sqrt_2); + // float result = sqrt_2d - term2 + term3 - term4; + // return fputil::cast(result); + // } // Special optimization for large input values. - if (LIBC_UNLIKELY(xf32 >= 32.0f)) { - float result = static_cast(log_eval(2.0f * xf32)); - return fputil::cast(result); - } + // if (LIBC_UNLIKELY(xf32 >= 32.0f)) { + // float result = static_cast(log_eval(2.0f * xf32)); + // return fputil::cast(result); + // } // Standard computation for general case. float sqrt_term = fputil::sqrt(fputil::multiply_add(xf32, xf32, -1.0f)); @@ -95,4 +95,4 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { return fputil::cast(result); } -} // namespace LIBC_NAMESPACE_DECL \ No newline at end of file +} // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index d13e8f693f80b..2e457815debcd 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -17,6 +17,9 @@ using LlvmLibcAcoshf16Test = LIBC_NAMESPACE::testing::FPTest; namespace mpfr = LIBC_NAMESPACE::testing::mpfr; +static constexpr uint16_t START = 0x3c00U; +static constexpr uint16_t STOP = 0x7bffU; + TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { LIBC_NAMESPACE::libc_errno = 0; @@ -40,21 +43,12 @@ TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { EXPECT_MATH_ERRNO(EDOM); } -TEST_F(LlvmLibcAcoshf16Test, InFloat16Range) { - constexpr uint32_t COUNT = 100'000; - constexpr uint32_t STEP = UINT32_MAX / COUNT; - - for (uint32_t i = 0, v = 0; i <= COUNT; ++i, v += STEP) { - LIBC_NAMESPACE::fputil::FPBits bits(v); - float xf32 = bits.get_val(); - if (bits.is_nan() || bits.is_inf()) - continue; - if (xf32 < 1.0f) - continue; - float16 xh = LIBC_NAMESPACE::fputil::cast(xf32); +TEST_F(LlvmLibcAcoshf16Test, PositiveRange) { + for (uint16_t v = START; v <= STOP; ++v) { + float16 x = FPBits(v).get_val(); - EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, xh, - LIBC_NAMESPACE::acoshf16(xh), 3.0); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, x, + LIBC_NAMESPACE::acoshf16(x), 3.0); } } From c6582d9b5e0f20848814c25c375b669c31491327 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 11 Mar 2025 20:35:32 +0800 Subject: [PATCH 05/17] [libc][math][c23] Do some change for comments again. --- libc/src/math/generic/CMakeLists.txt | 1 - libc/src/math/generic/acoshf16.cpp | 4 ++-- libc/test/src/math/acoshf16_test.cpp | 2 +- libc/test/src/math/smoke/acoshf16_test.cpp | 2 +- 4 files changed, 4 insertions(+), 5 deletions(-) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index db5a23c199eea..fd5ecd72b765f 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3956,7 +3956,6 @@ 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/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index 8ebaabdf670ad..be596873579d1 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -55,7 +55,7 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { } // Domain error for inputs less than 1.0. - if (LIBC_UNLIKELY(x_abs < 0x3c00)) { + if (LIBC_UNLIKELY(x_abs < 0x3c00U)) { fputil::set_errno_if_required(EDOM); fputil::raise_except_if_required(FE_INVALID); return FPBits::quiet_nan().get_val(); @@ -65,7 +65,7 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { if (LIBC_UNLIKELY(x_u == 0x3c00U)) return float16(0.0f); - float xf32 = static_cast(x); + float xf32 = x; // High precision for inputs very close to 1.0 // if (LIBC_UNLIKELY(xf32 < 1.25f)) { diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index 2e457815debcd..d2883abc0f350 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -1,4 +1,4 @@ -//===-- Unittests for acoshf16 ----------------------------------------------===// +//===-- Unittests for acoshf16 --------------------------------------------===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. diff --git a/libc/test/src/math/smoke/acoshf16_test.cpp b/libc/test/src/math/smoke/acoshf16_test.cpp index e197f3f305834..b8d13b9fe88c3 100644 --- a/libc/test/src/math/smoke/acoshf16_test.cpp +++ b/libc/test/src/math/smoke/acoshf16_test.cpp @@ -38,4 +38,4 @@ TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(float16(0.5f))); EXPECT_MATH_ERRNO(EDOM); -} \ No newline at end of file +} From 83ea70e7b87a91273eb00c9623f867ad102f8454 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 11 Mar 2025 20:50:15 +0800 Subject: [PATCH 06/17] [libc][math][c23] Fix clang format issue. --- libc/src/math/generic/acoshf16.cpp | 21 +++++++++++---------- libc/test/src/math/acoshf16_test.cpp | 6 +++--- 2 files changed, 14 insertions(+), 13 deletions(-) diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index be596873579d1..162ec9cd1bbb4 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -26,15 +26,15 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { uint16_t x_u = xbits.uintval(); uint16_t x_abs = x_u & 0x7fff; -// if (LIBC_UNLIKELY(x <= 1.0f)) { -// if (x == 1.0f) -// return 0.0f; -// // x < 1. -// fputil::set_errno_if_required(EDOM); -// fputil::raise_except_if_required(FE_INVALID); -// return FPBits::quiet_nan().get_val(); -// } - + // if (LIBC_UNLIKELY(x <= 1.0f)) { + // if (x == 1.0f) + // return 0.0f; + // // x < 1. + // fputil::set_errno_if_required(EDOM); + // fputil::raise_except_if_required(FE_INVALID); + // return FPBits::quiet_nan().get_val(); + // } + // Check for NaN input first. if (LIBC_UNLIKELY(xbits.is_nan())) { if (xbits.is_signaling_nan()) { @@ -89,7 +89,8 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { // } // Standard computation for general case. - float sqrt_term = fputil::sqrt(fputil::multiply_add(xf32, xf32, -1.0f)); + float sqrt_term = + fputil::sqrt(fputil::multiply_add(xf32, xf32, -1.0f)); float result = static_cast(log_eval(xf32 + sqrt_term)); return fputil::cast(result); diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index d2883abc0f350..84574a565b9be 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -1,15 +1,15 @@ //===-- Unittests for acoshf16 --------------------------------------------===// // -// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// 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/__support/FPUtil/cast.h" #include "src/errno/libc_errno.h" #include "src/math/acoshf16.h" -#include "src/__support/FPUtil/FPBits.h" #include "test/UnitTest/FPMatcher.h" #include "utils/MPFRWrapper/MPFRUtils.h" #include @@ -18,7 +18,7 @@ using LlvmLibcAcoshf16Test = LIBC_NAMESPACE::testing::FPTest; namespace mpfr = LIBC_NAMESPACE::testing::mpfr; static constexpr uint16_t START = 0x3c00U; -static constexpr uint16_t STOP = 0x7bffU; +static constexpr uint16_t STOP = 0x7bffU; TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { LIBC_NAMESPACE::libc_errno = 0; From 59a1dc3812e72d9f6ad259dc5e1a4433965904a8 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Wed, 12 Mar 2025 21:57:24 +0800 Subject: [PATCH 07/17] [libc][math][c23] Reduce MPFR to 1. --- libc/src/math/generic/acoshf16.cpp | 51 +++++++++++----------------- libc/test/src/math/acoshf16_test.cpp | 2 +- 2 files changed, 21 insertions(+), 32 deletions(-) diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index 162ec9cd1bbb4..c133026f8f445 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -7,34 +7,28 @@ //===----------------------------------------------------------------------===// #include "src/math/acoshf16.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/cast.h" +#include "src/__support/FPUtil/except_value_utils.h" +#include "src/__support/FPUtil/generic/sqrt.h" #include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/FPUtil/sqrt.h" #include "src/__support/macros/optimization.h" -#include "src/math/generic/common_constants.h" #include "src/math/generic/explogxf.h" namespace LIBC_NAMESPACE_DECL { +static constexpr size_t N_EXCEPTS = 1; +static constexpr fputil::ExceptValues ACOSHF16_EXCEPTS{ + {// (input, RZ output, RU offset, RD offset, RN offset) + {0x41B7, 0x3ED8, 0, 1, 0}}}; + LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { using FPBits = fputil::FPBits; FPBits xbits(x); uint16_t x_u = xbits.uintval(); uint16_t x_abs = x_u & 0x7fff; - // if (LIBC_UNLIKELY(x <= 1.0f)) { - // if (x == 1.0f) - // return 0.0f; - // // x < 1. - // fputil::set_errno_if_required(EDOM); - // fputil::raise_except_if_required(FE_INVALID); - // return FPBits::quiet_nan().get_val(); - // } - // Check for NaN input first. if (LIBC_UNLIKELY(xbits.is_nan())) { if (xbits.is_signaling_nan()) { @@ -68,25 +62,20 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { float xf32 = x; // High precision for inputs very close to 1.0 - // if (LIBC_UNLIKELY(xf32 < 1.25f)) { - // float delta = xf32 - 1.0f; - // float sqrt_2 = fputil::sqrt(2.0f * delta); - // float sqrt_2d = fputil::sqrt(2.0f * delta); - // float d32 = delta * fputil::sqrt(delta); - // float term2 = d32 / (6.0f * fputil::sqrt(2.0f)); - // float d52 = d32 * delta; - // float term3 = 3.0f * d52 / (80.0f * sqrt_2); - // float d72 = d52 * delta; - // float term4 = 5.0f * d72 / (1792.0f * sqrt_2); - // float result = sqrt_2d - term2 + term3 - term4; - // return fputil::cast(result); - // } + if (LIBC_UNLIKELY(xf32 < 1.25f)) { + float delta = xf32 - 1.0f; + float sqrt_2_delta = fputil::sqrt(2.0 * delta); + float x2 = delta; + float pe = fputil::polyeval(x2, 0x1.0000000000000p+0f, + -0x1.55551a83a9472p-4f, 0x1.331601c4b8ecfp-6f, + -0x1.6890f49eb0acbp-8f, 0x1.8f3a617040a6ap-10f); + float approx = sqrt_2_delta * pe; + return fputil::cast(approx); + } - // Special optimization for large input values. - // if (LIBC_UNLIKELY(xf32 >= 32.0f)) { - // float result = static_cast(log_eval(2.0f * xf32)); - // return fputil::cast(result); - // } + if (auto r = ACOSHF16_EXCEPTS.lookup(xbits.uintval()); + LIBC_UNLIKELY(r.has_value())) + return r.value(); // Standard computation for general case. float sqrt_term = diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index 84574a565b9be..a68b0742be03e 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -48,7 +48,7 @@ TEST_F(LlvmLibcAcoshf16Test, PositiveRange) { float16 x = FPBits(v).get_val(); EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, x, - LIBC_NAMESPACE::acoshf16(x), 3.0); + LIBC_NAMESPACE::acoshf16(x), 1.0); } } From 1fe434ec31b6ccc7e7e11f706e34539a1810acfa Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Thu, 13 Mar 2025 12:51:55 +0000 Subject: [PATCH 08/17] [libc][math][c23] Reduce MPFR to 0.5. --- libc/src/math/generic/acoshf16.cpp | 14 +++++++------- libc/test/src/math/acoshf16_test.cpp | 2 +- 2 files changed, 8 insertions(+), 8 deletions(-) diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index c133026f8f445..1a4f68987ee6b 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -18,10 +18,11 @@ namespace LIBC_NAMESPACE_DECL { -static constexpr size_t N_EXCEPTS = 1; +static constexpr size_t N_EXCEPTS = 2; static constexpr fputil::ExceptValues ACOSHF16_EXCEPTS{ {// (input, RZ output, RU offset, RD offset, RN offset) - {0x41B7, 0x3ED8, 0, 1, 0}}}; + {0x41B7, 0x3ED8, 1, 0, 0}, + {0x3CE4, 0x393E, 1, 0, 1}}}; LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { using FPBits = fputil::FPBits; @@ -59,8 +60,11 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { if (LIBC_UNLIKELY(x_u == 0x3c00U)) return float16(0.0f); - float xf32 = x; + if (auto r = ACOSHF16_EXCEPTS.lookup(xbits.uintval()); + LIBC_UNLIKELY(r.has_value())) + return r.value(); + float xf32 = x; // High precision for inputs very close to 1.0 if (LIBC_UNLIKELY(xf32 < 1.25f)) { float delta = xf32 - 1.0f; @@ -73,10 +77,6 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { return fputil::cast(approx); } - if (auto r = ACOSHF16_EXCEPTS.lookup(xbits.uintval()); - LIBC_UNLIKELY(r.has_value())) - return r.value(); - // Standard computation for general case. float sqrt_term = fputil::sqrt(fputil::multiply_add(xf32, xf32, -1.0f)); diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index a68b0742be03e..a779d3110f329 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -48,7 +48,7 @@ TEST_F(LlvmLibcAcoshf16Test, PositiveRange) { float16 x = FPBits(v).get_val(); EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, x, - LIBC_NAMESPACE::acoshf16(x), 1.0); + LIBC_NAMESPACE::acoshf16(x), 0.5); } } From 02ecb340dba375eedce71a338a84f8dc9cf8e6ea Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Fri, 14 Mar 2025 15:13:23 +0000 Subject: [PATCH 09/17] [libc][math][c23] Update accuracy. --- libc/src/math/generic/acoshf16.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index 1a4f68987ee6b..c585c5079e33c 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -70,9 +70,9 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { float delta = xf32 - 1.0f; float sqrt_2_delta = fputil::sqrt(2.0 * delta); float x2 = delta; - float pe = fputil::polyeval(x2, 0x1.0000000000000p+0f, - -0x1.55551a83a9472p-4f, 0x1.331601c4b8ecfp-6f, - -0x1.6890f49eb0acbp-8f, 0x1.8f3a617040a6ap-10f); + float pe = + fputil::polyeval(x2, 0x1.000000p+0f, -0x1.55551ap-4f, 0x1.33160cp-6f, + -0x1.6890f4p-8f, 0x1.8f3a62p-10f); float approx = sqrt_2_delta * pe; return fputil::cast(approx); } From 9420e19f37020f40a7355e2e142c3c565ec4a706 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Mon, 17 Mar 2025 16:04:47 +0000 Subject: [PATCH 10/17] [libc][math][c23] Update for comments. --- libc/src/math/generic/CMakeLists.txt | 1 - libc/src/math/generic/acoshf16.cpp | 58 +++++++++++++++------- libc/test/src/math/smoke/acoshf16_test.cpp | 2 +- 3 files changed, 42 insertions(+), 19 deletions(-) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index fd5ecd72b765f..cc2062af21f71 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3951,7 +3951,6 @@ add_entrypoint_object( HDRS ../acoshf16.h DEPENDS - .common_constants .explogxf libc.hdr.errno_macros libc.hdr.fenv_macros diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index c585c5079e33c..588de147c9cb0 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -7,22 +7,27 @@ //===----------------------------------------------------------------------===// #include "src/math/acoshf16.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/cast.h" #include "src/__support/FPUtil/except_value_utils.h" -#include "src/__support/FPUtil/generic/sqrt.h" #include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/sqrt.h" #include "src/__support/macros/optimization.h" #include "src/math/generic/explogxf.h" namespace LIBC_NAMESPACE_DECL { static constexpr size_t N_EXCEPTS = 2; -static constexpr fputil::ExceptValues ACOSHF16_EXCEPTS{ - {// (input, RZ output, RU offset, RD offset, RN offset) - {0x41B7, 0x3ED8, 1, 0, 0}, - {0x3CE4, 0x393E, 1, 0, 1}}}; +static constexpr fputil::ExceptValues ACOSHF16_EXCEPTS{{ + // (input, RZ output, RU offset, RD offset, RN offset) + // x = 0x1.6ep+1, acoshf16(x) = 0x1.dbp+0 (RZ) + {0x41B7, 0x3ED8, 1, 0, 0}, + // x = 0x1.c8p+0, acoshf16(x) = 0x1.27cp-1 (RZ) + {0x3CE4, 0x393E, 1, 0, 1}, +}}; LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { using FPBits = fputil::FPBits; @@ -58,29 +63,48 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { // acosh(1.0) exactly equals 0.0 if (LIBC_UNLIKELY(x_u == 0x3c00U)) - return float16(0.0f); + return FPBits::zero().get_val(); if (auto r = ACOSHF16_EXCEPTS.lookup(xbits.uintval()); LIBC_UNLIKELY(r.has_value())) return r.value(); - float xf32 = x; + float xf = x; // High precision for inputs very close to 1.0 - if (LIBC_UNLIKELY(xf32 < 1.25f)) { - float delta = xf32 - 1.0f; + // For inputs close to 1 (1 <= x < 1.25), use polynomial approximation: + // + // Step-by-step derivation: + // 1. Let y = acosh(x), thus x = cosh(y). + // + // 2. Rewrite cosh identity using exponential form: + // cosh(y) = (e^y + e^{-y}) / 2 + // For y close to 0, let y = sqrt(2 * delta), thus: + // x = cosh(y) ≈ 1 + delta (since cosh(0) = 1, and delta is small) + // thus delta = x - 1. + // + // 3: Express y in terms of delta (for small delta): + // y ≈ sqrt(2 * delta) + // + // 4: Expand acosh(1 + delta) using a Taylor expansion around delta = 0: + // acosh(1 + delta) ≈ sqrt(2 * delta) * P(delta), where P(delta) + // is a polynomial approximation obtained by fitting the function + // precisely in the interval [0, 0.25]. + // + // Because delta = x - 1 and 0 <= delta < 0.25, the polynomial approximation + // remains numerically stable and accurate in this domain, ensuring high + // precision. + if (LIBC_UNLIKELY(xf < 1.25f)) { + float delta = xf - 1.0f; float sqrt_2_delta = fputil::sqrt(2.0 * delta); - float x2 = delta; - float pe = - fputil::polyeval(x2, 0x1.000000p+0f, -0x1.55551ap-4f, 0x1.33160cp-6f, - -0x1.6890f4p-8f, 0x1.8f3a62p-10f); + float pe = fputil::polyeval(delta, 0x1p+0f, -0x1.55551ap-4f, 0x1.33160cp-6f, + -0x1.6890f4p-8f, 0x1.8f3a62p-10f); float approx = sqrt_2_delta * pe; return fputil::cast(approx); } - // Standard computation for general case. - float sqrt_term = - fputil::sqrt(fputil::multiply_add(xf32, xf32, -1.0f)); - float result = static_cast(log_eval(xf32 + sqrt_term)); + // acosh(x) = log(x + sqrt(x^2 - 1)) + float sqrt_term = fputil::sqrt(fputil::multiply_add(xf, xf, -1.0f)); + float result = static_cast(log_eval(xf + sqrt_term)); return fputil::cast(result); } diff --git a/libc/test/src/math/smoke/acoshf16_test.cpp b/libc/test/src/math/smoke/acoshf16_test.cpp index b8d13b9fe88c3..083dc3aa63b7d 100644 --- a/libc/test/src/math/smoke/acoshf16_test.cpp +++ b/libc/test/src/math/smoke/acoshf16_test.cpp @@ -21,7 +21,7 @@ TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acoshf16(sNaN), FE_INVALID); EXPECT_MATH_ERRNO(0); - EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(zero)); + EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acoshf16(zero), FE_INVALID); EXPECT_MATH_ERRNO(EDOM); EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(neg_zero)); From 573e9d05e4b0be40309f2bc6da3b137fd781d9e0 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Mon, 17 Mar 2025 16:29:55 +0000 Subject: [PATCH 11/17] [libc][math][c23] Update for comments again. --- libc/src/math/generic/acoshf16.cpp | 40 +++++++++++++---------- libc/test/src/math/acoshf16_test.cpp | 47 ---------------------------- 2 files changed, 24 insertions(+), 63 deletions(-) diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index 588de147c9cb0..6f030e46e8532 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -70,29 +70,37 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { return r.value(); float xf = x; - // High precision for inputs very close to 1.0 - // For inputs close to 1 (1 <= x < 1.25), use polynomial approximation: + // High precision polynomial approximation for inputs very close to 1.0 + // Specifically, for inputs within the range [1, 1.25), we employ the + // following step-by-step Taylor expansion derivation to maintain numerical + // accuracy: // // Step-by-step derivation: - // 1. Let y = acosh(x), thus x = cosh(y). + // 1. Define y = acosh(x), thus by definition x = cosh(y). // - // 2. Rewrite cosh identity using exponential form: + // 2. Expand cosh(y) using exponential identities: // cosh(y) = (e^y + e^{-y}) / 2 - // For y close to 0, let y = sqrt(2 * delta), thus: - // x = cosh(y) ≈ 1 + delta (since cosh(0) = 1, and delta is small) - // thus delta = x - 1. + // For small y, let us set y ≈ sqrt(2 * delta), thus: + // x ≈ cosh(y) ≈ 1 + delta, for small delta + // hence delta = x - 1. // - // 3: Express y in terms of delta (for small delta): - // y ≈ sqrt(2 * delta) + // 3. Express y explicitly in terms of delta (for small delta): + // y = acosh(1 + delta) ≈ sqrt(2 * delta) for very small delta. // - // 4: Expand acosh(1 + delta) using a Taylor expansion around delta = 0: - // acosh(1 + delta) ≈ sqrt(2 * delta) * P(delta), where P(delta) - // is a polynomial approximation obtained by fitting the function - // precisely in the interval [0, 0.25]. + // 4. Use Taylor expansion around delta = 0 to obtain a more accurate + // polynomial: + // acosh(1 + delta) ≈ sqrt(2 * delta) * [1 - delta/12 + 3*delta^2/160 - + // 5*delta^3/896 + 35*delta^4/18432 + ...] For practical computation and + // precision, truncate and fit the polynomial precisely in the range [0, + // 0.25]. // - // Because delta = x - 1 and 0 <= delta < 0.25, the polynomial approximation - // remains numerically stable and accurate in this domain, ensuring high - // precision. + // 5. The implemented polynomial approximation (coefficients obtained from + // careful numerical fitting) is: + // P(delta) ≈ 1 - 0x1.55551ap-4 * delta + 0x1.33160cp-6 * delta^2 - + // 0x1.6890f4p-8 * delta^3 + 0x1.8f3a62p-10 * delta^4 + // + // Since delta = x - 1, and 0 <= delta < 0.25, this approximation achieves + // high precision and numerical stability. if (LIBC_UNLIKELY(xf < 1.25f)) { float delta = xf - 1.0f; float sqrt_2_delta = fputil::sqrt(2.0 * delta); diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index a779d3110f329..5df6e0315211e 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -20,56 +20,9 @@ namespace mpfr = LIBC_NAMESPACE::testing::mpfr; static constexpr uint16_t START = 0x3c00U; static constexpr uint16_t STOP = 0x7bffU; -TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { - LIBC_NAMESPACE::libc_errno = 0; - - EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(aNaN)); - EXPECT_MATH_ERRNO(0); - - EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acoshf16(sNaN), FE_INVALID); - EXPECT_MATH_ERRNO(0); - - EXPECT_FP_EQ_ALL_ROUNDING(float16(0.0f), - LIBC_NAMESPACE::acoshf16(float16(1.0f))); - EXPECT_MATH_ERRNO(0); - - EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(float16(0.5f))); - EXPECT_MATH_ERRNO(EDOM); - - EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::acoshf16(inf)); - EXPECT_MATH_ERRNO(0); - - EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::acoshf16(neg_inf)); - EXPECT_MATH_ERRNO(EDOM); -} - TEST_F(LlvmLibcAcoshf16Test, PositiveRange) { for (uint16_t v = START; v <= STOP; ++v) { float16 x = FPBits(v).get_val(); - - EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, x, - LIBC_NAMESPACE::acoshf16(x), 0.5); - } -} - -TEST_F(LlvmLibcAcoshf16Test, SpecificBitPatterns) { - constexpr int N = 12; - constexpr uint16_t INPUTS[N] = { - 0x3C00, // x = 1.0 - 0x3C01, // x = just above 1.0 (minimally larger than 1) - 0x3E00, // x = 1.5 - 0x4200, // x = 3.0 - 0x4500, // x = 5.0 - 0x4900, // x = 10.0 - 0x51FF, // x = ~47.94 - 0x5CB0, // x = ~300.0 - 0x643F, // x = ~1087.6 - 0x77FF, // x = just below next exponent interval (max for exponent 0x1D) - 0x7801, // x = just above previous value (min for exponent 0x1E) - 0x7BFF // x = 65504.0 (max finite half) - }; - for (int i = 0; i < N; ++i) { - float16 x = FPBits(INPUTS[i]).get_val(); EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Acosh, x, LIBC_NAMESPACE::acoshf16(x), 0.5); } From 33c2efa971e12af2da8923b29684a11086e79200 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 18 Mar 2025 14:14:10 +0000 Subject: [PATCH 12/17] [libc][math][c23] Do some change. --- libc/src/math/generic/acoshf16.cpp | 54 ++++++++-------------- libc/test/src/math/acoshf16_test.cpp | 3 +- libc/test/src/math/smoke/CMakeLists.txt | 3 +- libc/test/src/math/smoke/acoshf16_test.cpp | 7 ++- 4 files changed, 26 insertions(+), 41 deletions(-) diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index 6f030e46e8532..c16a0d91ffff9 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -1,4 +1,4 @@ -//===-- Half-precision acoshf16 function ----------------------------------===// +//===-- Half-precision acosh(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. @@ -7,6 +7,7 @@ //===----------------------------------------------------------------------===// #include "src/math/acoshf16.h" +#include "explogxf.h" #include "hdr/errno_macros.h" #include "hdr/fenv_macros.h" #include "src/__support/FPUtil/FEnvImpl.h" @@ -16,16 +17,15 @@ #include "src/__support/FPUtil/multiply_add.h" #include "src/__support/FPUtil/sqrt.h" #include "src/__support/macros/optimization.h" -#include "src/math/generic/explogxf.h" namespace LIBC_NAMESPACE_DECL { static constexpr size_t N_EXCEPTS = 2; static constexpr fputil::ExceptValues ACOSHF16_EXCEPTS{{ // (input, RZ output, RU offset, RD offset, RN offset) - // x = 0x1.6ep+1, acoshf16(x) = 0x1.dbp+0 (RZ) + // x = 0x1.6dcp+1, acoshf16(x) = 0x1.b6p+0 (RZ) {0x41B7, 0x3ED8, 1, 0, 0}, - // x = 0x1.c8p+0, acoshf16(x) = 0x1.27cp-1 (RZ) + // x = 0x1.39p+0, acoshf16(x) = 0x1.4f8p-1 (RZ) {0x3CE4, 0x393E, 1, 0, 1}, }}; @@ -70,42 +70,24 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { return r.value(); float xf = x; - // High precision polynomial approximation for inputs very close to 1.0 - // Specifically, for inputs within the range [1, 1.25), we employ the - // following step-by-step Taylor expansion derivation to maintain numerical - // accuracy: + // High-precision polynomial approximation for inputs close to 1.0 + // ([1, 1.25)). // - // Step-by-step derivation: - // 1. Define y = acosh(x), thus by definition x = cosh(y). - // - // 2. Expand cosh(y) using exponential identities: - // cosh(y) = (e^y + e^{-y}) / 2 - // For small y, let us set y ≈ sqrt(2 * delta), thus: - // x ≈ cosh(y) ≈ 1 + delta, for small delta - // hence delta = x - 1. - // - // 3. Express y explicitly in terms of delta (for small delta): - // y = acosh(1 + delta) ≈ sqrt(2 * delta) for very small delta. - // - // 4. Use Taylor expansion around delta = 0 to obtain a more accurate - // polynomial: - // acosh(1 + delta) ≈ sqrt(2 * delta) * [1 - delta/12 + 3*delta^2/160 - - // 5*delta^3/896 + 35*delta^4/18432 + ...] For practical computation and - // precision, truncate and fit the polynomial precisely in the range [0, - // 0.25]. - // - // 5. The implemented polynomial approximation (coefficients obtained from - // careful numerical fitting) is: - // P(delta) ≈ 1 - 0x1.55551ap-4 * delta + 0x1.33160cp-6 * delta^2 - - // 0x1.6890f4p-8 * delta^3 + 0x1.8f3a62p-10 * delta^4 - // - // Since delta = x - 1, and 0 <= delta < 0.25, this approximation achieves - // high precision and numerical stability. + // Brief derivation: + // 1. Start from the definition: acosh(x) = y means x = cosh(y). + // 2. Set x = 1 + delta. For small delta, y ≈ sqrt(2 * delta). + // 3. Expand acosh(1 + delta) using Taylor series around delta=0: + // acosh(1 + delta) ≈ sqrt(2 * delta) * [1 - delta/12 + 3*delta^2/160 + // - 5*delta^3/896 + 35*delta^4/18432 + ...] + // 4. Truncate the series to fit accurately for delta in [0, 0.25]. + // 5. Polynomial coefficients (from sollya) used here are: + // P(delta) ≈ 1 - 0x1.555556p-4 * delta + 0x1.333334p-6 * delta^2 + // - 0x1.6db6dcp-8 * delta^3 + 0x1.f1c71cp-10 * delta^4 if (LIBC_UNLIKELY(xf < 1.25f)) { float delta = xf - 1.0f; float sqrt_2_delta = fputil::sqrt(2.0 * delta); - float pe = fputil::polyeval(delta, 0x1p+0f, -0x1.55551ap-4f, 0x1.33160cp-6f, - -0x1.6890f4p-8f, 0x1.8f3a62p-10f); + float pe = fputil::polyeval(delta, 0x1p+0f, -0x1.555556p-4f, 0x1.333334p-6f, + -0x1.6db6dcp-8f, 0x1.f1c71cp-10f); float approx = sqrt_2_delta * pe; return fputil::cast(approx); } diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index 5df6e0315211e..cfce364b1846d 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -6,11 +6,10 @@ // //===----------------------------------------------------------------------===// -#include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/cast.h" #include "src/errno/libc_errno.h" #include "src/math/acoshf16.h" #include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" #include "utils/MPFRWrapper/MPFRUtils.h" #include diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index 2c152e39e1c39..13cdbd6c3f8a8 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -3953,7 +3953,8 @@ add_fp_unittest( acoshf16_test.cpp DEPENDS libc.src.errno.errno - libc.src.math.acoshf16 + libc.src.math.acoshf16 + libc.src.__support.FPUtil.cast ) add_fp_unittest( diff --git a/libc/test/src/math/smoke/acoshf16_test.cpp b/libc/test/src/math/smoke/acoshf16_test.cpp index 083dc3aa63b7d..b8615721f329d 100644 --- a/libc/test/src/math/smoke/acoshf16_test.cpp +++ b/libc/test/src/math/smoke/acoshf16_test.cpp @@ -6,6 +6,7 @@ // //===----------------------------------------------------------------------===// +#include "src/__support/FPUtil/cast.h" #include "src/errno/libc_errno.h" #include "src/math/acoshf16.h" #include "test/UnitTest/FPMatcher.h" @@ -33,9 +34,11 @@ TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(neg_inf)); EXPECT_MATH_ERRNO(EDOM); - EXPECT_FP_EQ(float16(0.0f), LIBC_NAMESPACE::acoshf16(float16(1.0f))); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::acoshf16( + LIBC_NAMESPACE::fputil::cast(1.0))); EXPECT_MATH_ERRNO(0); - EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(float16(0.5f))); + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16( + LIBC_NAMESPACE::fputil::cast(0.5))); EXPECT_MATH_ERRNO(EDOM); } From 4a22df94e3eb72ecb722c2048c5fc67e3ed64d33 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 18 Mar 2025 14:28:41 +0000 Subject: [PATCH 13/17] [libc][math][c23] Update test. --- libc/test/src/math/smoke/acoshf16_test.cpp | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/libc/test/src/math/smoke/acoshf16_test.cpp b/libc/test/src/math/smoke/acoshf16_test.cpp index b8615721f329d..4450b605a1226 100644 --- a/libc/test/src/math/smoke/acoshf16_test.cpp +++ b/libc/test/src/math/smoke/acoshf16_test.cpp @@ -28,12 +28,20 @@ TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(neg_zero)); EXPECT_MATH_ERRNO(EDOM); + EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acoshf16(neg_zero), + FE_INVALID); + EXPECT_MATH_ERRNO(EDOM); + EXPECT_FP_EQ(inf, LIBC_NAMESPACE::acoshf16(inf)); EXPECT_MATH_ERRNO(0); EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16(neg_inf)); EXPECT_MATH_ERRNO(EDOM); + EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::acoshf16(neg_inf), + FE_INVALID); + EXPECT_MATH_ERRNO(EDOM); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::acoshf16( LIBC_NAMESPACE::fputil::cast(1.0))); EXPECT_MATH_ERRNO(0); From c8f5517a59b49279f224e14aeba7affc0362e908 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Tue, 18 Mar 2025 16:16:45 +0000 Subject: [PATCH 14/17] [libc][math][c23] Update range. --- libc/test/src/math/acoshf16_test.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/test/src/math/acoshf16_test.cpp b/libc/test/src/math/acoshf16_test.cpp index cfce364b1846d..7348018396bd7 100644 --- a/libc/test/src/math/acoshf16_test.cpp +++ b/libc/test/src/math/acoshf16_test.cpp @@ -17,7 +17,7 @@ using LlvmLibcAcoshf16Test = LIBC_NAMESPACE::testing::FPTest; namespace mpfr = LIBC_NAMESPACE::testing::mpfr; static constexpr uint16_t START = 0x3c00U; -static constexpr uint16_t STOP = 0x7bffU; +static constexpr uint16_t STOP = 0x7c00; TEST_F(LlvmLibcAcoshf16Test, PositiveRange) { for (uint16_t v = START; v <= STOP; ++v) { From ce34bb75ebcd4d6f40d6eecfe87c98c640969542 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Thu, 20 Mar 2025 16:36:39 +0000 Subject: [PATCH 15/17] [libc][math][c23] Do some change and add some tests. --- libc/config/linux/aarch64/entrypoints.txt | 1 + libc/config/linux/riscv/entrypoints.txt | 1 + libc/src/math/generic/acoshf16.cpp | 53 ++++++++++++---------- libc/test/src/math/smoke/acoshf16_test.cpp | 18 ++++++++ 4 files changed, 49 insertions(+), 24 deletions(-) diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt index c1e688ea7e86b..1558209055de6 100644 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -644,6 +644,7 @@ endif() if(LIBC_TYPES_HAS_FLOAT16) list(APPEND TARGET_LIBM_ENTRYPOINTS # math.h C23 _Float16 entrypoints + libc.src.math.acoshf16 libc.src.math.canonicalizef16 libc.src.math.ceilf16 libc.src.math.copysignf16 diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt index ff238c3c32d98..7246da7114471 100644 --- a/libc/config/linux/riscv/entrypoints.txt +++ b/libc/config/linux/riscv/entrypoints.txt @@ -397,6 +397,7 @@ set(TARGET_LIBM_ENTRYPOINTS # math.h entrypoints libc.src.math.acosf libc.src.math.acoshf + libc.src.math.acoshf16 libc.src.math.asinf libc.src.math.asinhf libc.src.math.atan2 diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index c16a0d91ffff9..a51cb37f5ed64 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -33,19 +33,22 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { using FPBits = fputil::FPBits; FPBits xbits(x); uint16_t x_u = xbits.uintval(); - uint16_t x_abs = x_u & 0x7fff; + + // Domain error for inputs less than 1.0. + if (LIBC_UNLIKELY(x <= 1.0f)) { + if (x == 1.0f) + return FPBits::zero().get_val(); + fputil::set_errno_if_required(EDOM); + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } // Check for NaN input first. - if (LIBC_UNLIKELY(xbits.is_nan())) { + if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) { if (xbits.is_signaling_nan()) { fputil::raise_except_if_required(FE_INVALID); return FPBits::quiet_nan().get_val(); } - return x; - } - - // Check for infinite inputs. - if (LIBC_UNLIKELY(xbits.is_inf())) { if (xbits.is_neg()) { fputil::set_errno_if_required(EDOM); fputil::raise_except_if_required(FE_INVALID); @@ -54,17 +57,6 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { return x; } - // Domain error for inputs less than 1.0. - if (LIBC_UNLIKELY(x_abs < 0x3c00U)) { - fputil::set_errno_if_required(EDOM); - fputil::raise_except_if_required(FE_INVALID); - return FPBits::quiet_nan().get_val(); - } - - // acosh(1.0) exactly equals 0.0 - if (LIBC_UNLIKELY(x_u == 0x3c00U)) - return FPBits::zero().get_val(); - if (auto r = ACOSHF16_EXCEPTS.lookup(xbits.uintval()); LIBC_UNLIKELY(r.has_value())) return r.value(); @@ -74,16 +66,29 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { // ([1, 1.25)). // // Brief derivation: - // 1. Start from the definition: acosh(x) = y means x = cosh(y). - // 2. Set x = 1 + delta. For small delta, y ≈ sqrt(2 * delta). - // 3. Expand acosh(1 + delta) using Taylor series around delta=0: + // 1. Expand acosh(1 + delta) using Taylor series around delta=0: // acosh(1 + delta) ≈ sqrt(2 * delta) * [1 - delta/12 + 3*delta^2/160 // - 5*delta^3/896 + 35*delta^4/18432 + ...] - // 4. Truncate the series to fit accurately for delta in [0, 0.25]. - // 5. Polynomial coefficients (from sollya) used here are: + // 2. Truncate the series to fit accurately for delta in [0, 0.25]. + // 3. Polynomial coefficients (from sollya) used here are: // P(delta) ≈ 1 - 0x1.555556p-4 * delta + 0x1.333334p-6 * delta^2 // - 0x1.6db6dcp-8 * delta^3 + 0x1.f1c71cp-10 * delta^4 - if (LIBC_UNLIKELY(xf < 1.25f)) { + // 4. The Sollya commands used to generate these coefficients were: + // > display = hexadecimal; + // > round(1/12, SG, RN); + // > round(3/160, SG, RN); + // > round(5/896, SG, RN); + // > round(35/18432, SG, RN); + // With hexadecimal display mode enabled, the outputs were: + // 0x1.555556p-4 + // 0x1.333334p-6 + // 0x1.6db6dcp-8 + // 0x1.f1c71cp-10 + // 5. The maximum absolute error, estimated using: + // dirtyinfnorm(acosh(1 + x) - sqrt(2*x) * P(x), [0, 0.25]) + // is: + // 0x1.d84281p-22 + if (LIBC_UNLIKELY(x_u < 1.25f)) { float delta = xf - 1.0f; float sqrt_2_delta = fputil::sqrt(2.0 * delta); float pe = fputil::polyeval(delta, 0x1p+0f, -0x1.555556p-4f, 0x1.333334p-6f, diff --git a/libc/test/src/math/smoke/acoshf16_test.cpp b/libc/test/src/math/smoke/acoshf16_test.cpp index 4450b605a1226..7681c2a4e7fbc 100644 --- a/libc/test/src/math/smoke/acoshf16_test.cpp +++ b/libc/test/src/math/smoke/acoshf16_test.cpp @@ -49,4 +49,22 @@ TEST_F(LlvmLibcAcoshf16Test, SpecialNumbers) { EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::acoshf16( LIBC_NAMESPACE::fputil::cast(0.5))); EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ_WITH_EXCEPTION( + aNaN, + LIBC_NAMESPACE::acoshf16(LIBC_NAMESPACE::fputil::cast(-1.0)), + FE_INVALID); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ_WITH_EXCEPTION( + aNaN, + LIBC_NAMESPACE::acoshf16(LIBC_NAMESPACE::fputil::cast(-2.0)), + FE_INVALID); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ_WITH_EXCEPTION( + aNaN, + LIBC_NAMESPACE::acoshf16(LIBC_NAMESPACE::fputil::cast(-3.0)), + FE_INVALID); + EXPECT_MATH_ERRNO(EDOM); } From 45386afc707333bfc9298416158fc50bc1ae8ef1 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Sat, 22 Mar 2025 18:08:13 +0800 Subject: [PATCH 16/17] [libc][math][c23] Update for comments. --- libc/src/math/generic/acoshf16.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index a51cb37f5ed64..4f169ca7c3bf2 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -88,7 +88,7 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { // dirtyinfnorm(acosh(1 + x) - sqrt(2*x) * P(x), [0, 0.25]) // is: // 0x1.d84281p-22 - if (LIBC_UNLIKELY(x_u < 1.25f)) { + if (LIBC_UNLIKELY(x_u < 0x3D00U)) { float delta = xf - 1.0f; float sqrt_2_delta = fputil::sqrt(2.0 * delta); float pe = fputil::polyeval(delta, 0x1p+0f, -0x1.555556p-4f, 0x1.333334p-6f, From 73b34db959a0310e0e9a5326e2a798256f45dfa9 Mon Sep 17 00:00:00 2001 From: Harrison Hao Date: Sat, 22 Mar 2025 16:02:22 +0000 Subject: [PATCH 17/17] [libc][math][c23] Update for comments. --- libc/config/linux/riscv/entrypoints.txt | 1 - libc/src/math/generic/CMakeLists.txt | 1 + libc/src/math/generic/acoshf16.cpp | 21 ++++++++++++--------- 3 files changed, 13 insertions(+), 10 deletions(-) diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt index 7246da7114471..ff238c3c32d98 100644 --- a/libc/config/linux/riscv/entrypoints.txt +++ b/libc/config/linux/riscv/entrypoints.txt @@ -397,7 +397,6 @@ set(TARGET_LIBM_ENTRYPOINTS # math.h entrypoints libc.src.math.acosf libc.src.math.acoshf - libc.src.math.acoshf16 libc.src.math.asinf libc.src.math.asinhf libc.src.math.atan2 diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index cc2062af21f71..d9650f9fd6b35 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3955,6 +3955,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/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index 4f169ca7c3bf2..44783a8749ac2 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -12,10 +12,13 @@ #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" #include "src/__support/FPUtil/except_value_utils.h" #include "src/__support/FPUtil/multiply_add.h" #include "src/__support/FPUtil/sqrt.h" +#include "src/__support/common.h" +#include "src/__support/macros/config.h" #include "src/__support/macros/optimization.h" namespace LIBC_NAMESPACE_DECL { @@ -34,15 +37,6 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { FPBits xbits(x); uint16_t x_u = xbits.uintval(); - // Domain error for inputs less than 1.0. - if (LIBC_UNLIKELY(x <= 1.0f)) { - if (x == 1.0f) - return FPBits::zero().get_val(); - fputil::set_errno_if_required(EDOM); - fputil::raise_except_if_required(FE_INVALID); - return FPBits::quiet_nan().get_val(); - } - // Check for NaN input first. if (LIBC_UNLIKELY(xbits.is_inf_or_nan())) { if (xbits.is_signaling_nan()) { @@ -57,6 +51,15 @@ LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { return x; } + // Domain error for inputs less than 1.0. + if (LIBC_UNLIKELY(x <= 1.0f)) { + if (x == 1.0f) + return FPBits::zero().get_val(); + fputil::set_errno_if_required(EDOM); + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + if (auto r = ACOSHF16_EXCEPTS.lookup(xbits.uintval()); LIBC_UNLIKELY(r.has_value())) return r.value();