diff --git a/libc/shared/math.h b/libc/shared/math.h index 663900b1d3519..1c3b46e439c1d 100644 --- a/libc/shared/math.h +++ b/libc/shared/math.h @@ -56,6 +56,7 @@ #include "math/expf16.h" #include "math/expm1.h" #include "math/expm1f.h" +#include "math/expm1f16.h" #include "math/frexpf.h" #include "math/frexpf128.h" #include "math/frexpf16.h" diff --git a/libc/shared/math/expm1f16.h b/libc/shared/math/expm1f16.h new file mode 100644 index 0000000000000..5698400d7066a --- /dev/null +++ b/libc/shared/math/expm1f16.h @@ -0,0 +1,29 @@ +//===-- Shared expm1f16 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_EXPM1F16_H +#define LLVM_LIBC_SHARED_MATH_EXPM1F16_H + +#include "include/llvm-libc-macros/float16-macros.h" +#include "shared/libc_common.h" + +#ifdef LIBC_TYPES_HAS_FLOAT16 + +#include "src/__support/math/expm1f16.h" + +namespace LIBC_NAMESPACE_DECL { +namespace shared { + +using math::expm1f16; + +} // namespace shared +} // namespace LIBC_NAMESPACE_DECL + +#endif // LIBC_TYPES_HAS_FLOAT16 + +#endif // LLVM_LIBC_SHARED_MATH_EXPM1F16_H diff --git a/libc/src/__support/math/CMakeLists.txt b/libc/src/__support/math/CMakeLists.txt index 00b44c5725565..7fe029ddcfb65 100644 --- a/libc/src/__support/math/CMakeLists.txt +++ b/libc/src/__support/math/CMakeLists.txt @@ -923,6 +923,22 @@ add_header_library( libc.src.errno.errno ) +add_header_library( + expm1f16 + HDRS + expm1f16.h + DEPENDS + .expxf16_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.rounding_mode + libc.src.__support.macros.optimization +) + add_header_library( range_reduction_double HDRS diff --git a/libc/src/__support/math/expm1f16.h b/libc/src/__support/math/expm1f16.h new file mode 100644 index 0000000000000..79547b62b0892 --- /dev/null +++ b/libc/src/__support/math/expm1f16.h @@ -0,0 +1,153 @@ +//===-- Implementation header for expm1f16 ----------------------*- 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_EXPM1F16_H +#define LLVM_LIBC_SRC___SUPPORT_MATH_EXPM1F16_H + +#include "include/llvm-libc-macros/float16-macros.h" + +#ifdef LIBC_TYPES_HAS_FLOAT16 + +#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/rounding_mode.h" +#include "src/__support/common.h" +#include "src/__support/macros/config.h" +#include "src/__support/macros/optimization.h" +#include "src/__support/math/expxf16_utils.h" + +namespace LIBC_NAMESPACE_DECL { + +namespace math { + +LIBC_INLINE static constexpr float16 expm1f16(float16 x) { +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + constexpr fputil::ExceptValues EXPM1F16_EXCEPTS_LO = {{ + // (input, RZ output, RU offset, RD offset, RN offset) + // x = 0x1.564p-5, expm1f16(x) = 0x1.5d4p-5 (RZ) + {0x2959U, 0x2975U, 1U, 0U, 1U}, + }}; + +#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT + constexpr size_t N_EXPM1F16_EXCEPTS_HI = 2; +#else + constexpr size_t N_EXPM1F16_EXCEPTS_HI = 3; +#endif + + constexpr fputil::ExceptValues + EXPM1F16_EXCEPTS_HI = {{ + // (input, RZ output, RU offset, RD offset, RN offset) + // x = 0x1.c34p+0, expm1f16(x) = 0x1.34cp+2 (RZ) + {0x3f0dU, 0x44d3U, 1U, 0U, 1U}, + // x = -0x1.e28p-3, expm1f16(x) = -0x1.adcp-3 (RZ) + {0xb38aU, 0xb2b7U, 0U, 1U, 1U}, +#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT + // x = 0x1.a08p-3, exp10m1f(x) = 0x1.cdcp-3 (RZ) + {0x3282U, 0x3337U, 1U, 0U, 0U}, +#endif + }}; +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + using namespace math::expxf16_internal; + using FPBits = fputil::FPBits; + FPBits x_bits(x); + + uint16_t x_u = x_bits.uintval(); + uint16_t x_abs = x_u & 0x7fffU; + + // When |x| <= 2^(-3), or |x| >= -11 * log(2), or x is NaN. + if (LIBC_UNLIKELY(x_abs <= 0x3000U || x_abs >= 0x47a0U)) { + // expm1(NaN) = NaN + if (x_bits.is_nan()) { + if (x_bits.is_signaling_nan()) { + fputil::raise_except_if_required(FE_INVALID); + return FPBits::quiet_nan().get_val(); + } + + return x; + } + + // expm1(+/-0) = +/-0 + if (x_abs == 0) + return x; + + // When x >= 16 * log(2). + if (x_bits.is_pos() && x_abs >= 0x498cU) { + // expm1(+inf) = +inf + if (x_bits.is_inf()) + return FPBits::inf().get_val(); + + switch (fputil::quick_get_round()) { + case FE_TONEAREST: + case FE_UPWARD: + fputil::set_errno_if_required(ERANGE); + fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT); + return FPBits::inf().get_val(); + default: + return FPBits::max_normal().get_val(); + } + } + + // When x <= -11 * log(2). + if (x_u >= 0xc7a0U) { + // expm1(-inf) = -1 + if (x_bits.is_inf()) + return FPBits::one(Sign::NEG).get_val(); + + // When x > -0x1.0ap+3, round(expm1(x), HP, RN) = -1. + if (x_u > 0xc828U) + return fputil::round_result_slightly_up( + FPBits::one(Sign::NEG).get_val()); + // When x <= -0x1.0ap+3, round(expm1(x), HP, RN) = -0x1.ffcp-1. + return fputil::round_result_slightly_down( + fputil::cast(-0x1.ffcp-1)); + } + + // When 0 < |x| <= 2^(-3). + if (x_abs <= 0x3000U && !x_bits.is_zero()) { + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + if (auto r = EXPM1F16_EXCEPTS_LO.lookup(x_u); + LIBC_UNLIKELY(r.has_value())) + return r.value(); +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + float xf = x; + // Degree-5 minimax polynomial generated by Sollya with the following + // commands: + // > display = hexadecimal; + // > P = fpminimax(expm1(x)/x, 4, [|SG...|], [-2^-3, 2^-3]); + // > x * P; + return fputil::cast( + xf * fputil::polyeval(xf, 0x1p+0f, 0x1.fffff8p-2f, 0x1.555556p-3f, + 0x1.55905ep-5f, 0x1.1124c2p-7f)); + } + } + +#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS + if (auto r = EXPM1F16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value())) + return r.value(); +#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS + + // exp(x) = exp(hi + mid) * exp(lo) + auto [exp_hi_mid, exp_lo] = exp_range_reduction(x); + // expm1(x) = exp(hi + mid) * exp(lo) - 1 + return fputil::cast(fputil::multiply_add(exp_hi_mid, exp_lo, -1.0f)); +} + +} // namespace math + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LIBC_TYPES_HAS_FLOAT16 + +#endif // LLVM_LIBC_SRC___SUPPORT_MATH_EXPM1F16_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 2a2f4e7567c71..c14f356b39476 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -1571,17 +1571,7 @@ add_entrypoint_object( HDRS ../expm1f16.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.rounding_mode - libc.src.__support.macros.optimization - libc.src.__support.math.expxf16_utils + libc.src.__support.math.expm1f16 ) add_entrypoint_object( diff --git a/libc/src/math/generic/expm1f16.cpp b/libc/src/math/generic/expm1f16.cpp index c2231f0aca715..68bf21df1721e 100644 --- a/libc/src/math/generic/expm1f16.cpp +++ b/libc/src/math/generic/expm1f16.cpp @@ -7,135 +7,9 @@ //===----------------------------------------------------------------------===// #include "src/math/expm1f16.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/rounding_mode.h" -#include "src/__support/common.h" -#include "src/__support/macros/config.h" -#include "src/__support/macros/optimization.h" -#include "src/__support/math/expxf16_utils.h" +#include "src/__support/math/expm1f16.h" namespace LIBC_NAMESPACE_DECL { - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS -static constexpr fputil::ExceptValues EXPM1F16_EXCEPTS_LO = {{ - // (input, RZ output, RU offset, RD offset, RN offset) - // x = 0x1.564p-5, expm1f16(x) = 0x1.5d4p-5 (RZ) - {0x2959U, 0x2975U, 1U, 0U, 1U}, -}}; - -#ifdef LIBC_TARGET_CPU_HAS_FMA_FLOAT -static constexpr size_t N_EXPM1F16_EXCEPTS_HI = 2; -#else -static constexpr size_t N_EXPM1F16_EXCEPTS_HI = 3; -#endif - -static constexpr fputil::ExceptValues - EXPM1F16_EXCEPTS_HI = {{ - // (input, RZ output, RU offset, RD offset, RN offset) - // x = 0x1.c34p+0, expm1f16(x) = 0x1.34cp+2 (RZ) - {0x3f0dU, 0x44d3U, 1U, 0U, 1U}, - // x = -0x1.e28p-3, expm1f16(x) = -0x1.adcp-3 (RZ) - {0xb38aU, 0xb2b7U, 0U, 1U, 1U}, -#ifndef LIBC_TARGET_CPU_HAS_FMA_FLOAT - // x = 0x1.a08p-3, exp10m1f(x) = 0x1.cdcp-3 (RZ) - {0x3282U, 0x3337U, 1U, 0U, 0U}, -#endif - }}; -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - -LLVM_LIBC_FUNCTION(float16, expm1f16, (float16 x)) { - using namespace math::expxf16_internal; - using FPBits = fputil::FPBits; - FPBits x_bits(x); - - uint16_t x_u = x_bits.uintval(); - uint16_t x_abs = x_u & 0x7fffU; - - // When |x| <= 2^(-3), or |x| >= -11 * log(2), or x is NaN. - if (LIBC_UNLIKELY(x_abs <= 0x3000U || x_abs >= 0x47a0U)) { - // expm1(NaN) = NaN - if (x_bits.is_nan()) { - if (x_bits.is_signaling_nan()) { - fputil::raise_except_if_required(FE_INVALID); - return FPBits::quiet_nan().get_val(); - } - - return x; - } - - // expm1(+/-0) = +/-0 - if (x_abs == 0) - return x; - - // When x >= 16 * log(2). - if (x_bits.is_pos() && x_abs >= 0x498cU) { - // expm1(+inf) = +inf - if (x_bits.is_inf()) - return FPBits::inf().get_val(); - - switch (fputil::quick_get_round()) { - case FE_TONEAREST: - case FE_UPWARD: - fputil::set_errno_if_required(ERANGE); - fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT); - return FPBits::inf().get_val(); - default: - return FPBits::max_normal().get_val(); - } - } - - // When x <= -11 * log(2). - if (x_u >= 0xc7a0U) { - // expm1(-inf) = -1 - if (x_bits.is_inf()) - return FPBits::one(Sign::NEG).get_val(); - - // When x > -0x1.0ap+3, round(expm1(x), HP, RN) = -1. - if (x_u > 0xc828U) - return fputil::round_result_slightly_up( - FPBits::one(Sign::NEG).get_val()); - // When x <= -0x1.0ap+3, round(expm1(x), HP, RN) = -0x1.ffcp-1. - return fputil::round_result_slightly_down( - fputil::cast(-0x1.ffcp-1)); - } - - // When 0 < |x| <= 2^(-3). - if (x_abs <= 0x3000U && !x_bits.is_zero()) { - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - if (auto r = EXPM1F16_EXCEPTS_LO.lookup(x_u); - LIBC_UNLIKELY(r.has_value())) - return r.value(); -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - - float xf = x; - // Degree-5 minimax polynomial generated by Sollya with the following - // commands: - // > display = hexadecimal; - // > P = fpminimax(expm1(x)/x, 4, [|SG...|], [-2^-3, 2^-3]); - // > x * P; - return fputil::cast( - xf * fputil::polyeval(xf, 0x1p+0f, 0x1.fffff8p-2f, 0x1.555556p-3f, - 0x1.55905ep-5f, 0x1.1124c2p-7f)); - } - } - -#ifndef LIBC_MATH_HAS_SKIP_ACCURATE_PASS - if (auto r = EXPM1F16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value())) - return r.value(); -#endif // !LIBC_MATH_HAS_SKIP_ACCURATE_PASS - - // exp(x) = exp(hi + mid) * exp(lo) - auto [exp_hi_mid, exp_lo] = exp_range_reduction(x); - // expm1(x) = exp(hi + mid) * exp(lo) - 1 - return fputil::cast(fputil::multiply_add(exp_hi_mid, exp_lo, -1.0f)); -} +LLVM_LIBC_FUNCTION(float16, expm1f16, (float16 x)) { return math::expm1f16(x); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/shared/CMakeLists.txt b/libc/test/shared/CMakeLists.txt index 64d000d6fc445..aadbd4ad980f0 100644 --- a/libc/test/shared/CMakeLists.txt +++ b/libc/test/shared/CMakeLists.txt @@ -47,6 +47,7 @@ add_fp_unittest( libc.src.__support.math.exp2m1f16 libc.src.__support.math.expm1 libc.src.__support.math.expm1f + libc.src.__support.math.expm1f16 libc.src.__support.math.exp10 libc.src.__support.math.exp10f libc.src.__support.math.exp10f16 diff --git a/libc/test/shared/shared_math_test.cpp b/libc/test/shared/shared_math_test.cpp index bbe1fde87592b..cf75672774105 100644 --- a/libc/test/shared/shared_math_test.cpp +++ b/libc/test/shared/shared_math_test.cpp @@ -31,6 +31,7 @@ TEST(LlvmLibcSharedMathTest, AllFloat16) { EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::exp2f16(0.0f16)); EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::exp2m1f16(0.0f16)); EXPECT_FP_EQ(0x1p+0f16, LIBC_NAMESPACE::shared::expf16(0.0f16)); + EXPECT_FP_EQ(0x0p+0f16, LIBC_NAMESPACE::shared::expm1f16(0.0f16)); ASSERT_FP_EQ(float16(8 << 5), LIBC_NAMESPACE::shared::ldexpf16(8.0f16, 5)); ASSERT_FP_EQ(float16(-1 * (8 << 5)), diff --git a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel index 7c563ff109fb8..3ccddbacf833d 100644 --- a/utils/bazel/llvm-project-overlay/libc/BUILD.bazel +++ b/utils/bazel/llvm-project-overlay/libc/BUILD.bazel @@ -3053,6 +3053,22 @@ libc_support_library( ], ) +libc_support_library( + name = "__support_math_expm1f16", + hdrs = ["src/__support/math/expm1f16.h"], + deps = [ + ":__support_fputil_fma", + ":__support_fputil_multiply_add", + ":__support_fputil_nearest_integer", + ":__support_fputil_polyeval", + ":__support_fputil_rounding_mode", + ":__support_fputil_except_value_utils", + ":__support_macros_optimization", + ":__support_macros_properties_cpu_features", + ":__support_math_expxf16_utils" + ], +) + libc_support_library( name = "__support_range_reduction_double", hdrs = [ @@ -3799,7 +3815,7 @@ libc_math_function( libc_math_function( name = "expm1f16", additional_deps = [ - ":__support_math_expxf16_utils", + ":__support_math_expm1f16", ], )