[libc-commits] [libc] [libc][math][c23] Add tanf16 function (PR #121018)
via libc-commits
libc-commits at lists.llvm.org
Tue Dec 31 04:05:44 PST 2024
https://github.com/wldfngrs updated https://github.com/llvm/llvm-project/pull/121018
>From f8a5f5a13f51cd817c5a827b6178a0413dab6287 Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Tue, 24 Dec 2024 02:20:44 +0100
Subject: [PATCH 1/2] Add tanf16 function
---
libc/config/linux/x86_64/entrypoints.txt | 1 +
libc/hdrgen/yaml/math.yaml | 7 ++
libc/src/math/CMakeLists.txt | 1 +
libc/src/math/generic/CMakeLists.txt | 19 ++++
libc/src/math/generic/sincosf16_utils.h | 7 +-
libc/src/math/generic/tanf16.cpp | 112 +++++++++++++++++++++++
libc/src/math/generic/tanpif16.cpp | 2 +-
libc/src/math/tanf16.h | 21 +++++
libc/test/src/math/CMakeLists.txt | 11 +++
libc/test/src/math/cosf16_test.cpp | 2 +-
libc/test/src/math/smoke/CMakeLists.txt | 11 +++
libc/test/src/math/smoke/tanf16_test.cpp | 34 +++++++
libc/test/src/math/tanf16_test.cpp | 40 ++++++++
13 files changed, 263 insertions(+), 5 deletions(-)
create mode 100644 libc/src/math/generic/tanf16.cpp
create mode 100644 libc/src/math/tanf16.h
create mode 100644 libc/test/src/math/smoke/tanf16_test.cpp
create mode 100644 libc/test/src/math/tanf16_test.cpp
diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt
index 08d8559d8c81a2..782f4ef9b99ba1 100644
--- a/libc/config/linux/x86_64/entrypoints.txt
+++ b/libc/config/linux/x86_64/entrypoints.txt
@@ -714,6 +714,7 @@ if(LIBC_TYPES_HAS_FLOAT16)
libc.src.math.sinhf16
libc.src.math.sinpif16
libc.src.math.sqrtf16
+ libc.src.math.tanf16
libc.src.math.tanhf16
libc.src.math.tanpif16
libc.src.math.totalorderf16
diff --git a/libc/hdrgen/yaml/math.yaml b/libc/hdrgen/yaml/math.yaml
index 3b8caec66bbfd2..14e527f4a6f3d4 100644
--- a/libc/hdrgen/yaml/math.yaml
+++ b/libc/hdrgen/yaml/math.yaml
@@ -2417,6 +2417,13 @@ functions:
return_type: float
arguments:
- type: float
+ - name: tanf16
+ standards:
+ - stdc
+ return_type: _Float16
+ arguments:
+ - type: _Float16
+ guard: LIBC_TYPES_HAS_FLOAT16
- name: tanhf
standards:
- stdc
diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt
index e4e2c49642f2d0..fe5ebd793b40af 100644
--- a/libc/src/math/CMakeLists.txt
+++ b/libc/src/math/CMakeLists.txt
@@ -501,6 +501,7 @@ add_math_entrypoint_object(sqrtf128)
add_math_entrypoint_object(tan)
add_math_entrypoint_object(tanf)
+add_math_entrypoint_object(tanf16)
add_math_entrypoint_object(tanh)
add_math_entrypoint_object(tanhf)
diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt
index b3d46129151974..824a161df1b19f 100644
--- a/libc/src/math/generic/CMakeLists.txt
+++ b/libc/src/math/generic/CMakeLists.txt
@@ -662,6 +662,25 @@ add_entrypoint_object(
${libc_opt_high_flag}
)
+add_entrypoint_object(
+ tanf16
+ SRCS
+ tanf16.cpp
+ HDRS
+ ../tanf16.h
+ DEPENDS
+ .sincosf16_utils
+ libc.hdr.errno_macros
+ libc.hdr.fenv_macros
+ libc.src.__support.FPUtil.cast
+ libc.src.__support.FPUtil.fenv_impl
+ libc.src.__support.FPUtil.fp_bits
+ libc.src.__support.FPUtil.except_value_utils
+ libc.src.__support.FPUtil.multiply_add
+ libc.src.__support.macros.optimization
+ libc.src.__support.macros.properties.types
+)
+
add_entrypoint_object(
tanpif16
SRCS
diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h
index 5e5edd4a8c85bd..48feee590937fd 100644
--- a/libc/src/math/generic/sincosf16_utils.h
+++ b/libc/src/math/generic/sincosf16_utils.h
@@ -63,10 +63,11 @@ LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) {
// further intermediate computation.
LIBC_INLINE int32_t range_reduction_sincosf16(float x, float &y) {
double prod = x * 0x1.45f306dc9c883p3;
- double kf = fputil::nearest_integer(prod);
- y = static_cast<float>(prod - kf);
+ double kd = fputil::nearest_integer(prod);
- return static_cast<int32_t>(kf);
+ y = static_cast<float>(prod - kd);
+
+ return static_cast<int32_t>(kd);
}
static LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k,
diff --git a/libc/src/math/generic/tanf16.cpp b/libc/src/math/generic/tanf16.cpp
new file mode 100644
index 00000000000000..bae77f9f331871
--- /dev/null
+++ b/libc/src/math/generic/tanf16.cpp
@@ -0,0 +1,112 @@
+//===-- Half-precision tan(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/tanf16.h"
+#include "hdr/errno_macros.h"
+#include "hdr/fenv_macros.h"
+#include "sincosf16_utils.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/multiply_add.h"
+#include "src/__support/macros/optimization.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+constexpr size_t N_EXCEPTS = 9;
+
+constexpr fputil::ExceptValues<float16, N_EXCEPTS> TANF16_EXCEPTS{{
+ // (input, RZ output, RU offset, RD offset, RN offset)
+ {0x2894, 0x2894, 1, 0, 1},
+ {0x3091, 0x3099, 1, 0, 0},
+ {0x3098, 0x30a0, 1, 0, 0},
+ {0x55ed, 0x3911, 1, 0, 0},
+ {0x607b, 0xc638, 0, 1, 1},
+ {0x674e, 0x3b7d, 1, 0, 0},
+ {0x6807, 0x4014, 1, 0, 1},
+ {0x6f4d, 0xbe19, 0, 1, 1},
+ {0x7330, 0xcb62, 0, 1, 0},
+}};
+
+LLVM_LIBC_FUNCTION(float16, tanf16, (float16 x)) {
+ using FPBits = fputil::FPBits<float16>;
+ FPBits xbits(x);
+
+ uint16_t x_u = xbits.uintval();
+ uint16_t x_abs = x_u & 0x7fff;
+ bool x_sign = x_u >> 15;
+ float xf = x;
+
+ // Handle exceptional values
+ if (auto r = TANF16_EXCEPTS.lookup_odd(x_abs, x_sign);
+ LIBC_UNLIKELY(r.has_value()))
+ return r.value();
+
+ // |x| <= 0x1.d1p-5
+ if (LIBC_UNLIKELY(x_abs <= 0x2b44)) {
+ if (LIBC_UNLIKELY(x_abs <= 0x10e6)) {
+ // tan(+/-0) = +/-0
+ if (LIBC_UNLIKELY(x_abs == 0U))
+ return x;
+
+ int rounding = fputil::quick_get_round();
+
+ // Exhaustive tests show that, when:
+ // x > 0, and rounding upward or
+ // x < 0, and rounding downward then,
+ // tan(x) = x * 2^-11 + x
+ if ((xbits.is_pos() && rounding == FE_UPWARD) ||
+ (xbits.is_neg() && rounding == FE_DOWNWARD))
+ return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
+ else
+ return x;
+ }
+
+ float xsq = xf * xf;
+
+ float result = fputil::polyeval(xsq, 0x1p0f, 0x1.555556p-2f, 0x1.110ee4p-3f,
+ 0x1.be80f6p-5f);
+
+ return fputil::cast<float16>(xf * result);
+ }
+
+ // tan(+/-inf)= NaN, and tan(NaN) = NaN
+ if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
+ if (x_abs == 0x7c00) {
+ fputil::set_errno_if_required(EDOM);
+ fputil::raise_except_if_required(FE_INVALID);
+ }
+
+ return x + FPBits::quiet_nan().get_val();
+ }
+
+ // Range reduction:
+ // For |x| > pi/32, we perform range reduction as follows:
+ // Find k and y such that:
+ // x = (k + y) * pi/32;
+ // k is an integer, |y| < 0.5
+ //
+ // This is done by performing:
+ // k = round(x * 32/pi)
+ // y = x * 32/pi - k
+ //
+ // Once k and y are computed, we then deduce the answer by the formula:
+ // tan(x) = sin(x) / cos(x)
+ // = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)
+ float sin_k, cos_k, sin_y, cosm1_y;
+ sincosf16_eval(xf, sin_k, cos_k, sin_y, cosm1_y);
+
+ // Note that, cosm1_y = cos_y - 1:
+ using fputil::multiply_add;
+ return fputil::cast<float16>(
+ multiply_add(sin_y, cos_k, multiply_add(cosm1_y, sin_k, sin_k)) /
+ multiply_add(sin_y, -sin_k, multiply_add(cosm1_y, cos_k, cos_k)));
+}
+
+} // namespace LIBC_NAMESPACE_DECL
diff --git a/libc/src/math/generic/tanpif16.cpp b/libc/src/math/generic/tanpif16.cpp
index 67635536ee3193..cf4f9917d45375 100644
--- a/libc/src/math/generic/tanpif16.cpp
+++ b/libc/src/math/generic/tanpif16.cpp
@@ -79,7 +79,7 @@ LLVM_LIBC_FUNCTION(float16, tanpif16, (float16 x)) {
// k = round(x * 32)
// y = x * 32 - k
//
- // Once k and y are computed, we then deduce the answer by tthe formula:
+ // Once k and y are computed, we then deduce the answer by the formula:
// tan(x) = sin(x) / cos(x)
// = (sin_y * cos_k + cos_y * sin_k) / (cos_y * cos_k - sin_y * sin_k)
float xf = x;
diff --git a/libc/src/math/tanf16.h b/libc/src/math/tanf16.h
new file mode 100644
index 00000000000000..bf1b61e9837f72
--- /dev/null
+++ b/libc/src/math/tanf16.h
@@ -0,0 +1,21 @@
+//===-- Implementation header for tanf16 ------------------------*- 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_TANF16_H
+#define LLVM_LIBC_SRC_MATH_TANF16_H
+
+#include "src/__support/macros/config.h"
+#include "src/__support/macros/properties/types.h"
+
+namespace LIBC_NAMESPACE_DECL {
+
+float16 tanf16(float16 x);
+
+} // namespace LIBC_NAMESPACE_DECL
+
+#endif // LLVM_LIBC_SRC_MATH_TANF16_H
diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt
index 16e7d4957ba114..ae8518ee4b4cc1 100644
--- a/libc/test/src/math/CMakeLists.txt
+++ b/libc/test/src/math/CMakeLists.txt
@@ -190,6 +190,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ tanf16_test
+ NEED_MPFR
+ SUITE
+ libc-math-unittests
+ SRCS
+ tanf16_test.cpp
+ DEPENDS
+ libc.src.math.tanf16
+)
+
add_fp_unittest(
tanpif16_test
NEED_MPFR
diff --git a/libc/test/src/math/cosf16_test.cpp b/libc/test/src/math/cosf16_test.cpp
index 9e4687f0325c49..b744e7817e4ba9 100644
--- a/libc/test/src/math/cosf16_test.cpp
+++ b/libc/test/src/math/cosf16_test.cpp
@@ -17,7 +17,7 @@ namespace mpfr = LIBC_NAMESPACE::testing::mpfr;
// Range: [0, Inf]
static constexpr uint16_t POS_START = 0x0000U;
-static constexpr uint16_t POS_STOP = 0x7c00u;
+static constexpr uint16_t POS_STOP = 0x7c00U;
// Range: [-Inf, 0]
static constexpr uint16_t NEG_START = 0x8000U;
diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt
index 31f85a3ecfd27b..e23e7f41222d4a 100644
--- a/libc/test/src/math/smoke/CMakeLists.txt
+++ b/libc/test/src/math/smoke/CMakeLists.txt
@@ -121,6 +121,17 @@ add_fp_unittest(
libc.src.__support.FPUtil.fp_bits
)
+add_fp_unittest(
+ tanf16_test
+ SUITE
+ libc-math-smoke-tests
+ SRCS
+ tanf16_test.cpp
+ DEPENDS
+ libc.src.errno.errno
+ libc.src.math.tanf16
+)
+
add_fp_unittest(
tanpif16_test
SUITE
diff --git a/libc/test/src/math/smoke/tanf16_test.cpp b/libc/test/src/math/smoke/tanf16_test.cpp
new file mode 100644
index 00000000000000..39d1182ba891e5
--- /dev/null
+++ b/libc/test/src/math/smoke/tanf16_test.cpp
@@ -0,0 +1,34 @@
+//===-- Unittests for tanf16 ----------------------------------------------===//
+//
+//
+// 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/tanf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+
+using LlvmLibcTanf16Test = LIBC_NAMESPACE::testing::FPTest<float16>;
+
+TEST_F(LlvmLibcTanf16Test, SpecialNumbers) {
+ LIBC_NAMESPACE::libc_errno = 0;
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::tanf16(aNaN));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(zero, LIBC_NAMESPACE::tanf16(zero));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::tanf16(neg_zero));
+ EXPECT_MATH_ERRNO(0);
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::tanf16(inf));
+ EXPECT_MATH_ERRNO(EDOM);
+
+ EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::tanf16(neg_inf));
+ EXPECT_MATH_ERRNO(EDOM);
+}
diff --git a/libc/test/src/math/tanf16_test.cpp b/libc/test/src/math/tanf16_test.cpp
new file mode 100644
index 00000000000000..f2e874182efc1d
--- /dev/null
+++ b/libc/test/src/math/tanf16_test.cpp
@@ -0,0 +1,40 @@
+//===-- Exhaustive test for tanf16 ----------------------------------------===//
+//
+// 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/tanf16.h"
+#include "test/UnitTest/FPMatcher.h"
+#include "test/UnitTest/Test.h"
+#include "utils/MPFRWrapper/MPFRUtils.h"
+
+using LlvmLibcTanf16Test = 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(LlvmLibcTanf16Test, PositiveRange) {
+ for (uint16_t v = POS_START; v <= POS_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tan, x,
+ LIBC_NAMESPACE::tanf16(x), 0.5);
+ }
+}
+
+TEST_F(LlvmLibcTanf16Test, NegativeRange) {
+ for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) {
+ float16 x = FPBits(v).get_val();
+ EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Tan, x,
+ LIBC_NAMESPACE::tanf16(x), 0.5);
+ }
+}
>From 1111a4ed3d512c356a527f447ce70a81884de808 Mon Sep 17 00:00:00 2001
From: wldfngrs <wldfngrs at gmail.com>
Date: Tue, 31 Dec 2024 13:03:29 +0100
Subject: [PATCH 2/2] Add comments, update docs
---
libc/docs/headers/math/index.rst | 2 +-
libc/src/math/generic/sincosf16_utils.h | 20 +++++++++-----------
libc/src/math/generic/tanf16.cpp | 11 +++++++----
3 files changed, 17 insertions(+), 16 deletions(-)
diff --git a/libc/docs/headers/math/index.rst b/libc/docs/headers/math/index.rst
index 2808165ad539b6..8548e4a5773bc0 100644
--- a/libc/docs/headers/math/index.rst
+++ b/libc/docs/headers/math/index.rst
@@ -346,7 +346,7 @@ Higher Math Functions
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| sqrt | |check| | |check| | |check| | |check| | |check| | 7.12.7.10 | F.10.4.10 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
-| tan | |check| | |check| | | | | 7.12.4.7 | F.10.1.7 |
+| tan | |check| | |check| | | |check| | | 7.12.4.7 | F.10.1.7 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
| tanh | |check| | | | |check| | | 7.12.5.6 | F.10.2.6 |
+-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+
diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h
index 48feee590937fd..79fe901b66d41d 100644
--- a/libc/src/math/generic/sincosf16_utils.h
+++ b/libc/src/math/generic/sincosf16_utils.h
@@ -49,22 +49,20 @@ LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) {
// Recall, range reduction:
// k = round(x * 32/pi)
-// y = x * 32/pi - k
//
-// The constant 0x1.45f306dc9c883p3 is 32/pi rounded to double-precision.
-// 32/pi is generated by Sollya with the following commands:
-// > display = hexadecimal;
-// > round(32/pi, D, RN);
-//
-// The precision choice of 'double' is to minimize rounding errors
-// in this initial scaling step, preserving enough bits so errors accumulated
-// while computing the subtraction: y = x * 32/pi - round(x * 32/pi)
+// The precision choice of 'double' in the following function is to minimize
+// rounding errors in this initial scaling step,
+// preserving enough bits so errors accumulated while computing the subtraction:
+// y = x * 32/pi - round(x * 32/pi)
// are beyond the least-significant bit of single-precision used during
// further intermediate computation.
LIBC_INLINE int32_t range_reduction_sincosf16(float x, float &y) {
- double prod = x * 0x1.45f306dc9c883p3;
- double kd = fputil::nearest_integer(prod);
+ // Generated by Sollya with:
+ // > D(32/pi);
+ constexpr double THIRTYTWO_OVER_PI = 0x1.45f306dc9c883p3;
+ double prod = x * THIRTYTWO_OVER_PI;
+ double kd = fputil::nearest_integer(prod);
y = static_cast<float>(prod - kd);
return static_cast<int32_t>(kd);
diff --git a/libc/src/math/generic/tanf16.cpp b/libc/src/math/generic/tanf16.cpp
index bae77f9f331871..b88651edef45de 100644
--- a/libc/src/math/generic/tanf16.cpp
+++ b/libc/src/math/generic/tanf16.cpp
@@ -50,9 +50,10 @@ LLVM_LIBC_FUNCTION(float16, tanf16, (float16 x)) {
// |x| <= 0x1.d1p-5
if (LIBC_UNLIKELY(x_abs <= 0x2b44)) {
+ // |x| <= 1.398p-11
if (LIBC_UNLIKELY(x_abs <= 0x10e6)) {
// tan(+/-0) = +/-0
- if (LIBC_UNLIKELY(x_abs == 0U))
+ if (LIBC_UNLIKELY(x_abs == 0))
return x;
int rounding = fputil::quick_get_round();
@@ -64,20 +65,22 @@ LLVM_LIBC_FUNCTION(float16, tanf16, (float16 x)) {
if ((xbits.is_pos() && rounding == FE_UPWARD) ||
(xbits.is_neg() && rounding == FE_DOWNWARD))
return fputil::cast<float16>(fputil::multiply_add(xf, 0x1.0p-11f, xf));
- else
- return x;
+ return x;
}
float xsq = xf * xf;
+ // Degree-6 minimax odd polynomial of tan(x) generated by Sollya with:
+ // > P = fpminimax(tan(x)/x, [|0, 2, 4, 6|], [|1, SG...|], [0, pi/32]);
float result = fputil::polyeval(xsq, 0x1p0f, 0x1.555556p-2f, 0x1.110ee4p-3f,
0x1.be80f6p-5f);
return fputil::cast<float16>(xf * result);
}
- // tan(+/-inf)= NaN, and tan(NaN) = NaN
+ // tan(+/-inf) = NaN, and tan(NaN) = NaN
if (LIBC_UNLIKELY(x_abs >= 0x7c00)) {
+ // x = +/-inf
if (x_abs == 0x7c00) {
fputil::set_errno_if_required(EDOM);
fputil::raise_except_if_required(FE_INVALID);
More information about the libc-commits
mailing list