From 49b868b06308000df1c8bf50ccf41ea0d4c4dcfb Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 18 Nov 2024 18:57:45 +0100 Subject: [PATCH 01/12] Add sinf16 function --- libc/config/linux/x86_64/entrypoints.txt | 1 + libc/newhdrgen/yaml/math.yaml | 7 +++ libc/src/math/CMakeLists.txt | 1 + libc/src/math/generic/CMakeLists.txt | 19 ++++++ libc/src/math/generic/sincosf16_utils.h | 35 +++++++++-- libc/src/math/generic/sinf16.cpp | 77 ++++++++++++++++++++++++ libc/src/math/sinf16.h | 21 +++++++ libc/test/src/math/CMakeLists.txt | 11 ++++ libc/test/src/math/sinf16_test.cpp | 38 ++++++++++++ libc/test/src/math/smoke/CMakeLists.txt | 11 ++++ libc/test/src/math/smoke/sinf16_test.cpp | 34 +++++++++++ 11 files changed, 251 insertions(+), 4 deletions(-) create mode 100644 libc/src/math/generic/sinf16.cpp create mode 100644 libc/src/math/sinf16.h create mode 100644 libc/test/src/math/sinf16_test.cpp create mode 100644 libc/test/src/math/smoke/sinf16_test.cpp diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index 957e28bd66cc4..99f4e34262a8d 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -699,6 +699,7 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.scalbnf16 libc.src.math.setpayloadf16 libc.src.math.setpayloadsigf16 + libc.src.math.sinf16 libc.src.math.sinhf16 libc.src.math.sinpif16 libc.src.math.sqrtf16 diff --git a/libc/newhdrgen/yaml/math.yaml b/libc/newhdrgen/yaml/math.yaml index e09f0929e45f8..00efc34789667 100644 --- a/libc/newhdrgen/yaml/math.yaml +++ b/libc/newhdrgen/yaml/math.yaml @@ -2339,6 +2339,13 @@ functions: return_type: float arguments: - type: float + - name: sinf16 + standards: + - stdc + return_type: _Float16 + arguments: + - type: _Float16 + guard: LIBC_TYPES_HAS_FLOAT16 - name: sinhf standards: - stdc diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index 76a5e491effa0..390a59d07a28b 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -484,6 +484,7 @@ add_math_entrypoint_object(sincosf) add_math_entrypoint_object(sin) add_math_entrypoint_object(sinf) +add_math_entrypoint_object(sinf16) add_math_entrypoint_object(sinpif) add_math_entrypoint_object(sinpif16) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 34b0f1424e8fd..51396865b99d4 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -496,6 +496,25 @@ add_entrypoint_object( -O3 ) +add_entrypoint_object( + sinf16 + SRCS + sinf16.cpp + HDRS + ../sinf16.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.multiply_add + libc.src.__support.macros.optimization + COMPILE_OPTIONS + -O3 +) + add_entrypoint_object( sincos SRCS diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h index 83511755a56c4..fe6e58c1f919c 100644 --- a/libc/src/math/generic/sincosf16_utils.h +++ b/libc/src/math/generic/sincosf16_utils.h @@ -39,6 +39,11 @@ constexpr float SIN_K_PI_OVER_32[64] = { -0x1.6a09e6p-1, -0x1.44cf32p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2, -0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4}; +static constexpr float THIRTYTWO_OVER_PI[4] = { +//0x1.45fp3, 0x1.837p-12, -0x1.b1b8p-28, -0x1.f568p-43 + 0x1.46p3, -0x1.9f4p-10, 0x1.b94p-22, -0x1.bcp-36 +}; + LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) { float kf = fputil::nearest_integer(x * 32); y = fputil::multiply_add(x, 32.0, -kf); @@ -46,11 +51,18 @@ LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) { return static_cast(kf); } -LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, - float &sin_y, float &cosm1_y) { - float y; - int32_t k = range_reduction_sincospif16(xf, y); +LIBC_INLINE int32_t range_reduction_sincosf16(float x, float &y) { + double xd = x; + double prod = xd * 0x1.45f306dc9c883p3; + double kf = fputil::nearest_integer(prod); + double yd = prod - kf; + y = static_cast(yd); + return static_cast(kf); +} + +static LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k, float &cos_k, float &sin_y, float &cosm1_y) { + sin_k = SIN_K_PI_OVER_32[k & 63]; cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; @@ -72,6 +84,21 @@ LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, 0x1.a6f7a2p-29f); } +LIBC_INLINE void sincosf16_eval(float xf, float &sin_k, float &cos_k, float &sin_y, float &cosm1_y) { + float y; + int32_t k = range_reduction_sincosf16(xf, y); + + sincosf16_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y); +} + +LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, + float &sin_y, float &cosm1_y) { + float y; + int32_t k = range_reduction_sincospif16(xf, y); + + sincosf16_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y); +} + } // namespace LIBC_NAMESPACE_DECL #endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp new file mode 100644 index 0000000000000..7f1c8c7cf6521 --- /dev/null +++ b/libc/src/math/generic/sinf16.cpp @@ -0,0 +1,77 @@ +//===-- Half-precision sin 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/sinf16.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/multiply_add.h" +#include "src/__support/macros/optimization.h" + +namespace LIBC_NAMESPACE_DECL { + +LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + + uint16_t x_u = xbits.uintval(); + uint16_t x_abs = x_u & 0x7fff; + float xf = x; + + // 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 sine of sum + // formula: + // sin(x) = sin((k + y) * pi/32) + // = sin(k * pi/32) * cos(y * pi/32) + + // sin(y * pi/32) * cos(k * pi/32) + + if (LIBC_UNLIKELY(x_abs <= 0x13d0)) { + int rounding = fputil::quick_get_round(); + + // For signed zeros + if ((LIBC_UNLIKELY(x_abs == 0U)) || + (rounding == FE_UPWARD && xbits.is_pos()) || + (rounding == FE_DOWNWARD && xbits.is_neg())) + return x; + } + + if (LIBC_UNLIKELY(x_abs >= 0x7c00)) { + // If value is equal to infinity + if (x_abs == 0x7c00) { + fputil::set_errno_if_required(EDOM); + fputil::raise_except_if_required(FE_INVALID); + } + + return x + FPBits::quiet_nan().get_val(); + } + + + float sin_k, cos_k, sin_y, cosm1_y; + sincosf16_eval(xf, sin_k, cos_k, sin_y, cosm1_y); + + if (LIBC_UNLIKELY(sin_y == 0 && sin_k == 0)) + return FPBits::zero(xbits.sign()).get_val(); + + // Since, cosm1_y = cos_y - 1, therfore: + // sin(x) = cos_k * sin_y + sin_k + (cosm1_y * sin_k) + return fputil::cast(fputil::multiply_add(sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k))); +} + +} // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/sinf16.h b/libc/src/math/sinf16.h new file mode 100644 index 0000000000000..23f1aa99b6233 --- /dev/null +++ b/libc/src/math/sinf16.h @@ -0,0 +1,21 @@ +//===-- Implementation header for sinf16 ------------------------*- 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_SINF16_H +#define LLVM_LIBC_SRC_MATH_SINF16_H + +#include "src/__support/macros/config.h" +#include "src/__support/macros/properties/types.h" + +namespace LIBC_NAMESPACE_DECL { + +float16 sinf16(float16 x); + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC_MATH_SINF16_H diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index 610f4d9fc1a3b..ea75720df4f43 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -85,6 +85,17 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + sinf16_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + sinf16_test.cpp + DEPENDS + libc.src.math.sinf16 +) + add_fp_unittest( sinpif_test NEED_MPFR diff --git a/libc/test/src/math/sinf16_test.cpp b/libc/test/src/math/sinf16_test.cpp new file mode 100644 index 0000000000000..cdb69d388f7f7 --- /dev/null +++ b/libc/test/src/math/sinf16_test.cpp @@ -0,0 +1,38 @@ +//===-- Exhaustive test for sinf16 ----------------------------------------===// +// +// 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/sinf16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +using LlvmLibcSinf16Test = 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(LlvmLibcSinf16Test, PositiveRange) { + for (uint16_t v = POS_START; v <= POS_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, LIBC_NAMESPACE::sinf16(x), 0.5); + } +} + +TEST_F(LlvmLibcSinf16Test, NegativeRange) { + for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, LIBC_NAMESPACE::sinf16(x), 0.5); + } +} diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index e9c785f7d9330..2c1c4dba73846 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -49,6 +49,17 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + sinf16_test + SUITE + libc-math-smoke-tests + SRCS + sinf16_test.cpp + DEPENDS + libc.src.errno.errno + libc.src.math.sinf16 +) + add_fp_unittest( sinpif_test SUITE diff --git a/libc/test/src/math/smoke/sinf16_test.cpp b/libc/test/src/math/smoke/sinf16_test.cpp new file mode 100644 index 0000000000000..301e91a683a3e --- /dev/null +++ b/libc/test/src/math/smoke/sinf16_test.cpp @@ -0,0 +1,34 @@ +//===-- Unittests for sinf16 ----------------------------------------------===// +// +// 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/sinf16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" + +using LlvmLibcSinf16Test = LIBC_NAMESPACE::testing::FPTest; + +TEST_F(LlvmLibcSinf16Test, SpecialNumbers) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinf16(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinf16(zero)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinf16(neg_zero)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinf16(inf)); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinf16(neg_inf)); + EXPECT_MATH_ERRNO(EDOM); +} From 1d5f809cf8f22a33d01038c407f9fa800bb9b05c Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Fri, 22 Nov 2024 20:19:13 +0100 Subject: [PATCH 02/12] Fixed rounding errors --- libc/src/math/generic/CMakeLists.txt | 1 + libc/src/math/generic/sincosf16_utils.h | 23 +++++------ libc/src/math/generic/sinf16.cpp | 49 +++++++++++++++++++----- libc/src/math/generic/tanpif16.cpp | 4 +- libc/test/src/math/sinf16_test.cpp | 6 ++- libc/test/src/math/smoke/sinf16_test.cpp | 3 +- 6 files changed, 57 insertions(+), 29 deletions(-) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 51396865b99d4..97a6784ef3c5f 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -509,6 +509,7 @@ add_entrypoint_object( 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 COMPILE_OPTIONS diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h index fe6e58c1f919c..13d27fb9828a3 100644 --- a/libc/src/math/generic/sincosf16_utils.h +++ b/libc/src/math/generic/sincosf16_utils.h @@ -11,6 +11,7 @@ #include "src/__support/FPUtil/FPBits.h" #include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/cast.h" #include "src/__support/FPUtil/nearest_integer.h" #include "src/__support/common.h" #include "src/__support/macros/config.h" @@ -39,11 +40,6 @@ constexpr float SIN_K_PI_OVER_32[64] = { -0x1.6a09e6p-1, -0x1.44cf32p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2, -0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4}; -static constexpr float THIRTYTWO_OVER_PI[4] = { -//0x1.45fp3, 0x1.837p-12, -0x1.b1b8p-28, -0x1.f568p-43 - 0x1.46p3, -0x1.9f4p-10, 0x1.b94p-22, -0x1.bcp-36 -}; - LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) { float kf = fputil::nearest_integer(x * 32); y = fputil::multiply_add(x, 32.0, -kf); @@ -52,17 +48,17 @@ LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) { } LIBC_INLINE int32_t range_reduction_sincosf16(float x, float &y) { - double xd = x; - double prod = xd * 0x1.45f306dc9c883p3; + double prod = x * 0x1.45f306dc9c883p3; double kf = fputil::nearest_integer(prod); - double yd = prod - kf; + y = static_cast(prod - kf); - y = static_cast(yd); return static_cast(kf); } -static LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k, float &cos_k, float &sin_y, float &cosm1_y) { - +static LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k, + float &cos_k, float &sin_y, + float &cosm1_y) { + sin_k = SIN_K_PI_OVER_32[k & 63]; cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; @@ -84,7 +80,8 @@ static LIBC_INLINE void sincosf16_poly_eval(int32_t k, float y, float &sin_k, fl 0x1.a6f7a2p-29f); } -LIBC_INLINE void sincosf16_eval(float xf, float &sin_k, float &cos_k, float &sin_y, float &cosm1_y) { +LIBC_INLINE void sincosf16_eval(float xf, float &sin_k, float &cos_k, + float &sin_y, float &cosm1_y) { float y; int32_t k = range_reduction_sincosf16(xf, y); @@ -95,7 +92,7 @@ LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, float &sin_y, float &cosm1_y) { float y; int32_t k = range_reduction_sincospif16(xf, y); - + sincosf16_poly_eval(k, y, sin_k, cos_k, sin_y, cosm1_y); } diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp index 7f1c8c7cf6521..f307b3ea330c5 100644 --- a/libc/src/math/generic/sinf16.cpp +++ b/libc/src/math/generic/sinf16.cpp @@ -13,11 +13,20 @@ #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 = 2; + +constexpr fputil::ExceptValues SINF16_EXCEPTS{{ + // (input, RZ output, RU offset, RD offset, RN offset) + {0x585c, 0x3ba3, 1, 0, 1}, + {0x5cb0, 0xbbff, 0, 1, 0}, +}}; + LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { using FPBits = typename fputil::FPBits; FPBits xbits(x); @@ -25,7 +34,7 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { uint16_t x_u = xbits.uintval(); uint16_t x_abs = x_u & 0x7fff; float xf = x; - + // Range reduction: // For !x| > pi/32, we perform range reduction as follows: // Find k and y such that: @@ -41,15 +50,35 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { // sin(x) = sin((k + y) * pi/32) // = sin(k * pi/32) * cos(y * pi/32) + // sin(y * pi/32) * cos(k * pi/32) - + + // Handle exceptional values + if (LIBC_UNLIKELY(x_abs == 0x585C || x_abs == 0x5cb0)) { + bool x_sign = x_u >> 15; + if (auto r = SINF16_EXCEPTS.lookup_odd(x_abs, x_sign); + LIBC_UNLIKELY(r.has_value())) + return r.value(); + } + + int rounding = fputil::quick_get_round(); if (LIBC_UNLIKELY(x_abs <= 0x13d0)) { - int rounding = fputil::quick_get_round(); - - // For signed zeros - if ((LIBC_UNLIKELY(x_abs == 0U)) || - (rounding == FE_UPWARD && xbits.is_pos()) || - (rounding == FE_DOWNWARD && xbits.is_neg())) + if (LIBC_UNLIKELY(x_abs == 0U)) return x; + + if ((rounding == FE_UPWARD && xbits.is_pos()) || + (rounding == FE_DOWNWARD && xbits.is_neg())) + return x; + + if (rounding == FE_UPWARD && xbits.is_neg()) { + x_u--; + return FPBits(x_u).get_val(); + } + } + + if (LIBC_UNLIKELY(x_abs == 0x2b45)) { + if (rounding == FE_DOWNWARD && xbits.is_neg()) { + x_u--; + return FPBits(x_u).get_val(); + } } if (LIBC_UNLIKELY(x_abs >= 0x7c00)) { @@ -61,7 +90,6 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { return x + FPBits::quiet_nan().get_val(); } - float sin_k, cos_k, sin_y, cosm1_y; sincosf16_eval(xf, sin_k, cos_k, sin_y, cosm1_y); @@ -71,7 +99,8 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { // Since, cosm1_y = cos_y - 1, therfore: // sin(x) = cos_k * sin_y + sin_k + (cosm1_y * sin_k) - return fputil::cast(fputil::multiply_add(sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k))); + return fputil::cast(fputil::multiply_add( + sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k))); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/tanpif16.cpp b/libc/src/math/generic/tanpif16.cpp index ab3c9cb2122ba..67635536ee319 100644 --- a/libc/src/math/generic/tanpif16.cpp +++ b/libc/src/math/generic/tanpif16.cpp @@ -21,7 +21,7 @@ namespace LIBC_NAMESPACE_DECL { constexpr size_t N_EXCEPTS = 21; -constexpr fputil::ExceptValues TANF16_EXCEPTS{{ +constexpr fputil::ExceptValues TANPIF16_EXCEPTS{{ // (input, RZ output, RU offset, RD offset, RN offset) {0x07f2, 0x0e3d, 1, 0, 0}, {0x086a, 0x0eee, 1, 0, 1}, {0x08db, 0x0fa0, 1, 0, 0}, {0x094c, 0x1029, 1, 0, 0}, @@ -49,7 +49,7 @@ LLVM_LIBC_FUNCTION(float16, tanpif16, (float16 x)) { return x; bool x_sign = x_u >> 15; - if (auto r = TANF16_EXCEPTS.lookup_odd(x_abs, x_sign); + if (auto r = TANPIF16_EXCEPTS.lookup_odd(x_abs, x_sign); LIBC_UNLIKELY(r.has_value())) return r.value(); } diff --git a/libc/test/src/math/sinf16_test.cpp b/libc/test/src/math/sinf16_test.cpp index cdb69d388f7f7..b469e1dcb8d33 100644 --- a/libc/test/src/math/sinf16_test.cpp +++ b/libc/test/src/math/sinf16_test.cpp @@ -26,13 +26,15 @@ static constexpr uint16_t NEG_STOP = 0xfc00U; TEST_F(LlvmLibcSinf16Test, PositiveRange) { for (uint16_t v = POS_START; v <= POS_STOP; ++v) { float16 x = FPBits(v).get_val(); - EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, LIBC_NAMESPACE::sinf16(x), 0.5); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, + LIBC_NAMESPACE::sinf16(x), 0.5); } } TEST_F(LlvmLibcSinf16Test, NegativeRange) { for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) { float16 x = FPBits(v).get_val(); - EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, LIBC_NAMESPACE::sinf16(x), 0.5); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sin, x, + LIBC_NAMESPACE::sinf16(x), 0.5); } } diff --git a/libc/test/src/math/smoke/sinf16_test.cpp b/libc/test/src/math/smoke/sinf16_test.cpp index 301e91a683a3e..2966c3c952fd2 100644 --- a/libc/test/src/math/smoke/sinf16_test.cpp +++ b/libc/test/src/math/smoke/sinf16_test.cpp @@ -6,7 +6,6 @@ // //===----------------------------------------------------------------------===// - #include "src/errno/libc_errno.h" #include "src/math/sinf16.h" #include "test/UnitTest/FPMatcher.h" @@ -16,7 +15,7 @@ using LlvmLibcSinf16Test = LIBC_NAMESPACE::testing::FPTest; TEST_F(LlvmLibcSinf16Test, SpecialNumbers) { LIBC_NAMESPACE::libc_errno = 0; - + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinf16(aNaN)); EXPECT_MATH_ERRNO(0); From 09b79304b02c29a1cf5756a3765c7e97836ffdea Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Fri, 22 Nov 2024 21:29:20 +0100 Subject: [PATCH 03/12] small fix --- libc/src/math/generic/sinf16.cpp | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp index f307b3ea330c5..75d96659da054 100644 --- a/libc/src/math/generic/sinf16.cpp +++ b/libc/src/math/generic/sinf16.cpp @@ -19,12 +19,13 @@ namespace LIBC_NAMESPACE_DECL { -constexpr size_t N_EXCEPTS = 2; +constexpr size_t N_EXCEPTS = 3; constexpr fputil::ExceptValues SINF16_EXCEPTS{{ // (input, RZ output, RU offset, RD offset, RN offset) {0x585c, 0x3ba3, 1, 0, 1}, {0x5cb0, 0xbbff, 0, 1, 0}, + {0x51f5, 0xb80f, 0, 1, 0}, }}; LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { @@ -52,7 +53,7 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { // sin(y * pi/32) * cos(k * pi/32) // Handle exceptional values - if (LIBC_UNLIKELY(x_abs == 0x585C || x_abs == 0x5cb0)) { + if (LIBC_UNLIKELY(x_abs == 0x585C || x_abs == 0x5cb0 || x_abs == 0x51f5)) { bool x_sign = x_u >> 15; if (auto r = SINF16_EXCEPTS.lookup_odd(x_abs, x_sign); LIBC_UNLIKELY(r.has_value())) From 0eefcc29e067087ddce6d24323d559346e67e923 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Fri, 22 Nov 2024 21:46:04 +0100 Subject: [PATCH 04/12] Update Docs --- libc/docs/math/index.rst | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst index 2b86f49a3619e..4934e93ccb164 100644 --- a/libc/docs/math/index.rst +++ b/libc/docs/math/index.rst @@ -336,7 +336,7 @@ Higher Math Functions +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | rsqrt | | | | | | 7.12.7.9 | F.10.4.9 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| sin | |check| | |check| | | | | 7.12.4.6 | F.10.1.6 | +| sin | |check| | |check| | | |check| | | 7.12.4.6 | F.10.1.6 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | sincos | |check| | |check| | | | | | | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ From fd28fcd44ffbf3bcb0ab1247860816ea7fbeb775 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sun, 24 Nov 2024 00:55:08 +0100 Subject: [PATCH 05/12] nit formatting --- libc/src/math/generic/sinf16.cpp | 2 +- libc/test/src/math/sinf16_test.cpp | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp index 75d96659da054..dcee9be748770 100644 --- a/libc/src/math/generic/sinf16.cpp +++ b/libc/src/math/generic/sinf16.cpp @@ -1,4 +1,4 @@ -//===-- Half-precision sin function ------------------------------------===// +//===-- Half-precision sin(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. diff --git a/libc/test/src/math/sinf16_test.cpp b/libc/test/src/math/sinf16_test.cpp index b469e1dcb8d33..b05501cb0f145 100644 --- a/libc/test/src/math/sinf16_test.cpp +++ b/libc/test/src/math/sinf16_test.cpp @@ -4,7 +4,7 @@ // See https://llvm.org/LICENSE.txt for license information. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // -//===---------------------------------------------------------------------===// +//===----------------------------------------------------------------------===// #include "src/math/sinf16.h" #include "test/UnitTest/FPMatcher.h" From e234abd17bf426f228caf46bfafc5369bc821cf7 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sun, 24 Nov 2024 00:57:02 +0100 Subject: [PATCH 06/12] formatting --- libc/src/math/generic/sinf16.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp index dcee9be748770..bb36b20c69658 100644 --- a/libc/src/math/generic/sinf16.cpp +++ b/libc/src/math/generic/sinf16.cpp @@ -1,4 +1,4 @@ -//===-- Half-precision sin(x) function -------------------------------------===// +//===-- Half-precision sin(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. From c6bc8571fffb30a5f5f3d31d09f1ffe0f612beb4 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Wed, 27 Nov 2024 23:45:02 +0100 Subject: [PATCH 07/12] Added comments, improved code readability --- libc/src/math/generic/sincosf16_utils.h | 14 +++++++++++++ libc/src/math/generic/sinf16.cpp | 27 ++++++++++++------------- 2 files changed, 27 insertions(+), 14 deletions(-) diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h index 13d27fb9828a3..b3d262519e1ad 100644 --- a/libc/src/math/generic/sincosf16_utils.h +++ b/libc/src/math/generic/sincosf16_utils.h @@ -47,6 +47,20 @@ LIBC_INLINE int32_t range_reduction_sincospif16(float x, float &y) { return static_cast(kf); } +// 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) +// are beyond the least-significant bit of single-precision used during +// intermediate computation. LIBC_INLINE int32_t range_reduction_sincosf16(float x, float &y) { double prod = x * 0x1.45f306dc9c883p3; double kf = fputil::nearest_integer(prod); diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp index bb36b20c69658..7bc7daed2eb73 100644 --- a/libc/src/math/generic/sinf16.cpp +++ b/libc/src/math/generic/sinf16.cpp @@ -19,17 +19,18 @@ namespace LIBC_NAMESPACE_DECL { -constexpr size_t N_EXCEPTS = 3; +constexpr size_t N_EXCEPTS = 4; constexpr fputil::ExceptValues SINF16_EXCEPTS{{ // (input, RZ output, RU offset, RD offset, RN offset) + {0x2b45, 0x2b43, 1, 0, 1}, {0x585c, 0x3ba3, 1, 0, 1}, {0x5cb0, 0xbbff, 0, 1, 0}, {0x51f5, 0xb80f, 0, 1, 0}, }}; LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { - using FPBits = typename fputil::FPBits; + using FPBits = fputil::FPBits; FPBits xbits(x); uint16_t x_u = xbits.uintval(); @@ -37,7 +38,7 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { float xf = x; // Range reduction: - // For !x| > pi/32, we perform range reduction as follows: + // 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 @@ -53,13 +54,19 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { // sin(y * pi/32) * cos(k * pi/32) // Handle exceptional values - if (LIBC_UNLIKELY(x_abs == 0x585C || x_abs == 0x5cb0 || x_abs == 0x51f5)) { + if (LIBC_UNLIKELY(x_abs == 0x585c || x_abs == 0x5cb0 || x_abs == 0x51f5 || + x_abs == 0x2b45)) { bool x_sign = x_u >> 15; if (auto r = SINF16_EXCEPTS.lookup_odd(x_abs, x_sign); LIBC_UNLIKELY(r.has_value())) return r.value(); } + // Exhaustive tests show that for |x| <= 0x13d0, 1ULP rounding errors occur. + // To fix this, the following apply: + // - When x >= 0, and rounding upward, sin(x) == x. + // - When x < 0, and rounding downward, sin(x) == x. + // - When x < 0, and rounding upward, sin(x) == (x - 1ULP) int rounding = fputil::quick_get_round(); if (LIBC_UNLIKELY(x_abs <= 0x13d0)) { if (LIBC_UNLIKELY(x_abs == 0U)) @@ -75,16 +82,8 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { } } - if (LIBC_UNLIKELY(x_abs == 0x2b45)) { - if (rounding == FE_DOWNWARD && xbits.is_neg()) { - x_u--; - return FPBits(x_u).get_val(); - } - } - - if (LIBC_UNLIKELY(x_abs >= 0x7c00)) { - // If value is equal to infinity - if (x_abs == 0x7c00) { + if (xbits.is_inf_or_nan()) { + if (xbits.is_inf()){ fputil::set_errno_if_required(EDOM); fputil::raise_except_if_required(FE_INVALID); } From 92bdef59839511be287ffe1e3a236eac25fa0cf6 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Wed, 27 Nov 2024 23:47:43 +0100 Subject: [PATCH 08/12] clang-format --- libc/src/math/generic/sincosf16_utils.h | 8 ++++---- libc/src/math/generic/sinf16.cpp | 6 +++--- 2 files changed, 7 insertions(+), 7 deletions(-) diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h index b3d262519e1ad..5e5edd4a8c85b 100644 --- a/libc/src/math/generic/sincosf16_utils.h +++ b/libc/src/math/generic/sincosf16_utils.h @@ -50,17 +50,17 @@ 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) -// are beyond the least-significant bit of single-precision used during -// intermediate computation. +// 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 kf = fputil::nearest_integer(prod); diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp index 7bc7daed2eb73..f8c90264f17b6 100644 --- a/libc/src/math/generic/sinf16.cpp +++ b/libc/src/math/generic/sinf16.cpp @@ -55,7 +55,7 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { // Handle exceptional values if (LIBC_UNLIKELY(x_abs == 0x585c || x_abs == 0x5cb0 || x_abs == 0x51f5 || - x_abs == 0x2b45)) { + x_abs == 0x2b45)) { bool x_sign = x_u >> 15; if (auto r = SINF16_EXCEPTS.lookup_odd(x_abs, x_sign); LIBC_UNLIKELY(r.has_value())) @@ -63,7 +63,7 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { } // Exhaustive tests show that for |x| <= 0x13d0, 1ULP rounding errors occur. - // To fix this, the following apply: + // To fix this, the following apply: // - When x >= 0, and rounding upward, sin(x) == x. // - When x < 0, and rounding downward, sin(x) == x. // - When x < 0, and rounding upward, sin(x) == (x - 1ULP) @@ -83,7 +83,7 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { } if (xbits.is_inf_or_nan()) { - if (xbits.is_inf()){ + if (xbits.is_inf()) { fputil::set_errno_if_required(EDOM); fputil::raise_except_if_required(FE_INVALID); } From 7a72ef9e7ab58fefa5f6a99000a9db62488748dc Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Tue, 3 Dec 2024 01:33:13 +0100 Subject: [PATCH 09/12] minor change --- libc/src/math/generic/sinf16.cpp | 11 ++++++----- 1 file changed, 6 insertions(+), 5 deletions(-) diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp index f8c90264f17b6..6be352a8b6adb 100644 --- a/libc/src/math/generic/sinf16.cpp +++ b/libc/src/math/generic/sinf16.cpp @@ -62,20 +62,21 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { return r.value(); } - // Exhaustive tests show that for |x| <= 0x13d0, 1ULP rounding errors occur. - // To fix this, the following apply: - // - When x >= 0, and rounding upward, sin(x) == x. - // - When x < 0, and rounding downward, sin(x) == x. - // - When x < 0, and rounding upward, sin(x) == (x - 1ULP) int rounding = fputil::quick_get_round(); + + // Exhaustive tests show that for |x| <= 0x1.f4p-11, 1ULP rounding errors occur. + // To fix this, the following apply: if (LIBC_UNLIKELY(x_abs <= 0x13d0)) { + // When x >= 0, and rounding upward, sin(x) == x. if (LIBC_UNLIKELY(x_abs == 0U)) return x; + // When x < 0, and rounding downward, sin(x) == x. if ((rounding == FE_UPWARD && xbits.is_pos()) || (rounding == FE_DOWNWARD && xbits.is_neg())) return x; + // When x < 0, and rounding upward, sin(x) == (x - 1ULP) if (rounding == FE_UPWARD && xbits.is_neg()) { x_u--; return FPBits(x_u).get_val(); From 2a6e3a53f538e4fbcdec73db9b6552c8175c01aa Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Tue, 3 Dec 2024 01:35:12 +0100 Subject: [PATCH 10/12] clang format --- libc/src/math/generic/sinf16.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp index 6be352a8b6adb..6dadc37eeace5 100644 --- a/libc/src/math/generic/sinf16.cpp +++ b/libc/src/math/generic/sinf16.cpp @@ -63,9 +63,9 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { } int rounding = fputil::quick_get_round(); - - // Exhaustive tests show that for |x| <= 0x1.f4p-11, 1ULP rounding errors occur. - // To fix this, the following apply: + + // Exhaustive tests show that for |x| <= 0x1.f4p-11, 1ULP rounding errors + // occur. To fix this, the following apply: if (LIBC_UNLIKELY(x_abs <= 0x13d0)) { // When x >= 0, and rounding upward, sin(x) == x. if (LIBC_UNLIKELY(x_abs == 0U)) From ec08e2c9aeef046c50562c00ddb78ea747152d59 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Tue, 3 Dec 2024 13:02:43 +0100 Subject: [PATCH 11/12] Update libc/src/math/generic/sinf16.cpp Fixed comments Co-authored-by: OverMighty --- libc/src/math/generic/sinf16.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/libc/src/math/generic/sinf16.cpp b/libc/src/math/generic/sinf16.cpp index 6dadc37eeace5..86546348ba739 100644 --- a/libc/src/math/generic/sinf16.cpp +++ b/libc/src/math/generic/sinf16.cpp @@ -67,10 +67,11 @@ LLVM_LIBC_FUNCTION(float16, sinf16, (float16 x)) { // Exhaustive tests show that for |x| <= 0x1.f4p-11, 1ULP rounding errors // occur. To fix this, the following apply: if (LIBC_UNLIKELY(x_abs <= 0x13d0)) { - // When x >= 0, and rounding upward, sin(x) == x. + // sin(+/-0) = +/-0 if (LIBC_UNLIKELY(x_abs == 0U)) return x; + // When x > 0, and rounding upward, sin(x) == x. // When x < 0, and rounding downward, sin(x) == x. if ((rounding == FE_UPWARD && xbits.is_pos()) || (rounding == FE_DOWNWARD && xbits.is_neg())) From 9d21c3f1d45e0444a17700ddac4fb0fc46adcf4b Mon Sep 17 00:00:00 2001 From: OverMighty Date: Tue, 3 Dec 2024 20:09:16 +0100 Subject: [PATCH 12/12] Update libc/src/math/generic/CMakeLists.txt --- libc/src/math/generic/CMakeLists.txt | 1 + 1 file changed, 1 insertion(+) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 2dea88a1aa34e..aeb758d4a092d 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -514,6 +514,7 @@ add_entrypoint_object( libc.src.__support.FPUtil.except_value_utils libc.src.__support.FPUtil.multiply_add libc.src.__support.macros.optimization + libc.src.__support.macros.properties.types COMPILE_OPTIONS -O3 )