diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index af7f429561fe0..5e9cc71279ab1 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -700,6 +700,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/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| | | | | | | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ 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 a5d17ad023f52..aeb758d4a092d 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -498,6 +498,27 @@ add_entrypoint_object( ${libc_opt_high_flag} ) +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.except_value_utils + libc.src.__support.FPUtil.multiply_add + libc.src.__support.macros.optimization + libc.src.__support.macros.properties.types + 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..5e5edd4a8c85b 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" @@ -46,10 +47,31 @@ 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); +// 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 +// 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(prod - kf); + + 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 +94,22 @@ 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..86546348ba739 --- /dev/null +++ b/libc/src/math/generic/sinf16.cpp @@ -0,0 +1,108 @@ +//===-- 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. +// 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/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 = 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 = 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) + + // Handle exceptional values + 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(); + } + + 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)) { + // 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())) + 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(); + } + } + + if (xbits.is_inf_or_nan()) { + if (xbits.is_inf()) { + 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/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/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..b05501cb0f145 --- /dev/null +++ b/libc/test/src/math/sinf16_test.cpp @@ -0,0 +1,40 @@ +//===-- 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..2966c3c952fd2 --- /dev/null +++ b/libc/test/src/math/smoke/sinf16_test.cpp @@ -0,0 +1,33 @@ +//===-- 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); +}