[libc-commits] [libc] [libc][math][c23] Add expm1f16 C23 math function (PR #102387)
via libc-commits
libc-commits at lists.llvm.org
Mon Aug 12 09:24:33 PDT 2024
https://github.com/overmighty updated https://github.com/llvm/llvm-project/pull/102387
>From eff043c923d4b3d243c3e8ec56d38b067850275a Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Wed, 7 Aug 2024 23:50:26 +0200
Subject: [PATCH 1/3] [libc][math][c23] Add expm1f16 C23 math function
Part of #95250.
---
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/docs/math/index.rst | 2 +-
libc/spec/stdc.td | 1 +
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/expm1f16.h | 21 ++++
libc/src/math/generic/CMakeLists.txt | 28 ++++-
libc/src/math/generic/expf16.cpp | 61 +--------
libc/src/math/generic/expm1f16.cpp | 139 +++++++++++++++++++++
libc/src/math/generic/expxf16.h | 65 ++++++++++
libc/test/src/math/CMakeLists.txt | 31 +++--
libc/test/src/math/expm1f16_test.cpp | 40 ++++++
libc/test/src/math/smoke/CMakeLists.txt | 21 +++-
libc/test/src/math/smoke/expm1f16_test.cpp | 65 ++++++++++
13 files changed, 401 insertions(+), 75 deletions(-)
create mode 100644 libc/src/math/expm1f16.h
create mode 100644 libc/src/math/generic/expm1f16.cpp
create mode 100644 libc/test/src/math/expm1f16_test.cpp
create mode 100644 libc/test/src/math/smoke/expm1f16_test.cpp
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 39d90bf06a0cf4..332a90c84fbdc9 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -584,6 +584,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.exp10f16
libc.src.math.exp2f16
libc.src.math.expf16
+ libc.src.math.expm1f16
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 ce83571e2b3815..642776ddf49916 100644
--- a/libc/docs/math/index.rst
+++ b/libc/docs/math/index.rst
@@ -294,7 +294,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| exp2m1 | |check| | | | | | 7.12.6.5 | F.10.3.5 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| expm1 | |check| | |check| | | | | 7.12.6.6 | F.10.3.6 |
+| expm1 | |check| | |check| | | |check| | | 7.12.6.6 | F.10.3.6 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| fma | |check| | |check| | | | | 7.12.13.1 | F.10.10.1 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/spec/stdc.td b/libc/spec/stdc.td
index d24d43526fe5c7..76555439a2dc8a 100644
--- a/libc/spec/stdc.td
+++ b/libc/spec/stdc.td
@@ -594,6 +594,7 @@ def StdC : StandardSpec<"stdc"> {
FunctionSpec<"expm1", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"expm1f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
+ GuardedFunctionSpec<"expm1f16", RetValSpec<Float16Type>, [ArgSpec<Float16Type>], "LIBC_TYPES_HAS_FLOAT16">,
FunctionSpec<"exp10", RetValSpec<DoubleType>, [ArgSpec<DoubleType>]>,
FunctionSpec<"exp10f", RetValSpec<FloatType>, [ArgSpec<FloatType>]>,
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index a3c9b5473c290e..9ccc3f1a347f29 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -120,6 +120,7 @@ add_math_entrypoint_object(exp10f16)
add_math_entrypoint_object(expm1)
add_math_entrypoint_object(expm1f)
+add_math_entrypoint_object(expm1f16)
add_math_entrypoint_object(f16add)
add_math_entrypoint_object(f16addf)
diff --git a/libc/src/math/expm1f16.h b/libc/src/math/expm1f16.h
new file mode 100644
index 00000000000000..644e6cddd7666a
--- /dev/null
+++ b/libc/src/math/expm1f16.h
@@ -0,0 +1,21 @@
+//===-- 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_MATH_EXPM1F16_H
+#define LLVM_LIBC_SRC_MATH_EXPM1F16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 expm1f16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_EXPM1F16_H
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index 77f7f4fef007b9..469df67951ab33 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -1322,14 +1322,12 @@ add_entrypoint_object(
HDRS
../expf16.h
DEPENDS
+ .expxf16
libc.hdr.errno_macros
libc.hdr.fenv_macros
- libc.src.__support.CPP.array
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.nearest_integer
libc.src.__support.FPUtil.polyeval
libc.src.__support.FPUtil.rounding_mode
libc.src.__support.macros.optimization
@@ -1571,6 +1569,27 @@ add_entrypoint_object(
-O3
)
+add_entrypoint_object(
+ expm1f16
+ SRCS
+ expm1f16.cpp
+ HDRS
+ ../expm1f16.h
+ DEPENDS
+ .expxf16
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ 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
+ COMPILE_OPTIONS
+ -O3
+)
+
add_entrypoint_object(
powf
SRCS
@@ -4878,4 +4897,7 @@ add_header_library(
expxf16.h
DEPENDS
libc.src.__support.CPP.array
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.FPUtil.nearest_integer
+ libc.src.__support.FPUtil.polyeval
)
diff --git a/libc/src/math/generic/expf16.cpp b/libc/src/math/generic/expf16.cpp
index b198c559dfedb9..7ffdbd5191008a 100644
--- a/libc/src/math/generic/expf16.cpp
+++ b/libc/src/math/generic/expf16.cpp
@@ -7,15 +7,13 @@
//===----------------------------------------------------------------------===//
#include "src/math/expf16.h"
+#include "expxf16.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"
@@ -41,28 +39,6 @@ static constexpr fputil::ExceptValues<float16, 3> EXPF16_EXCEPTS_HI = {{
{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<float, 31> 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<float, 8> 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<float16>;
FPBits x_bits(x);
@@ -135,38 +111,9 @@ LLVM_LIBC_FUNCTION(float16, expf16, (float16 x)) {
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 = x;
- float kf = fputil::nearest_integer(xf * 0x1.0p+3f);
- int x_hi_mid = static_cast<int>(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<float16>(exp_hi * exp_mid * exp_lo);
+ // exp(x) = exp(hi + mid) * exp(lo)
+ auto [exp_hi_mid, exp_lo] = exp_range_reduction(x);
+ return static_cast<float16>(exp_hi_mid * exp_lo);
}
} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/expm1f16.cpp b/libc/src/math/generic/expm1f16.cpp
new file mode 100644
index 00000000000000..f13f6850db5244
--- /dev/null
+++ b/libc/src/math/generic/expm1f16.cpp
@@ -0,0 +1,139 @@
+//===-- Half-precision e^x - 1 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/expm1f16.h"
+#include "expxf16.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/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"
+
+namespace LIBC_NAMESPACE_DECL {
+
+static constexpr fputil::ExceptValues<float16, 1> 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
+static constexpr size_t N_EXPM1F16_EXCEPTS_HI = 2;
+#else
+static constexpr size_t N_EXPM1F16_EXCEPTS_HI = 3;
+#endif
+
+static constexpr fputil::ExceptValues<float16, N_EXPM1F16_EXCEPTS_HI>
+ 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
+ // x = 0x1.a08p-3, exp10m1f(x) = 0x1.cdcp-3 (RZ)
+ {0x3282U, 0x3337U, 1U, 0U, 0U},
+#endif
+ }};
+
+LLVM_LIBC_FUNCTION(float16, expm1f16, (float16 x)) {
+ using FPBits = fputil::FPBits<float16>;
+ 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 >= 12.
+ if (x_bits.is_pos() && x_abs >= 0x4a00U) {
+ // 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);
+ 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();
+
+ fputil::set_errno_if_required(ERANGE);
+ fputil::raise_except_if_required(FE_INEXACT);
+
+ switch (fputil::quick_get_round()) {
+ case FE_UPWARD:
+ case FE_TOWARDZERO:
+ return static_cast<float16>(-0x1.ffcp-1);
+ case FE_TONEAREST:
+ if (x_u > 0xc828U)
+ return FPBits::one(Sign::NEG).get_val();
+ return static_cast<float16>(-0x1.ffcp-1);
+ default:
+ return FPBits::one(Sign::NEG).get_val();
+ }
+ }
+
+ // When 0 < |x| <= 2^(-3).
+ if (x_abs <= 0x3000U && !x_bits.is_zero()) {
+ if (auto r = EXPM1F16_EXCEPTS_LO.lookup(x_u);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+
+ 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 static_cast<float16>(
+ xf * fputil::polyeval(xf, 0x1p+0f, 0x1.fffff8p-2f, 0x1.555556p-3f,
+ 0x1.55905ep-5f, 0x1.1124c2p-7f));
+ }
+ }
+
+ if (auto r = EXPM1F16_EXCEPTS_HI.lookup(x_u); LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+
+ // 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 static_cast<float16>(fputil::multiply_add(exp_hi_mid, exp_lo, -1.0f));
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/expxf16.h b/libc/src/math/generic/expxf16.h
index c33aca337b98dc..53815e0e275533 100644
--- a/libc/src/math/generic/expxf16.h
+++ b/libc/src/math/generic/expxf16.h
@@ -10,11 +10,76 @@
#define LLVM_LIBC_SRC_MATH_GENERIC_EXPXF16_H
#include "src/__support/CPP/array.h"
+#include "src/__support/FPUtil/PolyEval.h"
+#include "src/__support/FPUtil/multiply_add.h"
+#include "src/__support/FPUtil/nearest_integer.h"
#include "src/__support/macros/config.h"
#include <stdint.h>
namespace LIBC_NAMESPACE_DECL {
+// 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<float, 31> 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<float, 8> 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,
+};
+
+struct ExpRangeReduction {
+ float exp_hi_mid;
+ float exp_lo;
+};
+
+ExpRangeReduction exp_range_reduction(float16 x) {
+ // 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 = x;
+ float kf = fputil::nearest_integer(xf * 0x1.0p+3f);
+ int x_hi_mid = static_cast<int>(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 {exp_hi * exp_mid, exp_lo};
+}
+
// Generated by Sollya with the following commands:
// > display = hexadecimal;
// > for i from 0 to 7 do printsingle(round(2^(i * 2^-3), SG, RN));
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index d209c9a38499b6..b6b6953a688716 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -1663,6 +1663,19 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ expm1_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ expm1_test.cpp
+ DEPENDS
+ libc.src.errno.errno
+ libc.src.math.expm1
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
expm1f_test
NEED_MPFR
@@ -1677,16 +1690,14 @@ add_fp_unittest(
)
add_fp_unittest(
- expm1_test
- NEED_MPFR
- SUITE
- libc-math-unittests
- SRCS
- expm1_test.cpp
- DEPENDS
- libc.src.errno.errno
- libc.src.math.expm1
- libc.src.__support.FPUtil.fp_bits
+ expm1f16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ expm1f16_test.cpp
+ DEPENDS
+ libc.src.math.expm1f16
)
add_fp_unittest(
diff --git a/libc/test/src/math/expm1f16_test.cpp b/libc/test/src/math/expm1f16_test.cpp
new file mode 100644
index 00000000000000..a6a6fcf73d383f
--- /dev/null
+++ b/libc/test/src/math/expm1f16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for expm1f16 --------------------------------------===//
+//
+// 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/expm1f16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcExpm1f16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+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(LlvmLibcExpm1f16Test, PositiveRange) {
+ for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Expm1, x,
+ LIBC_NAMESPACE::expm1f16(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcExpm1f16Test, NegativeRange) {
+ for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Expm1, x,
+ LIBC_NAMESPACE::expm1f16(x), 0.5);
+ }
+}
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 9aa5399aae9d38..18d8eb27aec235 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -3294,6 +3294,18 @@ add_fp_unittest(
libc.src.math.fma
)
+add_fp_unittest(
+ expm1_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ expm1_test.cpp
+ DEPENDS
+ libc.src.errno.errno
+ libc.src.math.expm1
+ libc.src.__support.FPUtil.fp_bits
+)
+
add_fp_unittest(
expm1f_test
SUITE
@@ -3307,15 +3319,16 @@ add_fp_unittest(
)
add_fp_unittest(
- expm1_test
+ expm1f16_test
SUITE
libc-math-smoke-tests
SRCS
- expm1_test.cpp
+ expm1f16_test.cpp
DEPENDS
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
libc.src.errno.errno
- libc.src.math.expm1
- libc.src.__support.FPUtil.fp_bits
+ libc.src.math.expm1f16
)
add_fp_unittest(
diff --git a/libc/test/src/math/smoke/expm1f16_test.cpp b/libc/test/src/math/smoke/expm1f16_test.cpp
new file mode 100644
index 00000000000000..959fa517a5046d
--- /dev/null
+++ b/libc/test/src/math/smoke/expm1f16_test.cpp
@@ -0,0 +1,65 @@
+//===-- Unittests for expm1f16 --------------------------------------------===//
+//
+// 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/expm1f16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcExpm1f16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcExpm1f16Test, SpecialNumbers) {
+ LIBC_NAMESPACE::libc_errno = 0;
+
+ EXPECT_FP_EQ_ALL_ROUNDING(aNaN, LIBC_NAMESPACE::expm1f16(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(aNaN, LIBC_NAMESPACE::expm1f16(sNaN), FE_INVALID);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(inf, LIBC_NAMESPACE::expm1f16(inf));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(static_cast<float16>(-1.0),
+ LIBC_NAMESPACE::expm1f16(neg_inf));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(zero, LIBC_NAMESPACE::expm1f16(zero));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_ALL_ROUNDING(neg_zero, LIBC_NAMESPACE::expm1f16(neg_zero));
+ EXPECT_MATH_ERRNO(0);
+}
+
+TEST_F(LlvmLibcExpm1f16Test, Overflow) {
+ LIBC_NAMESPACE::libc_errno = 0;
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::expm1f16(max_normal),
+ FE_OVERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ inf, LIBC_NAMESPACE::expm1f16(static_cast<float16>(12.0)), FE_OVERFLOW);
+ EXPECT_MATH_ERRNO(ERANGE);
+}
+
+TEST_F(LlvmLibcExpm1f16Test, ResultNearNegOne) {
+ LIBC_NAMESPACE::libc_errno = 0;
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(static_cast<float16>(-1.0),
+ LIBC_NAMESPACE::expm1f16(neg_max_normal),
+ FE_INEXACT);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION(
+ static_cast<float16>(-1.0),
+ LIBC_NAMESPACE::expm1f16(static_cast<float16>(-9.0)), FE_INEXACT);
+ EXPECT_MATH_ERRNO(ERANGE);
+}
>From 0ada06a58fd1c8a06251f7f1473133a1b0e41ca6 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Thu, 8 Aug 2024 19:46:07 +0200
Subject: [PATCH 2/3] Fix ERANGE being set when output is close to -1, and
simplify rounding near -1
---
libc/src/math/generic/expm1f16.cpp | 21 ++++-------
libc/test/src/math/smoke/expm1f16_test.cpp | 41 +++++++++++++++++++---
2 files changed, 43 insertions(+), 19 deletions(-)
diff --git a/libc/src/math/generic/expm1f16.cpp b/libc/src/math/generic/expm1f16.cpp
index f13f6850db5244..a8e974bd70884d 100644
--- a/libc/src/math/generic/expm1f16.cpp
+++ b/libc/src/math/generic/expm1f16.cpp
@@ -93,20 +93,13 @@ LLVM_LIBC_FUNCTION(float16, expm1f16, (float16 x)) {
if (x_bits.is_inf())
return FPBits::one(Sign::NEG).get_val();
- fputil::set_errno_if_required(ERANGE);
- fputil::raise_except_if_required(FE_INEXACT);
-
- switch (fputil::quick_get_round()) {
- case FE_UPWARD:
- case FE_TOWARDZERO:
- return static_cast<float16>(-0x1.ffcp-1);
- case FE_TONEAREST:
- if (x_u > 0xc828U)
- return FPBits::one(Sign::NEG).get_val();
- return static_cast<float16>(-0x1.ffcp-1);
- default:
- 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(
+ static_cast<float16>(-0x1.ffcp-1));
}
// When 0 < |x| <= 2^(-3).
diff --git a/libc/test/src/math/smoke/expm1f16_test.cpp b/libc/test/src/math/smoke/expm1f16_test.cpp
index 959fa517a5046d..1270b85d2cd5ad 100644
--- a/libc/test/src/math/smoke/expm1f16_test.cpp
+++ b/libc/test/src/math/smoke/expm1f16_test.cpp
@@ -56,10 +56,41 @@ TEST_F(LlvmLibcExpm1f16Test, ResultNearNegOne) {
EXPECT_FP_EQ_WITH_EXCEPTION(static_cast<float16>(-1.0),
LIBC_NAMESPACE::expm1f16(neg_max_normal),
FE_INEXACT);
- EXPECT_MATH_ERRNO(ERANGE);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- static_cast<float16>(-1.0),
- LIBC_NAMESPACE::expm1f16(static_cast<float16>(-9.0)), FE_INEXACT);
- EXPECT_MATH_ERRNO(ERANGE);
+ // round(-11 * log(2), HP, RN);
+ float16 x = static_cast<float16>(-0x1.e8p+2);
+
+ EXPECT_FP_EQ_ROUNDING_NEAREST(static_cast<float16>(-0x1.ffcp-1),
+ LIBC_NAMESPACE::expm1f16(x));
+ EXPECT_FP_EXCEPTION(FE_INEXACT);
+
+ EXPECT_FP_EQ_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
+ LIBC_NAMESPACE::expm1f16(x));
+ EXPECT_FP_EXCEPTION(FE_INEXACT);
+
+ EXPECT_FP_EQ_ROUNDING_DOWNWARD(static_cast<float16>(-1.0),
+ LIBC_NAMESPACE::expm1f16(x));
+ EXPECT_FP_EXCEPTION(FE_INEXACT);
+
+ EXPECT_FP_EQ_ROUNDING_TOWARD_ZERO(static_cast<float16>(-0x1.ffcp-1),
+ LIBC_NAMESPACE::expm1f16(x));
+ EXPECT_FP_EXCEPTION(FE_INEXACT);
+
+ x = static_cast<float16>(-0x1.0a4p+3);
+
+ EXPECT_FP_EQ_ROUNDING_NEAREST(static_cast<float16>(-1.0),
+ LIBC_NAMESPACE::expm1f16(x));
+ EXPECT_FP_EXCEPTION(FE_INEXACT);
+
+ EXPECT_FP_EQ_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
+ LIBC_NAMESPACE::expm1f16(x));
+ EXPECT_FP_EXCEPTION(FE_INEXACT);
+
+ EXPECT_FP_EQ_ROUNDING_DOWNWARD(static_cast<float16>(-1.0),
+ LIBC_NAMESPACE::expm1f16(x));
+ EXPECT_FP_EXCEPTION(FE_INEXACT);
+
+ EXPECT_FP_EQ_ROUNDING_TOWARD_ZERO(static_cast<float16>(-0x1.ffcp-1),
+ LIBC_NAMESPACE::expm1f16(x));
+ EXPECT_FP_EXCEPTION(FE_INEXACT);
}
>From fdeeb7b3ccee43b00cdde0e136ed495409eb8cc2 Mon Sep 17 00:00:00 2001
From: OverMighty <its.overmighty at gmail.com>
Date: Mon, 12 Aug 2024 18:20:53 +0200
Subject: [PATCH 3/3] Fix overflow cases
---
libc/src/math/generic/expm1f16.cpp | 6 +-
libc/test/UnitTest/FPMatcher.h | 32 +++++++++++
libc/test/src/math/smoke/expm1f16_test.cpp | 66 +++++++++++++---------
3 files changed, 74 insertions(+), 30 deletions(-)
diff --git a/libc/src/math/generic/expm1f16.cpp b/libc/src/math/generic/expm1f16.cpp
index a8e974bd70884d..0facdc510e4287 100644
--- a/libc/src/math/generic/expm1f16.cpp
+++ b/libc/src/math/generic/expm1f16.cpp
@@ -70,8 +70,8 @@ LLVM_LIBC_FUNCTION(float16, expm1f16, (float16 x)) {
if (x_abs == 0)
return x;
- // When x >= 12.
- if (x_bits.is_pos() && x_abs >= 0x4a00U) {
+ // 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();
@@ -80,7 +80,7 @@ LLVM_LIBC_FUNCTION(float16, expm1f16, (float16 x)) {
case FE_TONEAREST:
case FE_UPWARD:
fputil::set_errno_if_required(ERANGE);
- fputil::raise_except_if_required(FE_OVERFLOW);
+ fputil::raise_except_if_required(FE_OVERFLOW | FE_INEXACT);
return FPBits::inf().get_val();
default:
return FPBits::max_normal().get_val();
diff --git a/libc/test/UnitTest/FPMatcher.h b/libc/test/UnitTest/FPMatcher.h
index 2749908ef18495..43752a4942ad56 100644
--- a/libc/test/UnitTest/FPMatcher.h
+++ b/libc/test/UnitTest/FPMatcher.h
@@ -234,4 +234,36 @@ template <typename T> struct FPTest : public Test {
#define EXPECT_FP_EQ_ROUNDING_TOWARD_ZERO(expected, actual) \
EXPECT_FP_EQ_ROUNDING_MODE((expected), (actual), RoundingMode::TowardZero)
+#define EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_MODE( \
+ expected, actual, expected_except, rounding_mode) \
+ do { \
+ using namespace LIBC_NAMESPACE::fputil::testing; \
+ ForceRoundingMode __r((rounding_mode)); \
+ if (__r.success) { \
+ LIBC_NAMESPACE::fputil::clear_except(FE_ALL_EXCEPT); \
+ EXPECT_FP_EQ((expected), (actual)); \
+ EXPECT_FP_EXCEPTION(expected_except); \
+ } \
+ } while (0)
+
+#define EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(expected, actual, \
+ expected_except) \
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_MODE( \
+ (expected), (actual), (expected_except), RoundingMode::Nearest)
+
+#define EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(expected, actual, \
+ expected_except) \
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_MODE( \
+ (expected), (actual), (expected_except), RoundingMode::Upward)
+
+#define EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(expected, actual, \
+ expected_except) \
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_MODE( \
+ (expected), (actual), (expected_except), RoundingMode::Downward)
+
+#define EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(expected, actual, \
+ expected_except) \
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_MODE( \
+ (expected), (actual), (expected_except), RoundingMode::TowardZero)
+
#endif // LLVM_LIBC_TEST_UNITTEST_FPMATCHER_H
diff --git a/libc/test/src/math/smoke/expm1f16_test.cpp b/libc/test/src/math/smoke/expm1f16_test.cpp
index 1270b85d2cd5ad..3bdbaad2279416 100644
--- a/libc/test/src/math/smoke/expm1f16_test.cpp
+++ b/libc/test/src/math/smoke/expm1f16_test.cpp
@@ -42,12 +42,27 @@ TEST_F(LlvmLibcExpm1f16Test, Overflow) {
LIBC_NAMESPACE::libc_errno = 0;
EXPECT_FP_EQ_WITH_EXCEPTION(inf, LIBC_NAMESPACE::expm1f16(max_normal),
- FE_OVERFLOW);
+ FE_OVERFLOW | FE_INEXACT);
EXPECT_MATH_ERRNO(ERANGE);
- EXPECT_FP_EQ_WITH_EXCEPTION(
- inf, LIBC_NAMESPACE::expm1f16(static_cast<float16>(12.0)), FE_OVERFLOW);
+ // round(16 * log(2), HP, RN);
+ float16 x = static_cast<float16>(0x1.63p+3);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(inf, LIBC_NAMESPACE::expm1f16(x),
+ FE_OVERFLOW | FE_INEXACT);
+ EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(inf, LIBC_NAMESPACE::expm1f16(x),
+ FE_OVERFLOW | FE_INEXACT);
EXPECT_MATH_ERRNO(ERANGE);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+ max_normal, LIBC_NAMESPACE::expm1f16(x), FE_INEXACT);
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+ max_normal, LIBC_NAMESPACE::expm1f16(x), FE_INEXACT);
+ EXPECT_MATH_ERRNO(0);
}
TEST_F(LlvmLibcExpm1f16Test, ResultNearNegOne) {
@@ -60,37 +75,34 @@ TEST_F(LlvmLibcExpm1f16Test, ResultNearNegOne) {
// round(-11 * log(2), HP, RN);
float16 x = static_cast<float16>(-0x1.e8p+2);
- EXPECT_FP_EQ_ROUNDING_NEAREST(static_cast<float16>(-0x1.ffcp-1),
- LIBC_NAMESPACE::expm1f16(x));
- EXPECT_FP_EXCEPTION(FE_INEXACT);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+ static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::expm1f16(x),
+ FE_INEXACT);
- EXPECT_FP_EQ_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
- LIBC_NAMESPACE::expm1f16(x));
- EXPECT_FP_EXCEPTION(FE_INEXACT);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
+ LIBC_NAMESPACE::expm1f16(x),
+ FE_INEXACT);
- EXPECT_FP_EQ_ROUNDING_DOWNWARD(static_cast<float16>(-1.0),
- LIBC_NAMESPACE::expm1f16(x));
- EXPECT_FP_EXCEPTION(FE_INEXACT);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+ static_cast<float16>(-1.0), LIBC_NAMESPACE::expm1f16(x), FE_INEXACT);
- EXPECT_FP_EQ_ROUNDING_TOWARD_ZERO(static_cast<float16>(-0x1.ffcp-1),
- LIBC_NAMESPACE::expm1f16(x));
- EXPECT_FP_EXCEPTION(FE_INEXACT);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+ static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::expm1f16(x),
+ FE_INEXACT);
x = static_cast<float16>(-0x1.0a4p+3);
- EXPECT_FP_EQ_ROUNDING_NEAREST(static_cast<float16>(-1.0),
- LIBC_NAMESPACE::expm1f16(x));
- EXPECT_FP_EXCEPTION(FE_INEXACT);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_NEAREST(
+ static_cast<float16>(-1.0), LIBC_NAMESPACE::expm1f16(x), FE_INEXACT);
- EXPECT_FP_EQ_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
- LIBC_NAMESPACE::expm1f16(x));
- EXPECT_FP_EXCEPTION(FE_INEXACT);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_UPWARD(static_cast<float16>(-0x1.ffcp-1),
+ LIBC_NAMESPACE::expm1f16(x),
+ FE_INEXACT);
- EXPECT_FP_EQ_ROUNDING_DOWNWARD(static_cast<float16>(-1.0),
- LIBC_NAMESPACE::expm1f16(x));
- EXPECT_FP_EXCEPTION(FE_INEXACT);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_DOWNWARD(
+ static_cast<float16>(-1.0), LIBC_NAMESPACE::expm1f16(x), FE_INEXACT);
- EXPECT_FP_EQ_ROUNDING_TOWARD_ZERO(static_cast<float16>(-0x1.ffcp-1),
- LIBC_NAMESPACE::expm1f16(x));
- EXPECT_FP_EXCEPTION(FE_INEXACT);
+ EXPECT_FP_EQ_WITH_EXCEPTION_ROUNDING_TOWARD_ZERO(
+ static_cast<float16>(-0x1.ffcp-1), LIBC_NAMESPACE::expm1f16(x),
+ FE_INEXACT);
}
More information about the libc-commits
mailing list