From 02710692541fe641035aaa158e11964d618ddf71 Mon Sep 17 00:00:00 2001 From: bassiounix Date: Mon, 21 Jul 2025 23:30:07 +0300 Subject: [PATCH] [libc][math] Refactor acoshf16 implementation to header-only in src/__support/math folder. --- libc/shared/math.h | 1 + libc/shared/math/acoshf16.h | 29 +++++ libc/src/__support/math/CMakeLists.txt | 16 +++ libc/src/__support/math/acoshf16.h | 123 ++++++++++++++++++ libc/src/math/generic/CMakeLists.txt | 14 +- libc/src/math/generic/acoshf16.cpp | 99 +------------- libc/test/shared/CMakeLists.txt | 1 + libc/test/shared/shared_math_test.cpp | 2 + .../llvm-project-overlay/libc/BUILD.bazel | 24 ++++ 9 files changed, 200 insertions(+), 109 deletions(-) create mode 100644 libc/shared/math/acoshf16.h create mode 100644 libc/src/__support/math/acoshf16.h diff --git a/libc/shared/math.h b/libc/shared/math.h index 33fe235c51824..e3c674c27ffaf 100644 --- a/libc/shared/math.h +++ b/libc/shared/math.h @@ -15,6 +15,7 @@ #include "math/acosf.h" #include "math/acosf16.h" #include "math/acoshf.h" +#include "math/acoshf16.h" #include "math/erff.h" #include "math/exp.h" #include "math/exp10.h" diff --git a/libc/shared/math/acoshf16.h b/libc/shared/math/acoshf16.h new file mode 100644 index 0000000000000..2f0bc6e80ab6d --- /dev/null +++ b/libc/shared/math/acoshf16.h @@ -0,0 +1,29 @@ +//===-- Shared acoshf16 function --------------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SHARED_MATH_ACOSHF16_H +#define LLVM_LIBC_SHARED_MATH_ACOSHF16_H + +#include "include/llvm-libc-macros/float16-macros.h" +#include "shared/libc_common.h" + +#ifdef LIBC_TYPES_HAS_FLOAT16 + +#include "src/__support/math/acoshf16.h" + +namespace LIBC_NAMESPACE_DECL { +namespace shared { + +using math::acoshf16; + +} // namespace shared +} // namespace LIBC_NAMESPACE_DECL + +#endif // LIBC_TYPES_HAS_FLOAT16 + +#endif // LLVM_LIBC_SHARED_MATH_ACOSHF16_H diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt index b0e562fb20d0d..9a8a4d16a4a0e 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -79,6 +79,22 @@ add_header_library( libc.src.__support.macros.optimization ) +add_header_library( + acoshf16 + HDRS + acoshf16.h + DEPENDS + .acoshf_utils + 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 +) + add_header_library( asin_utils HDRS diff --git a/libc/src/__support/math/acoshf16.h b/libc/src/__support/math/acoshf16.h new file mode 100644 index 0000000000000..15e7f6ae7e208 --- /dev/null +++ b/libc/src/__support/math/acoshf16.h @@ -0,0 +1,123 @@ +//===-- 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___SUPPORT_MATH_ACOSHF_H +#define LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_H + +#include "include/llvm-libc-macros/float16-macros.h" + +#ifdef LIBC_TYPES_HAS_FLOAT16 + +#include "acoshf_utils.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/macros/config.h" +#include "src/__support/macros/optimization.h" + +namespace LIBC_NAMESPACE_DECL { + +namespace math { + +static constexpr float16 acoshf16(float16 x) { + + using namespace acoshf_internal; + constexpr size_t N_EXCEPTS = 2; + constexpr fputil::ExceptValues ACOSHF16_EXCEPTS{{ + // (input, RZ output, RU offset, RD offset, RN offset) + // x = 0x1.6dcp+1, acoshf16(x) = 0x1.b6p+0 (RZ) + {0x41B7, 0x3ED8, 1, 0, 0}, + // x = 0x1.39p+0, acoshf16(x) = 0x1.4f8p-1 (RZ) + {0x3CE4, 0x393E, 1, 0, 1}, + }}; + + using FPBits = fputil::FPBits; + FPBits xbits(x); + uint16_t x_u = xbits.uintval(); + + // Check for NaN input first. + 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(); + } + if (xbits.is_neg()) { + fputil::set_errno_if_required(EDOM); + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + 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(); + + float xf = x; + // High-precision polynomial approximation for inputs close to 1.0 + // ([1, 1.25)). + // + // Brief derivation: + // 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 + ...] + // 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 + // 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 < 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, + -0x1.6db6dcp-8f, 0x1.f1c71cp-10f); + float approx = sqrt_2_delta * pe; + return fputil::cast(approx); + } + + // 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); +} + +} // namespace math + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LIBC_TYPES_HAS_FLOAT16 + +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_ACOSHF_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 86baf2b981bb2..408f99ef30760 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -3878,18 +3878,8 @@ add_entrypoint_object( HDRS ../acoshf16.h DEPENDS - .explogxf - 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 + libc.src.__support.math.acoshf16 + libc.src.errno.errno ) add_entrypoint_object( diff --git a/libc/src/math/generic/acoshf16.cpp b/libc/src/math/generic/acoshf16.cpp index 96ca956c82147..bb3a91f707080 100644 --- a/libc/src/math/generic/acoshf16.cpp +++ b/libc/src/math/generic/acoshf16.cpp @@ -7,105 +7,10 @@ //===----------------------------------------------------------------------===// #include "src/math/acoshf16.h" -#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" -#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" +#include "src/__support/math/acoshf16.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.6dcp+1, acoshf16(x) = 0x1.b6p+0 (RZ) - {0x41B7, 0x3ED8, 1, 0, 0}, - // x = 0x1.39p+0, acoshf16(x) = 0x1.4f8p-1 (RZ) - {0x3CE4, 0x393E, 1, 0, 1}, -}}; - -LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { - using namespace acoshf_internal; - using FPBits = fputil::FPBits; - FPBits xbits(x); - uint16_t x_u = xbits.uintval(); - - // Check for NaN input first. - 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(); - } - if (xbits.is_neg()) { - fputil::set_errno_if_required(EDOM); - fputil::raise_except_if_required(FE_INVALID); - return FPBits::quiet_nan().get_val(); - } - 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(); - - float xf = x; - // High-precision polynomial approximation for inputs close to 1.0 - // ([1, 1.25)). - // - // Brief derivation: - // 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 + ...] - // 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 - // 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 < 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, - -0x1.6db6dcp-8f, 0x1.f1c71cp-10f); - float approx = sqrt_2_delta * pe; - return fputil::cast(approx); - } - - // 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); -} +LLVM_LIBC_FUNCTION(float16, acoshf16, (float16 x)) { return math::acoshf16(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt index a402d6e39fa0b..89b607d7e5cc3 100644 --- a/libc/test/shared/CMakeLists.txt +++ b/libc/test/shared/CMakeLists.txt @@ -11,6 +11,7 @@ add_fp_unittest( libc.src.__support.math.acosf libc.src.__support.math.acosf16 libc.src.__support.math.acoshf + libc.src.__support.math.acoshf16 libc.src.__support.math.erff libc.src.__support.math.exp libc.src.__support.math.exp10 diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp index 25555b8caf534..8d3cebdf0745c 100644 --- a/libc/test/shared/shared_math_test.cpp +++ b/libc/test/shared/shared_math_test.cpp @@ -14,6 +14,8 @@ TEST(LlvmLibcSharedMathTest, AllFloat16) { int exponent; + EXPECT_FP_EQ(0x0p+0f, LIBC_NAMESPACE::shared::acoshf16(1.0f)); + EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::exp10f16(0.0f16)); EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::expf16(0.0f16)); diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel index 053d15869246b..88ccd05b54c1c 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -2141,6 +2141,22 @@ libc_support_library( ], ) +libc_support_library( + name = "__support_math_acoshf16", + hdrs = ["src/__support/math/acoshf16.h"], + deps = [ + ":__support_math_acoshf_utils", + ":__support_fputil_cast", + ":__support_fputil_except_value_utils", + ":__support_fputil_fenv_impl", + ":__support_fputil_fp_bits", + ":__support_fputil_multiply_add", + ":__support_fputil_polyeval", + ":__support_fputil_sqrt", + ":__support_macros_optimization", + ], +) + libc_support_library( name = "__support_math_asin_utils", hdrs = ["src/__support/math/asin_utils.h"], @@ -2682,6 +2698,14 @@ libc_math_function( ], ) +libc_math_function( + name = "acoshf16", + additional_deps = [ + ":__support_math_acoshf16", + ":errno", + ], +) + libc_math_function( name = "asinf", additional_deps = [