diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt index 2334fed773702..bbf49b7ad0e61 100644 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -534,6 +534,7 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.canonicalizef16 libc.src.math.ceilf16 libc.src.math.copysignf16 + libc.src.math.expf16 libc.src.math.f16add libc.src.math.f16addf libc.src.math.f16div diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index 035ceb8ca57bf..ac92f4cd567cb 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -562,6 +562,7 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.canonicalizef16 libc.src.math.ceilf16 libc.src.math.copysignf16 + libc.src.math.expf16 libc.src.math.f16add libc.src.math.f16addf libc.src.math.f16addl diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst index 64de548a1ea1c..b047abc8ba858 100644 --- a/libc/docs/math/index.rst +++ b/libc/docs/math/index.rst @@ -284,7 +284,7 @@ Higher Math Functions +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | erfc | | | | | | 7.12.8.2 | F.10.5.2 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| exp | |check| | |check| | | | | 7.12.6.1 | F.10.3.1 | +| exp | |check| | |check| | | |check| | | 7.12.6.1 | F.10.3.1 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | exp10 | |check| | |check| | | | | 7.12.6.2 | F.10.3.2 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td index d5a5cb6fedb4b..5f8490fe32ca6 100644 --- a/libc/spec/stdc.td +++ b/libc/spec/stdc.td @@ -573,6 +573,7 @@ def StdC : StandardSpec<"stdc"> { FunctionSpec<"exp", RetValSpec, [ArgSpec]>, FunctionSpec<"expf", RetValSpec, [ArgSpec]>, + GuardedFunctionSpec<"expf16", RetValSpec, [ArgSpec], "LIBC_TYPES_HAS_FLOAT16">, FunctionSpec<"exp2", RetValSpec, [ArgSpec]>, FunctionSpec<"exp2f", RetValSpec, [ArgSpec]>, diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index d70af33522d2b..323ce988f8613 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -97,6 +97,7 @@ add_math_entrypoint_object(erff) add_math_entrypoint_object(exp) add_math_entrypoint_object(expf) +add_math_entrypoint_object(expf16) add_math_entrypoint_object(exp2) add_math_entrypoint_object(exp2f) diff --git a/libc/src/math/expf16.h b/libc/src/math/expf16.h new file mode 100644 index 0000000000000..8547f65f942d7 --- /dev/null +++ b/libc/src/math/expf16.h @@ -0,0 +1,21 @@ +//===-- Implementation header for expf16 ------------------------*- 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_EXPF16_H +#define LLVM_LIBC_SRC_MATH_EXPF16_H + +#include "src/__support/macros/config.h" +#include "src/__support/macros/properties/types.h" + +namespace LIBC_NAMESPACE_DECL { + +float16 expf16(float16 x); + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC_MATH_EXPF16_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 80e862542e29c..73c19eeaa0094 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -1226,6 +1226,27 @@ add_entrypoint_object( -O3 ) +add_entrypoint_object( + expf16 + SRCS + expf16.cpp + HDRS + ../expf16.h + DEPENDS + libc.hdr.errno_macros + libc.hdr.fenv_macros + libc.src.__support.CPP.array + libc.src.__support.FPUtil.fenv_impl + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.polyeval + libc.src.__support.FPUtil.multiply_add + libc.src.__support.FPUtil.nearest_integer + libc.src.__support.FPUtil.rounding_mode + libc.src.__support.macros.optimization + COMPILE_OPTIONS + -O3 +) + add_entrypoint_object( exp2 SRCS diff --git a/libc/src/math/generic/expf16.cpp b/libc/src/math/generic/expf16.cpp new file mode 100644 index 0000000000000..b618edc36a046 --- /dev/null +++ b/libc/src/math/generic/expf16.cpp @@ -0,0 +1,173 @@ +//===-- Half-precision e^x function ---------------------------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "src/math/expf16.h" +#include "hdr/errno_macros.h" +#include "hdr/fenv_macros.h" +#include "src/__support/CPP/array.h" +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/except_value_utils.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/nearest_integer.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" + +namespace LIBC_NAMESPACE_DECL { + +static constexpr fputil::ExceptValues EXPF16_EXCEPTS_LO = {{ + // (input, RZ output, RU offset, RD offset, RN offset) + // x = 0x1.de4p-8, expf16(x) = 0x1.01cp+0 (RZ) + {0x1f79U, 0x3c07U, 1U, 0U, 0U}, + // x = 0x1.73cp-6, expf16(x) = 0x1.05cp+0 (RZ) + {0x25cfU, 0x3c17U, 1U, 0U, 0U}, +}}; + +static constexpr fputil::ExceptValues EXPF16_EXCEPTS_HI = {{ + // (input, RZ output, RU offset, RD offset, RN offset) + // x = 0x1.c34p+0, expf16(x) = 0x1.74cp+2 (RZ) + {0x3f0dU, 0x45d3U, 1U, 0U, 1U}, + // x = -0x1.488p-5, expf16(x) = 0x1.ebcp-1 (RZ) + {0xa922U, 0x3bafU, 1U, 0U, 0U}, + // x = -0x1.55p-5, expf16(x) = 0x1.ebp-1 (RZ) + {0xa954U, 0x3bacU, 1U, 0U, 0U}, +}}; + +// Generated by Sollya with the following commands: +// > display = hexadecimal; +// > for i from -18 to 12 do print(round(exp(i), SG, RN)); +static constexpr cpp::array EXP_HI = { + 0x1.05a628p-26f, 0x1.639e32p-25f, 0x1.e355bcp-24f, 0x1.4875cap-22f, + 0x1.be6c7p-21f, 0x1.2f6054p-19f, 0x1.9c54c4p-18f, 0x1.183542p-16f, + 0x1.7cd79cp-15f, 0x1.02cf22p-13f, 0x1.5fc21p-12f, 0x1.de16bap-11f, + 0x1.44e52p-9f, 0x1.b993fep-8f, 0x1.2c155cp-6f, 0x1.97db0cp-5f, + 0x1.152aaap-3f, 0x1.78b564p-2f, 0x1p+0f, 0x1.5bf0a8p+1f, + 0x1.d8e64cp+2f, 0x1.415e5cp+4f, 0x1.b4c902p+5f, 0x1.28d38ap+7f, + 0x1.936dc6p+8f, 0x1.122886p+10f, 0x1.749ea8p+11f, 0x1.fa7158p+12f, + 0x1.5829dcp+14f, 0x1.d3c448p+15f, 0x1.3de166p+17f, +}; + +// Generated by Sollya with the following commands: +// > display = hexadecimal; +// > for i from 0 to 7 do print(round(exp(i * 2^-3), SG, RN)); +static constexpr cpp::array EXP_MID = { + 0x1p+0f, 0x1.221604p+0f, 0x1.48b5e4p+0f, 0x1.747a52p+0f, + 0x1.a61298p+0f, 0x1.de455ep+0f, 0x1.0ef9dcp+1f, 0x1.330e58p+1f, +}; + +LLVM_LIBC_FUNCTION(float16, expf16, (float16 x)) { + using FPBits = fputil::FPBits; + FPBits x_bits(x); + + uint16_t x_u = x_bits.uintval(); + uint16_t x_abs = x_u & 0x7fffU; + + // When 0 < |x| <= 2^(-5), or |x| >= 12, or x is NaN. + if (LIBC_UNLIKELY(x_abs <= 0x2800U || x_abs >= 0x4a00U)) { + // exp(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; + } + + // When x >= 12. + if (x_bits.is_pos() && x_abs >= 0x4a00U) { + // exp(+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); + return FPBits::inf().get_val(); + default: + return FPBits::max_normal().get_val(); + } + } + + // When x <= -18. + if (x_u >= 0xcc80U) { + // exp(-inf) = +0 + if (x_bits.is_inf()) + return FPBits::zero().get_val(); + + fputil::set_errno_if_required(ERANGE); + fputil::raise_except_if_required(FE_UNDERFLOW | FE_INEXACT); + + switch (fputil::quick_get_round()) { + case FE_UPWARD: + return FPBits::min_subnormal().get_val(); + default: + return FPBits::zero().get_val(); + } + } + + // When 0 < |x| <= 2^(-5). + if (x_abs <= 0x2800U && !x_bits.is_zero()) { + if (auto r = EXPF16_EXCEPTS_LO.lookup(x_u); LIBC_UNLIKELY(r.has_value())) + return r.value(); + + float xf = x; + // Degree-3 minimax polynomial generated by Sollya with the following + // commands: + // > display = hexadecimal; + // > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-5, 2^-5]); + // > 1 + x * P; + float r = + fputil::polyeval(xf, 0x1p+0f, 0x1p+0f, 0x1.0004p-1f, 0x1.555778p-3f); + return static_cast(r); + } + } + + if (auto r = EXPF16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value())) + return r.value(); + + // For -18 < x < 12, to compute exp(x), we perform the following range + // reduction: find hi, mid, lo, such that: + // x = hi + mid + lo, in which + // hi is an integer, + // mid * 2^3 is an integer, + // -2^(-4) <= lo < 2^(-4). + // In particular, + // hi + mid = round(x * 2^3) * 2^(-3). + // Then, + // exp(x) = exp(hi + mid + lo) = exp(hi) * exp(mid) * exp(lo). + // We store exp(hi) and exp(mid) in the lookup tables EXP_HI and EXP_MID + // respectively. exp(lo) is computed using a degree-3 minimax polynomial + // generated by Sollya. + + float xf = static_cast(x); + float kf = fputil::nearest_integer(xf * 0x1.0p+3f); + int x_hi_mid = static_cast(kf); + int x_hi = x_hi_mid >> 3; + int x_mid = x_hi_mid & 0x7; + // lo = x - (hi + mid) = round(x * 2^3) * (-2^(-3)) + x + float lo = fputil::multiply_add(kf, -0x1.0p-3f, xf); + + float exp_hi = EXP_HI[x_hi + 18]; + float exp_mid = EXP_MID[x_mid]; + // Degree-3 minimax polynomial generated by Sollya with the following + // commands: + // > display = hexadecimal; + // > P = fpminimax(expm1(x)/x, 2, [|SG...|], [-2^-4, 2^-4]); + // > 1 + x * P; + float exp_lo = + fputil::polyeval(lo, 0x1p+0f, 0x1p+0f, 0x1.001p-1f, 0x1.555ddep-3f); + return static_cast(exp_hi * exp_mid * exp_lo); +} + +} // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/UnitTest/FPMatcher.h b/libc/test/UnitTest/FPMatcher.h index 6b50f3d23075a..2749908ef1849 100644 --- a/libc/test/UnitTest/FPMatcher.h +++ b/libc/test/UnitTest/FPMatcher.h @@ -74,7 +74,8 @@ template struct FPTest : public Test { static constexpr T inf = FPBits::inf(Sign::POS).get_val(); static constexpr T neg_inf = FPBits::inf(Sign::NEG).get_val(); static constexpr T min_normal = FPBits::min_normal().get_val(); - static constexpr T max_normal = FPBits::max_normal().get_val(); + static constexpr T max_normal = FPBits::max_normal(Sign::POS).get_val(); + static constexpr T neg_max_normal = FPBits::max_normal(Sign::NEG).get_val(); static constexpr T min_denormal = FPBits::min_subnormal().get_val(); static constexpr T max_denormal = FPBits::max_subnormal().get_val(); diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index 3ad5d98858165..d57e220f98b82 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -888,6 +888,19 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + exp_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + exp_test.cpp + DEPENDS + libc.src.errno.errno + libc.src.math.exp + libc.src.__support.FPUtil.fp_bits +) + add_fp_unittest( expf_test NEED_MPFR @@ -902,16 +915,14 @@ add_fp_unittest( ) add_fp_unittest( - exp_test - NEED_MPFR - SUITE - libc-math-unittests - SRCS - exp_test.cpp - DEPENDS - libc.src.errno.errno - libc.src.math.exp - libc.src.__support.FPUtil.fp_bits + expf16_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + expf16_test.cpp + DEPENDS + libc.src.math.expf16 ) add_fp_unittest( diff --git a/libc/test/src/math/expf16_test.cpp b/libc/test/src/math/expf16_test.cpp new file mode 100644 index 0000000000000..ee89a9c716e12 --- /dev/null +++ b/libc/test/src/math/expf16_test.cpp @@ -0,0 +1,40 @@ +//===-- Exhaustive test for expf16 ----------------------------------------===// +// +// 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/expf16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +using LlvmLibcExpf16Test = LIBC_NAMESPACE::testing::FPTest; + +namespace mpfr = LIBC_NAMESPACE::testing::mpfr; + +// Range: [0, Inf]; +static constexpr uint16_t POS_START = 0x0000U; +static constexpr uint16_t POS_STOP = 0x7c00U; + +// Range: [-Inf, 0]; +static constexpr uint16_t NEG_START = 0x8000U; +static constexpr uint16_t NEG_STOP = 0xfc00U; + +TEST_F(LlvmLibcExpf16Test, PositiveRange) { + for (uint16_t v = POS_START; v <= POS_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x, + LIBC_NAMESPACE::expf16(x), 0.5); + } +} + +TEST_F(LlvmLibcExpf16Test, NegativeRange) { + for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Exp, x, + LIBC_NAMESPACE::expf16(x), 0.5); + } +} diff --git a/libc/test/src/math/performance_testing/CMakeLists.txt b/libc/test/src/math/performance_testing/CMakeLists.txt index 536338bf2a56a..b43d21a242f35 100644 --- a/libc/test/src/math/performance_testing/CMakeLists.txt +++ b/libc/test/src/math/performance_testing/CMakeLists.txt @@ -175,6 +175,17 @@ add_perf_binary( -fno-builtin ) +add_perf_binary( + expf16_perf + SRCS + expf16_perf.cpp + DEPENDS + .single_input_single_output_diff + libc.src.math.expf16 + COMPILE_OPTIONS + -fno-builtin +) + add_perf_binary( fabsf_perf SRCS diff --git a/libc/test/src/math/performance_testing/expf16_perf.cpp b/libc/test/src/math/performance_testing/expf16_perf.cpp new file mode 100644 index 0000000000000..c1213689ff5e7 --- /dev/null +++ b/libc/test/src/math/performance_testing/expf16_perf.cpp @@ -0,0 +1,22 @@ +//===-- Performancel test for expf16 --------------------------------------===// +// +// 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 "SingleInputSingleOutputPerf.h" + +#include "src/math/expf16.h" + +// LLVM libc might be the only libc implementation with support for float16 math +// functions currently. We can't compare our float16 functions against the +// system libc, so we compare them against this placeholder function. +static float16 placeholderf16(float16 x) { return x; } + +int main() { + SINGLE_INPUT_SINGLE_OUTPUT_PERF_EX(float16, LIBC_NAMESPACE::expf16, + ::placeholderf16, 20'000, + "expf16_perf.log") +} diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index 5ddc88a38eb64..1f4e4ae483edf 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -940,6 +940,18 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + exp_test + SUITE + libc-math-smoke-tests + SRCS + exp_test.cpp + DEPENDS + libc.src.errno.errno + libc.src.math.exp + libc.src.__support.FPUtil.fp_bits +) + add_fp_unittest( expf_test SUITE @@ -953,15 +965,16 @@ add_fp_unittest( ) add_fp_unittest( - exp_test - SUITE - libc-math-smoke-tests - SRCS - exp_test.cpp - DEPENDS - libc.src.errno.errno - libc.src.math.exp - libc.src.__support.FPUtil.fp_bits + expf16_test + SUITE + libc-math-smoke-tests + SRCS + expf16_test.cpp + DEPENDS + libc.hdr.errno_macros + libc.hdr.fenv_macros + libc.src.errno.errno + libc.src.math.expf16 ) add_fp_unittest( diff --git a/libc/test/src/math/smoke/expf16_test.cpp b/libc/test/src/math/smoke/expf16_test.cpp new file mode 100644 index 0000000000000..f05ecd0dc4d0e --- /dev/null +++ b/libc/test/src/math/smoke/expf16_test.cpp @@ -0,0 +1,66 @@ +//===-- Unittests for expf16 ----------------------------------------------===// +// +// 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 "hdr/errno_macros.h" +#include "hdr/fenv_macros.h" +#include "src/errno/libc_errno.h" +#include "src/math/expf16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" + +using LlvmLibcExpf16Test = LIBC_NAMESPACE::testing::FPTest; + +TEST_F(LlvmLibcExpf16Test, SpecialNumbers) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::expf16(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ_WITH_EXCEPTION(sNaN, LIBC_NAMESPACE::expf16(sNaN), FE_INVALID); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::expf16(inf)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ_ALL_ROUNDING(static_cast(zero), + LIBC_NAMESPACE::expf16(neg_inf)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ_ALL_ROUNDING(static_cast(1.0f), + LIBC_NAMESPACE::expf16(zero)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ_ALL_ROUNDING(static_cast(1.0f), + LIBC_NAMESPACE::expf16(neg_zero)); + EXPECT_MATH_ERRNO(0); +} + +TEST_F(LlvmLibcExpf16Test, Overflow) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::expf16(max_normal), + FE_OVERFLOW); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ_WITH_EXCEPTION( + inf, LIBC_NAMESPACE::expf16(static_cast(12.0)), FE_OVERFLOW); + EXPECT_MATH_ERRNO(ERANGE); +} + +TEST_F(LlvmLibcExpf16Test, Underflow) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ_WITH_EXCEPTION(zero, LIBC_NAMESPACE::expf16(neg_max_normal), + FE_UNDERFLOW | FE_INEXACT); + EXPECT_MATH_ERRNO(ERANGE); + + EXPECT_FP_EQ_WITH_EXCEPTION( + zero, LIBC_NAMESPACE::expf16(static_cast(-18.0)), + FE_UNDERFLOW | FE_INEXACT); + EXPECT_MATH_ERRNO(ERANGE); +}