From 0d59443daf39e994757f372d9aeb30fb4c22418b Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 00:14:59 +0100 Subject: [PATCH 01/15] Add cospif16 function --- libc/config/linux/aarch64/entrypoints.txt | 1 + libc/config/linux/x86_64/entrypoints.txt | 1 + libc/newhdrgen/yaml/math.yaml | 7 + libc/src/math/CMakeLists.txt | 1 + libc/src/math/cospif16.h | 21 +++ libc/src/math/generic/CMakeLists.txt | 20 +++ libc/src/math/generic/cospif16.cpp | 141 +++++++++++++++++++++ libc/test/src/math/CMakeLists.txt | 11 ++ libc/test/src/math/cospif16_test.cpp | 38 ++++++ libc/test/src/math/smoke/CMakeLists.txt | 11 ++ libc/test/src/math/smoke/cospif16_test.cpp | 44 +++++++ 11 files changed, 296 insertions(+) create mode 100644 libc/src/math/cospif16.h create mode 100644 libc/src/math/generic/cospif16.cpp create mode 100644 libc/test/src/math/cospif16_test.cpp create mode 100644 libc/test/src/math/smoke/cospif16_test.cpp diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt index 885827d304efe..85bb5df358ec3 100644 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -607,6 +607,7 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.canonicalizef16 libc.src.math.ceilf16 libc.src.math.copysignf16 + libc.src.math.cospif16 # TODO: aarch64 bug # Please see https://github.com/llvm/llvm-project/pull/100632#issuecomment-2258772681 # libc.src.math.expf16 diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index 06ea7bba81f34..f40d752840b85 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -611,6 +611,7 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.ceilf16 libc.src.math.copysignf16 libc.src.math.coshf16 + libc.src.math.cospif16 libc.src.math.exp10f16 libc.src.math.exp10m1f16 libc.src.math.exp2f16 diff --git a/libc/newhdrgen/yaml/math.yaml b/libc/newhdrgen/yaml/math.yaml index 98ea1a0d25fbb..e0986f00a3b46 100644 --- a/libc/newhdrgen/yaml/math.yaml +++ b/libc/newhdrgen/yaml/math.yaml @@ -206,6 +206,13 @@ functions: return_type: float arguments: - type: float + - name: cospif16 + standards: + - stdc + return_type: _Float16 + arguments: + - type: _Float16 + guard: LIBC_TYPES_HAS_FLOAT16 - name: ddivl standards: - stdc diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index 516bed499b194..3836d6562a074 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -95,6 +95,7 @@ add_math_entrypoint_object(coshf) add_math_entrypoint_object(coshf16) add_math_entrypoint_object(cospif) +add_math_entrypoint_object(cospif16) add_math_entrypoint_object(daddl) add_math_entrypoint_object(daddf128) diff --git a/libc/src/math/cospif16.h b/libc/src/math/cospif16.h new file mode 100644 index 0000000000000..6779e67cdccae --- /dev/null +++ b/libc/src/math/cospif16.h @@ -0,0 +1,21 @@ +//===-- Implementation header for cospif16 ---------------------*- 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_COSPIF16_H +#define LLVM_LIBC_SRC_MATH_COSPIF16_H + +#include "src/__support/macros/config.h" +#include "src/__support/macros/properties/types.h" + +namespace LIBC_NAMESPACE_DECL { + +float16 cospif16(float16 x); + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC_MATH_SINPIF16_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index d7c7a3431d3d9..e535c9e50b9b8 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -422,6 +422,26 @@ add_entrypoint_object( -O3 ) + +add_entrypoint_object( + cospif16 + SRCS + cospif16.cpp + HDRS + ../cospif16.h + DEPENDS + libc.src.__support.common + 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.FPUtil.nearest_integer + libc.src.__support.FPUtil.polyeval + libc.src.__support.macros.properties.types + COMPILE_OPTIONS + -O3 +) + add_entrypoint_object( sin SRCS diff --git a/libc/src/math/generic/cospif16.cpp b/libc/src/math/generic/cospif16.cpp new file mode 100644 index 0000000000000..779aab781b9f0 --- /dev/null +++ b/libc/src/math/generic/cospif16.cpp @@ -0,0 +1,141 @@ +//===-- Half-precision cospif 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/cospif16.h" +#include "src/__support/FPUtil/FEnvImpl.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/cast.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/FPUtil/nearest_integer.h" +#include "src/__support/common.h" +#include "src/__support/macros/config.h" + +namespace LIBC_NAMESPACE_DECL { +// Lookup table for sin(k * pi / 32) with k = 0, ..., 63. +// Table is generated with Sollya as follows: +// > display = hexadecimal; +// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); }; +static constexpr float SIN_K_PI_OVER_32[64] = { + 0x0.0p0, 0x1.917a6cp-4, 0x1.8f8b84p-3, 0x1.294062p-2, + 0x1.87de2ap-2, 0x1.e2b5d4p-2, 0x1.1c73b4p-1, 0x1.44cf32p-1, + 0x1.6a09e6p-1, 0x1.8bc806p-1, 0x1.a9b662p-1, 0x1.c38b3p-1, + 0x1.d906bcp-1, 0x1.e9f416p-1, 0x1.f6297cp-1, 0x1.fd88dap-1, + 0x1p0, 0x1.fd88dap-1, 0x1.f6297cp-1, 0x1.e9f416p-1, + 0x1.d906bcp-1, 0x1.c38b3p-1, 0x1.a9b662p-1, 0x1.8bc806p-1, + 0x1.6a09e6p-1, 0x1.44cf32p-1, 0x1.1c73b4p-1, 0x1.e2b5d4p-2, + 0x1.87de2ap-2, 0x1.294062p-2, 0x1.8f8b84p-3, 0x1.917a6cp-4, + 0x0.0p0, -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294062p-2, + -0x1.87de2ap-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, -0x1.44cf32p-1, + -0x1.6a09e6p-1, -0x1.8bc806p-1, -0x1.a9b662p-1, -0x1.c38b3p-1, + -0x1.d906bcp-1, -0x1.e9f416p-1, -0x1.f6297ep-1, -0x1.fd88dap-1, + -0x1p0, -0x1.fd88dap-1, -0x1.f6297cp-1, -0x1.e9f416p-1, + -0x1.d906bcp-1, -0x1.c38b3p-1, -0x1.a9b662p-1, -0x1.8bc806p-1, + -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 LIBC_INLINE int32_t range_reduction(float x, float &y) { + float kf = fputil::nearest_integer(x * 32); + y = fputil::multiply_add(x, 32.0, -kf); + + return static_cast(kf); +} + +LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { + using FPBits = typename fputil::FPBits; + FPBits xbits(x); + + uint16_t x_u = xbits.uintval(); + uint16_t x_abs = x_u & 0x7fff; + + // Range reduction: + // For |x| > 1/32, we perform range reduction as follows: + // Find k and y such that: + // x = (k + y) * 1/32 + // k is an integer + // |y| < 0.5 + // + // This is done by performing: + // k = round(x * 32) + // y = x * 32 - k + // + // Once k and y are computed, we then deduce the answer by the sine of sum + // formula: + // sin(x * pi) = sin((k + y) * pi/32) + // = sin(k * pi/32) * cos(y * pi/32) + sin (y * pi/32) * cos (k * + // pi/32) + // The values of sin(k * pi/32) and cos (k * pi/32) for k = 0...63 are + // precomputed and stored using a vector of 64 single precision floats. sin(y + // * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev + // polynomials generated by Sollya. + + // For signed zeros + if (LIBC_UNLIKELY(x_abs == 0U)) return fputil::cast(1.0f); + + // Numbers greater or equal to 2^10 are integers, or infinity, or NaN + if (LIBC_UNLIKELY(x_abs >= 0x6400)) { + if (LIBC_UNLIKELY(x_abs <= 0x67FF)) { + return fputil::cast((x_abs & 0x1) ? -1.0f : 1.0f); + } + + // Check for NaN or infintiy values + 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(); + } + + return fputil::cast(1.0f); + } + + + float f32 = x; + float y; + int32_t k = range_reduction(f32, y); + + float sin_k = SIN_K_PI_OVER_32[k & 63]; + float cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; + + // Recall; + // cos(x * pi/32) = cos((k + y) * pi/32) + // = cos(y * pi/32) * cos(k * pi/32) + // - sin(y * pi/32) * sin(k * pi/32) + // Recall, after range reduction, -0.5 <= y <= 0.5. For very small + // values of y, calculating sin(y * p/32) can be inaccurate. Generating a + // polynomial for sin(y * p/32)/y instead significantly reduces the relative + // errors. + float ysq = y * y; + + // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya + // with: + // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]); + float sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f, + 0x1.a03354p-21f, -0x1.ad02d2p-20f); + + // Note that cosm1_y = cos(y*pi/32) - 1 = cos_y - 1 + // Derivation: + // cos(x * pi) = cos((k + y) * pi/32) + // = cos_k * cos_y + sin_k * sin_y + // = cos_k * (1 + cos_y - 1) + sin_k * sin_y + // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with: + // > P = fpminimax(cos(y*pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]); + float cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, + 0x1.a6f7a2p-29f); + + if (LIBC_UNLIKELY(sin_y == 0 && cos_k == 0)) + return fputil::cast(0.0f); + + // Since, cosm1_y = cos_y - 1, therefore: + // cos(x * pi) = cos_k(cosm1_y) + cos_k - sin_k * sin_y + return fputil::cast(fputil::multiply_add(cos_k, cosm1_y, fputil::multiply_add(-sin_k, sin_y, cos_k))); +} +} // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index 24a5abec898a8..c0209a1328702 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -45,6 +45,17 @@ add_fp_unittest( ) +add_fp_unittest( + cospif16_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + cospif16_test.cpp + DEPENDS + libc.src.math.cospif16 +) + add_fp_unittest( daddl_test NEED_MPFR diff --git a/libc/test/src/math/cospif16_test.cpp b/libc/test/src/math/cospif16_test.cpp new file mode 100644 index 0000000000000..1ba70c31cc151 --- /dev/null +++ b/libc/test/src/math/cospif16_test.cpp @@ -0,0 +1,38 @@ +//===-- Exhaustive test for cospif16 --------------------------------------===// +// +// 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/cospif16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +using LlvmLibcCospif16Test = 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(LlvmLibcCospif16Test, PositiveRange) { + for (uint16_t v = POS_START; v <= POS_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, x, LIBC_NAMESPACE::cospif16(x), 0.5); + } +} + +TEST_F(LlvmLibcCospif16Test, NegativeRange) { + for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, x, LIBC_NAMESPACE::cospif16(x), 0.5); + } +} diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index 3c077240356bd..8bba7a083da4d 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -25,6 +25,17 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + cospif16_test + SUITE + libc-math-smoke-tests + SRCS + cospif16_test.cpp + DEPENDS + libc.src.errno.errno + libc.src.math.cospif16 +) + add_fp_unittest( sinf_test SUITE diff --git a/libc/test/src/math/smoke/cospif16_test.cpp b/libc/test/src/math/smoke/cospif16_test.cpp new file mode 100644 index 0000000000000..7daa5d2956c10 --- /dev/null +++ b/libc/test/src/math/smoke/cospif16_test.cpp @@ -0,0 +1,44 @@ +//===-- Unittests for sinpif16 --------------------------------------------===// +// +// 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/cospif16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" + +using LlvmLibcCospif16Test = LIBC_NAMESPACE::testing::FPTest; + +TEST_F(LlvmLibcCospif16Test, SpecialNumbers) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::cospif16(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::cospif16(zero)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::cospif16(neg_zero)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::cospif16(inf)); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::cospif16(neg_inf)); + EXPECT_MATH_ERRNO(EDOM); +} + +TEST_F(LlvmLibcCospif16Test, Integers) { + EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::cospif16(-0x420)); + EXPECT_FP_EQ(1.0f, LIBC_NAMESPACE::cospif16(-0x1.4p+14)); + EXPECT_FP_EQ(-1.0f, LIBC_NAMESPACE::cospif16(0x421)); + EXPECT_FP_EQ(-1.0f, LIBC_NAMESPACE::cospif16(0x333)); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::cospif16(-0x1.28p4)); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::cospif16(-0x1.ffcp9)); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::cospif16(0x1.01p7)); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::cospif16(0x1.f6cp9)); +} From b6e0457809c4ce8033300f706802e3c79cbe8237 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 00:20:13 +0100 Subject: [PATCH 02/15] clang format --- libc/docs/math/index.rst | 2 +- libc/src/math/generic/cospif16.cpp | 23 ++++++++++++----------- libc/test/src/math/cospif16_test.cpp | 6 ++++-- 3 files changed, 17 insertions(+), 14 deletions(-) diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst index 6591cbbdc1558..84ed0e4135eba 100644 --- a/libc/docs/math/index.rst +++ b/libc/docs/math/index.rst @@ -280,7 +280,7 @@ Higher Math Functions +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | cosh | |check| | | | |check| | | 7.12.5.4 | F.10.2.4 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| cospi | |check| | | | | | 7.12.4.12 | F.10.1.12 | +| cospi | |check| | | | |check| | | 7.12.4.12 | F.10.1.12 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | dsqrt | N/A | N/A | |check| | N/A | |check|\* | 7.12.14.6 | F.10.11 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/src/math/generic/cospif16.cpp b/libc/src/math/generic/cospif16.cpp index 779aab781b9f0..90deb4a77dad9 100644 --- a/libc/src/math/generic/cospif16.cpp +++ b/libc/src/math/generic/cospif16.cpp @@ -75,28 +75,28 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { // polynomials generated by Sollya. // For signed zeros - if (LIBC_UNLIKELY(x_abs == 0U)) return fputil::cast(1.0f); + if (LIBC_UNLIKELY(x_abs == 0U)) + return fputil::cast(1.0f); // Numbers greater or equal to 2^10 are integers, or infinity, or NaN if (LIBC_UNLIKELY(x_abs >= 0x6400)) { if (LIBC_UNLIKELY(x_abs <= 0x67FF)) { - return fputil::cast((x_abs & 0x1) ? -1.0f : 1.0f); + return fputil::cast((x_abs & 0x1) ? -1.0f : 1.0f); } - + // Check for NaN or infintiy values 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); + fputil::raise_except_if_required(FE_INVALID); } return x + FPBits::quiet_nan().get_val(); } - - return fputil::cast(1.0f); + + return fputil::cast(1.0f); } - float f32 = x; float y; @@ -107,8 +107,8 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { // Recall; // cos(x * pi/32) = cos((k + y) * pi/32) - // = cos(y * pi/32) * cos(k * pi/32) - // - sin(y * pi/32) * sin(k * pi/32) + // = cos(y * pi/32) * cos(k * pi/32) + // - sin(y * pi/32) * sin(k * pi/32) // Recall, after range reduction, -0.5 <= y <= 0.5. For very small // values of y, calculating sin(y * p/32) can be inaccurate. Generating a // polynomial for sin(y * p/32)/y instead significantly reduces the relative @@ -116,7 +116,7 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { float ysq = y * y; // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya - // with: + // with: // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]); float sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f, 0x1.a03354p-21f, -0x1.ad02d2p-20f); @@ -136,6 +136,7 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { // Since, cosm1_y = cos_y - 1, therefore: // cos(x * pi) = cos_k(cosm1_y) + cos_k - sin_k * sin_y - return fputil::cast(fputil::multiply_add(cos_k, cosm1_y, fputil::multiply_add(-sin_k, sin_y, cos_k))); + return fputil::cast(fputil::multiply_add( + cos_k, cosm1_y, fputil::multiply_add(-sin_k, sin_y, cos_k))); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/src/math/cospif16_test.cpp b/libc/test/src/math/cospif16_test.cpp index 1ba70c31cc151..2c61086ef8735 100644 --- a/libc/test/src/math/cospif16_test.cpp +++ b/libc/test/src/math/cospif16_test.cpp @@ -26,13 +26,15 @@ static constexpr uint16_t NEG_STOP = 0xfc00U; TEST_F(LlvmLibcCospif16Test, PositiveRange) { for (uint16_t v = POS_START; v <= POS_STOP; ++v) { float16 x = FPBits(v).get_val(); - EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, x, LIBC_NAMESPACE::cospif16(x), 0.5); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, x, + LIBC_NAMESPACE::cospif16(x), 0.5); } } TEST_F(LlvmLibcCospif16Test, NegativeRange) { for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) { float16 x = FPBits(v).get_val(); - EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, x, LIBC_NAMESPACE::cospif16(x), 0.5); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Cospi, x, + LIBC_NAMESPACE::cospif16(x), 0.5); } } From 9a1e08becb07fc6a352c165bebc29825318e60fc Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 00:51:49 +0100 Subject: [PATCH 03/15] doc update --- libc/src/math/generic/cospif16.cpp | 6 +++--- libc/test/src/math/smoke/cospif16_test.cpp | 2 +- 2 files changed, 4 insertions(+), 4 deletions(-) diff --git a/libc/src/math/generic/cospif16.cpp b/libc/src/math/generic/cospif16.cpp index 90deb4a77dad9..662f1101862b9 100644 --- a/libc/src/math/generic/cospif16.cpp +++ b/libc/src/math/generic/cospif16.cpp @@ -66,9 +66,9 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { // // Once k and y are computed, we then deduce the answer by the sine of sum // formula: - // sin(x * pi) = sin((k + y) * pi/32) - // = sin(k * pi/32) * cos(y * pi/32) + sin (y * pi/32) * cos (k * - // pi/32) + // cos(x * pi) = cos((k + y) * pi/32) + // = cos(k * pi/32) * cos(y * pi/32) + // + sin(y * pi/32) * sin(k * pi/32) // The values of sin(k * pi/32) and cos (k * pi/32) for k = 0...63 are // precomputed and stored using a vector of 64 single precision floats. sin(y // * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev diff --git a/libc/test/src/math/smoke/cospif16_test.cpp b/libc/test/src/math/smoke/cospif16_test.cpp index 7daa5d2956c10..f6d7483393191 100644 --- a/libc/test/src/math/smoke/cospif16_test.cpp +++ b/libc/test/src/math/smoke/cospif16_test.cpp @@ -1,4 +1,4 @@ -//===-- Unittests for sinpif16 --------------------------------------------===// +//===-- Unittests for cospif16 --------------------------------------------===// // // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. // See https://llvm.org/LICENSE.txt for license information. From bf93bb460dde909e5a0b13dc2b72ed17c4cb0989 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 00:53:45 +0100 Subject: [PATCH 04/15] clang-format --- libc/src/math/generic/cospif16.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/src/math/generic/cospif16.cpp b/libc/src/math/generic/cospif16.cpp index 662f1101862b9..21c5722d244fa 100644 --- a/libc/src/math/generic/cospif16.cpp +++ b/libc/src/math/generic/cospif16.cpp @@ -67,7 +67,7 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { // Once k and y are computed, we then deduce the answer by the sine of sum // formula: // cos(x * pi) = cos((k + y) * pi/32) - // = cos(k * pi/32) * cos(y * pi/32) + // = cos(k * pi/32) * cos(y * pi/32) // + sin(y * pi/32) * sin(k * pi/32) // The values of sin(k * pi/32) and cos (k * pi/32) for k = 0...63 are // precomputed and stored using a vector of 64 single precision floats. sin(y From 697397e749fd9e05981cfeda4deb5aedc766cf91 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 08:06:43 +0100 Subject: [PATCH 05/15] added sincosf16_utils header, fixed formatting in existing pr --- libc/src/math/cospif16.h | 4 +- libc/src/math/generic/CMakeLists.txt | 27 ++++---- libc/src/math/generic/cospif16.cpp | 69 ++------------------ libc/src/math/generic/sincosf16_utils.h | 83 +++++++++++++++++++++++++ libc/src/math/generic/sinpif16.cpp | 72 ++------------------- libc/test/src/math/cospif16_test.cpp | 2 +- 6 files changed, 111 insertions(+), 146 deletions(-) create mode 100644 libc/src/math/generic/sincosf16_utils.h diff --git a/libc/src/math/cospif16.h b/libc/src/math/cospif16.h index 6779e67cdccae..ef9625dfed45f 100644 --- a/libc/src/math/cospif16.h +++ b/libc/src/math/cospif16.h @@ -1,10 +1,10 @@ -//===-- Implementation header for cospif16 ---------------------*- C++ -*-===// +//===-- Implementation header for cospif16 ----------------------*- 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_COSPIF16_H #define LLVM_LIBC_SRC_MATH_COSPIF16_H diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index e535c9e50b9b8..f45c88f1af755 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -351,6 +351,18 @@ add_header_library( libc.src.__support.common ) +add_header_library( + sincosf16_utils + HDRS + sincosf16_utils.h + DEPENDS + .range_reduction + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.polyeval + libc.src.__support.FPUtil.nearest_integer + libc.src.__support.common +) + add_header_library( sincos_eval HDRS @@ -422,7 +434,6 @@ add_entrypoint_object( -O3 ) - add_entrypoint_object( cospif16 SRCS @@ -430,14 +441,10 @@ add_entrypoint_object( HDRS ../cospif16.h DEPENDS - libc.src.__support.common + .sincosf16_utils 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.FPUtil.nearest_integer - libc.src.__support.FPUtil.polyeval - libc.src.__support.macros.properties.types + libc.src.__support.FPUtil.fp_bits COMPILE_OPTIONS -O3 ) @@ -555,14 +562,10 @@ add_entrypoint_object( HDRS ../sinpif16.h DEPENDS - libc.src.__support.common + .sincosf16_utils 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.FPUtil.nearest_integer - libc.src.__support.FPUtil.polyeval - libc.src.__support.macros.properties.types COMPILE_OPTIONS -O3 ) diff --git a/libc/src/math/generic/cospif16.cpp b/libc/src/math/generic/cospif16.cpp index 21c5722d244fa..90d781bb830fe 100644 --- a/libc/src/math/generic/cospif16.cpp +++ b/libc/src/math/generic/cospif16.cpp @@ -7,51 +7,20 @@ //===----------------------------------------------------------------------===// #include "src/math/cospif16.h" +#include "sincosf16_utils.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" #include "src/__support/FPUtil/cast.h" #include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/FPUtil/nearest_integer.h" -#include "src/__support/common.h" -#include "src/__support/macros/config.h" namespace LIBC_NAMESPACE_DECL { -// Lookup table for sin(k * pi / 32) with k = 0, ..., 63. -// Table is generated with Sollya as follows: -// > display = hexadecimal; -// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); }; -static constexpr float SIN_K_PI_OVER_32[64] = { - 0x0.0p0, 0x1.917a6cp-4, 0x1.8f8b84p-3, 0x1.294062p-2, - 0x1.87de2ap-2, 0x1.e2b5d4p-2, 0x1.1c73b4p-1, 0x1.44cf32p-1, - 0x1.6a09e6p-1, 0x1.8bc806p-1, 0x1.a9b662p-1, 0x1.c38b3p-1, - 0x1.d906bcp-1, 0x1.e9f416p-1, 0x1.f6297cp-1, 0x1.fd88dap-1, - 0x1p0, 0x1.fd88dap-1, 0x1.f6297cp-1, 0x1.e9f416p-1, - 0x1.d906bcp-1, 0x1.c38b3p-1, 0x1.a9b662p-1, 0x1.8bc806p-1, - 0x1.6a09e6p-1, 0x1.44cf32p-1, 0x1.1c73b4p-1, 0x1.e2b5d4p-2, - 0x1.87de2ap-2, 0x1.294062p-2, 0x1.8f8b84p-3, 0x1.917a6cp-4, - 0x0.0p0, -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294062p-2, - -0x1.87de2ap-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, -0x1.44cf32p-1, - -0x1.6a09e6p-1, -0x1.8bc806p-1, -0x1.a9b662p-1, -0x1.c38b3p-1, - -0x1.d906bcp-1, -0x1.e9f416p-1, -0x1.f6297ep-1, -0x1.fd88dap-1, - -0x1p0, -0x1.fd88dap-1, -0x1.f6297cp-1, -0x1.e9f416p-1, - -0x1.d906bcp-1, -0x1.c38b3p-1, -0x1.a9b662p-1, -0x1.8bc806p-1, - -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 LIBC_INLINE int32_t range_reduction(float x, float &y) { - float kf = fputil::nearest_integer(x * 32); - y = fputil::multiply_add(x, 32.0, -kf); - - return static_cast(kf); -} - LLVM_LIBC_FUNCTION(float16, cospif16, (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| > 1/32, we perform range reduction as follows: @@ -98,38 +67,8 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { return fputil::cast(1.0f); } - float f32 = x; - float y; - int32_t k = range_reduction(f32, y); - - float sin_k = SIN_K_PI_OVER_32[k & 63]; - float cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; - - // Recall; - // cos(x * pi/32) = cos((k + y) * pi/32) - // = cos(y * pi/32) * cos(k * pi/32) - // - sin(y * pi/32) * sin(k * pi/32) - // Recall, after range reduction, -0.5 <= y <= 0.5. For very small - // values of y, calculating sin(y * p/32) can be inaccurate. Generating a - // polynomial for sin(y * p/32)/y instead significantly reduces the relative - // errors. - float ysq = y * y; - - // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya - // with: - // > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]); - float sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f, - 0x1.a03354p-21f, -0x1.ad02d2p-20f); - - // Note that cosm1_y = cos(y*pi/32) - 1 = cos_y - 1 - // Derivation: - // cos(x * pi) = cos((k + y) * pi/32) - // = cos_k * cos_y + sin_k * sin_y - // = cos_k * (1 + cos_y - 1) + sin_k * sin_y - // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with: - // > P = fpminimax(cos(y*pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]); - float cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, - 0x1.a6f7a2p-29f); + float sin_k, cos_k, sin_y, cosm1_y; + sincospif16_eval(xf, sin_k, cos_k, sin_y, cosm1_y); if (LIBC_UNLIKELY(sin_y == 0 && cos_k == 0)) return fputil::cast(0.0f); diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h new file mode 100644 index 0000000000000..456e399e68f79 --- /dev/null +++ b/libc/src/math/generic/sincosf16_utils.h @@ -0,0 +1,83 @@ +//===-- Collection of utils for sinf16/cosf16 ------------------*- 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_GENERIC_SINCOSF16_UTILS_H +#define LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H + +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.h" +#include "src/__support/FPUtil/nearest_integer.h" +#include "src/__support/common.h" +#include "src/__support/macros/config.h" + + +namespace LIBC_NAMESPACE_DECL { + +// Lookup table for sin(k * pi / 32) with k = 0, ..., 63. +// Table is generated with Sollya as follows: +// > display = hexadecimmal; +// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); }; +constexpr float SIN_K_PI_OVER_32[64] = { + 0x0.0p0, 0x1.917a6cp-4, 0x1.8f8b84p-3, 0x1.294062p-2, + 0x1.87de2ap-2, 0x1.e2b5d4p-2, 0x1.1c73b4p-1, 0x1.44cf32p-1, + 0x1.6a09e6p-1, 0x1.8bc806p-1, 0x1.a9b662p-1, 0x1.c38b3p-1, + 0x1.d906bcp-1, 0x1.e9f416p-1, 0x1.f6297cp-1, 0x1.fd88dap-1, + 0x1p0, 0x1.fd88dap-1, 0x1.f6297cp-1, 0x1.e9f416p-1, + 0x1.d906bcp-1, 0x1.c38b3p-1, 0x1.a9b662p-1, 0x1.8bc806p-1, + 0x1.6a09e6p-1, 0x1.44cf32p-1, 0x1.1c73b4p-1, 0x1.e2b5d4p-2, + 0x1.87de2ap-2, 0x1.294062p-2, 0x1.8f8b84p-3, 0x1.917a6cp-4, + 0x0.0p0, -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294062p-2, + -0x1.87de2ap-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, -0x1.44cf32p-1, + -0x1.6a09e6p-1, -0x1.8bc806p-1, -0x1.a9b662p-1, -0x1.c38b3p-1, + -0x1.d906bcp-1, -0x1.e9f416p-1, -0x1.f6297ep-1, -0x1.fd88dap-1, + -0x1p0, -0x1.fd88dap-1, -0x1.f6297cp-1, -0x1.e9f416p-1, + -0x1.d906bcp-1, -0x1.c38b3p-1, -0x1.a9b662p-1, -0x1.8bc806p-1, + -0x1.6a09e6p-1, -0x1.44cf32p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2, + -0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4}; + +LIBC_INLINE int32_t range_reduction(float x, float &y) { + float kf = fputil::nearest_integer(x * 32); + y = fputil::multiply_add(x, 32.0, -kf); + + 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(xf, y); + + sin_k = SIN_K_PI_OVER_32[k & 63]; + cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; + + // Recall; + // sin(x * pi/32) = sin((k + y) * pi/32) + // = sin(y * pi/32) * cos(k * pi/32) + cos(y * pi/32) * sin(k * + // pi/32) Recall, after range reduction, -0.5 <= y <= 0.5. For very small + // values of y, calculating sin(y * p/32) can be inaccurate. Generating a + // polynomial for sin(y * p/32)/y instead significantly reduces the relative + // errors. + float ysq = y * y; + + // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya + // with: > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]); + sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f, + 0x1.a03354p-21f, -0x1.ad02d2p-20f); + + // Note that cosm1_y = cos(y*pi/32) - 1 = cos_y - 1 + // Derivation: + // sin(x * pi) = sin((k + y) * pi/32) + // = sin_y * cos_k + cos_y * sin_k + // = cos_k * sin_y + sin_k * (1 + cos_y - 1) + // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with: + // > P = fpminimax(cos(y*pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]); + cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, + 0x1.a6f7a2p-29f); +} + +} // namespace LIBC_NAMESPACE_DECL + +#endif // LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 17cca583e0c0e..2c52d4cd42688 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -7,52 +7,20 @@ //===----------------------------------------------------------------------===// #include "src/math/sinpif16.h" +#include "sincosf16_utils.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" -#include "src/__support/FPUtil/PolyEval.h" #include "src/__support/FPUtil/cast.h" #include "src/__support/FPUtil/multiply_add.h" -#include "src/__support/FPUtil/nearest_integer.h" -#include "src/__support/common.h" -#include "src/__support/macros/config.h" namespace LIBC_NAMESPACE_DECL { - -// Lookup table for sin(k * pi / 32) with k = 0, ..., 63. -// Table is generated with Sollya as follows: -// > display = hexadecimmal; -// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); }; -static constexpr float SIN_K_PI_OVER_32[64] = { - 0x0.0p0, 0x1.917a6cp-4, 0x1.8f8b84p-3, 0x1.294062p-2, - 0x1.87de2ap-2, 0x1.e2b5d4p-2, 0x1.1c73b4p-1, 0x1.44cf32p-1, - 0x1.6a09e6p-1, 0x1.8bc806p-1, 0x1.a9b662p-1, 0x1.c38b3p-1, - 0x1.d906bcp-1, 0x1.e9f416p-1, 0x1.f6297cp-1, 0x1.fd88dap-1, - 0x1p0, 0x1.fd88dap-1, 0x1.f6297cp-1, 0x1.e9f416p-1, - 0x1.d906bcp-1, 0x1.c38b3p-1, 0x1.a9b662p-1, 0x1.8bc806p-1, - 0x1.6a09e6p-1, 0x1.44cf32p-1, 0x1.1c73b4p-1, 0x1.e2b5d4p-2, - 0x1.87de2ap-2, 0x1.294062p-2, 0x1.8f8b84p-3, 0x1.917a6cp-4, - 0x0.0p0, -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294062p-2, - -0x1.87de2ap-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, -0x1.44cf32p-1, - -0x1.6a09e6p-1, -0x1.8bc806p-1, -0x1.a9b662p-1, -0x1.c38b3p-1, - -0x1.d906bcp-1, -0x1.e9f416p-1, -0x1.f6297ep-1, -0x1.fd88dap-1, - -0x1p0, -0x1.fd88dap-1, -0x1.f6297cp-1, -0x1.e9f416p-1, - -0x1.d906bcp-1, -0x1.c38b3p-1, -0x1.a9b662p-1, -0x1.8bc806p-1, - -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 LIBC_INLINE int32_t range_reduction(float x, float &y) { - float kf = fputil::nearest_integer(x * 32); - y = fputil::multiply_add(x, 32.0, -kf); - - return static_cast(kf); -} - LLVM_LIBC_FUNCTION(float16, sinpif16, (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| > 1/32, we perform range reduction as follows: @@ -68,8 +36,8 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // Once k and y are computed, we then deduce the answer by the sine of sum // formula: // sin(x * pi) = sin((k + y) * pi/32) - // = sin(k * pi/32) * cos(y * pi/32) + sin (y * pi/32) * cos (k * - // pi/32) + // = sin(k * pi/32) * cos(y * pi/32) + // + sin (y * pi/32) * cos (k * pi/32) // The values of sin(k * pi/32) and cos (k * pi/32) for k = 0...63 are // precomputed and stored using a vector of 64 single precision floats. sin(y // * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev @@ -94,36 +62,8 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { return FPBits::zero(xbits.sign()).get_val(); } - float f32 = x; - float y; - int32_t k = range_reduction(f32, y); - - float sin_k = SIN_K_PI_OVER_32[k & 63]; - float cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; - - // Recall; - // sin(x * pi/32) = sin((k + y) * pi/32) - // = sin(y * pi/32) * cos(k * pi/32) + cos(y * pi/32) * sin(k * - // pi/32) Recall, after range reduction, -0.5 <= y <= 0.5. For very small - // values of y, calculating sin(y * p/32) can be inaccurate. Generating a - // polynomial for sin(y * p/32)/y instead significantly reduces the relative - // errors. - float ysq = y * y; - - // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya - // with: > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]); - float sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f, - 0x1.a03354p-21f, -0x1.ad02d2p-20f); - - // Note that cosm1_y = cos(y*pi/32) - 1 = cos_y - 1 - // Derivation: - // sin(x * pi) = sin((k + y) * pi/32) - // = sin_y * cos_k + cos_y * sin_k - // = cos_k * sin_y + sin_k * (1 + cos_y - 1) - // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with: - // > P = fpminimax(cos(y*pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]); - float cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, - 0x1.a6f7a2p-29f); + float sin_k, cos_k, sin_y, cosm1_y; + sincospif16_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(); diff --git a/libc/test/src/math/cospif16_test.cpp b/libc/test/src/math/cospif16_test.cpp index 2c61086ef8735..6a32498b0570a 100644 --- a/libc/test/src/math/cospif16_test.cpp +++ b/libc/test/src/math/cospif16_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/cospif16.h" #include "test/UnitTest/FPMatcher.h" From a7f9f81f11f3c271b7151d2c0955bdb7744222f6 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 08:11:55 +0100 Subject: [PATCH 06/15] formatting --- libc/src/math/generic/sincosf16_utils.h | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h index 456e399e68f79..89f0d931fcf20 100644 --- a/libc/src/math/generic/sincosf16_utils.h +++ b/libc/src/math/generic/sincosf16_utils.h @@ -1,10 +1,10 @@ -//===-- Collection of utils for sinf16/cosf16 ------------------*- C++ -*-===// +//===-- Collection of utils for sinf16/cosf16 -------------------*- 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_GENERIC_SINCOSF16_UTILS_H #define LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H @@ -14,7 +14,6 @@ #include "src/__support/common.h" #include "src/__support/macros/config.h" - namespace LIBC_NAMESPACE_DECL { // Lookup table for sin(k * pi / 32) with k = 0, ..., 63. @@ -46,10 +45,11 @@ LIBC_INLINE int32_t range_reduction(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) { +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(xf, y); - + sin_k = SIN_K_PI_OVER_32[k & 63]; cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; @@ -65,7 +65,7 @@ LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, float &s // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya // with: > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]); sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f, - 0x1.a03354p-21f, -0x1.ad02d2p-20f); + 0x1.a03354p-21f, -0x1.ad02d2p-20f); // Note that cosm1_y = cos(y*pi/32) - 1 = cos_y - 1 // Derivation: @@ -75,7 +75,7 @@ LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, float &s // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with: // > P = fpminimax(cos(y*pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]); cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, - 0x1.a6f7a2p-29f); + 0x1.a6f7a2p-29f); } } // namespace LIBC_NAMESPACE_DECL From d63cf622540135ebcd30b970d4bf34c8c3241cf5 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 09:13:58 +0100 Subject: [PATCH 07/15] fixed mpfr_cospi() evaluation for MPFR version < 4.2 --- libc/utils/MPFRWrapper/MPFRUtils.cpp | 18 ++++-------------- 1 file changed, 4 insertions(+), 14 deletions(-) diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index bd4fbe294a622..9295280f6b7d0 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -255,19 +255,9 @@ class MPFRNumber { mpfr_cospi(result.value, value, mpfr_rounding); return result; #else - MPFRNumber value_frac(*this); - mpfr_frac(value_frac.value, value, MPFR_RNDN); - - if (mpfr_cmp_si(value_frac.value, 0.0) == 0) { - mpz_t integer_part; - mpz_init(integer_part); - mpfr_get_z(integer_part, value, MPFR_RNDN); - - if (mpz_tstbit(integer_part, 0)) { - mpfr_set_si(result.value, -1.0, MPFR_RNDN); // odd - } else { - mpfr_set_si(result.value, 1.0, MPFR_RNDN); // even - } + if (mpfr_integer_p(value)) { + auto d = mpfr_get_si(value, mpfr_rounding); + mpfr_set_si(result.value, (d & 1) ? -1 : 1, mpfr_rounding); return result; } @@ -277,7 +267,7 @@ class MPFRNumber { mpfr_cos(result.value, value_pi.value, mpfr_rounding); return result; -#endif +//#endif } MPFRNumber erf() const { From 5d620c2ecbf4e29f5198e21abbb21165d5cb9c2a Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 09:41:01 +0100 Subject: [PATCH 08/15] minor fix --- libc/utils/MPFRWrapper/MPFRUtils.cpp | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 9295280f6b7d0..92cac8c7ffe3c 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -267,7 +267,7 @@ class MPFRNumber { mpfr_cos(result.value, value_pi.value, mpfr_rounding); return result; -//#endif +#endif } MPFRNumber erf() const { From 6b428aab2ac6046074787b712d040fb80cf6a060 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 10:25:05 +0100 Subject: [PATCH 09/15] improved fix for mpfr_cospi() implementation for when MPFR version < 4.2 --- libc/utils/MPFRWrapper/MPFRUtils.cpp | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 92cac8c7ffe3c..3538416698f8a 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -256,8 +256,11 @@ class MPFRNumber { return result; #else if (mpfr_integer_p(value)) { - auto d = mpfr_get_si(value, mpfr_rounding); - mpfr_set_si(result.value, (d & 1) ? -1 : 1, mpfr_rounding); + mpz_t integer; + mpz_init(integer); + mpfr_get_z(integer, value, mpfr_rounding); + + mpfr_set_si(result.value, (mpz_tstbit(integer, 0)) ? -1 : 1, mpfr_rounding); return result; } From b606662b90951a2fc8e8db1d796ac4cd560e54dc Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sat, 19 Oct 2024 10:54:11 +0100 Subject: [PATCH 10/15] clang format...again --- libc/utils/MPFRWrapper/MPFRUtils.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 3538416698f8a..d9cfd00670746 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -260,7 +260,8 @@ class MPFRNumber { mpz_init(integer); mpfr_get_z(integer, value, mpfr_rounding); - mpfr_set_si(result.value, (mpz_tstbit(integer, 0)) ? -1 : 1, mpfr_rounding); + auto d = mpz_tstbit(integer, 0); + mpfr_set_si(result.value, d ? -1 : 1, mpfr_rounding); return result; } From c45bf585f062fd5ad75c233e0a7ac165074205e3 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Sun, 20 Oct 2024 14:36:31 +0100 Subject: [PATCH 11/15] nit --- libc/src/math/generic/CMakeLists.txt | 1 - libc/src/math/generic/cospif16.cpp | 7 +++++-- libc/src/math/generic/sincosf16_utils.h | 5 +++-- 3 files changed, 8 insertions(+), 5 deletions(-) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index f45c88f1af755..145c070879660 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -356,7 +356,6 @@ add_header_library( HDRS sincosf16_utils.h DEPENDS - .range_reduction libc.src.__support.FPUtil.fp_bits libc.src.__support.FPUtil.polyeval libc.src.__support.FPUtil.nearest_integer diff --git a/libc/src/math/generic/cospif16.cpp b/libc/src/math/generic/cospif16.cpp index 90d781bb830fe..eda1c20f8c17c 100644 --- a/libc/src/math/generic/cospif16.cpp +++ b/libc/src/math/generic/cospif16.cpp @@ -7,13 +7,17 @@ //===----------------------------------------------------------------------===// #include "src/math/cospif16.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, cospif16, (float16 x)) { using FPBits = typename fputil::FPBits; FPBits xbits(x); @@ -49,9 +53,8 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { // Numbers greater or equal to 2^10 are integers, or infinity, or NaN if (LIBC_UNLIKELY(x_abs >= 0x6400)) { - if (LIBC_UNLIKELY(x_abs <= 0x67FF)) { + if (LIBC_UNLIKELY(x_abs <= 0x67FF)) return fputil::cast((x_abs & 0x1) ? -1.0f : 1.0f); - } // Check for NaN or infintiy values if (LIBC_UNLIKELY(x_abs >= 0x7c00)) { diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h index 89f0d931fcf20..39121f58f2024 100644 --- a/libc/src/math/generic/sincosf16_utils.h +++ b/libc/src/math/generic/sincosf16_utils.h @@ -5,6 +5,7 @@ // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //===----------------------------------------------------------------------===// + #ifndef LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H #define LLVM_LIBC_SRC_MATH_GENERIC_SINCOSF16_UTILS_H @@ -38,7 +39,7 @@ 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}; -LIBC_INLINE int32_t range_reduction(float x, float &y) { +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); @@ -48,7 +49,7 @@ LIBC_INLINE int32_t range_reduction(float x, float &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(xf, y); + int32_t k = range_reduction_sincospif16(xf, y); sin_k = SIN_K_PI_OVER_32[k & 63]; cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; From a3d531597fe3a90dc011a853c641a46fc2203974 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Thu, 24 Oct 2024 01:45:26 +0100 Subject: [PATCH 12/15] Update libc/src/math/generic/CMakeLists.txt Co-authored-by: OverMighty --- libc/src/math/generic/CMakeLists.txt | 4 ++++ 1 file changed, 4 insertions(+) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 9491bd178f13d..2e1f4ca31fc8a 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -441,9 +441,13 @@ add_entrypoint_object( ../cospif16.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 ) From 5db2256ae7dbdf98acfc0d94021ffa89859452d4 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Tue, 29 Oct 2024 10:32:28 +0100 Subject: [PATCH 13/15] nit formatting and header updates --- libc/src/math/generic/CMakeLists.txt | 6 +++++- libc/src/math/generic/cospif16.cpp | 9 ++------- libc/src/math/generic/sincosf16_utils.h | 12 ++++++------ libc/src/math/generic/sinpif16.cpp | 15 +++++++-------- libc/utils/MPFRWrapper/MPFRUtils.cpp | 2 +- 5 files changed, 21 insertions(+), 23 deletions(-) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 2e1f4ca31fc8a..ca27759d3212f 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -566,9 +566,13 @@ add_entrypoint_object( ../sinpif16.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.fp_bits + libc.src.__support.FPUtil.multiply_add + libc.src.__support.macros.optimization COMPILE_OPTIONS -O3 ) diff --git a/libc/src/math/generic/cospif16.cpp b/libc/src/math/generic/cospif16.cpp index eda1c20f8c17c..af9bf73291ddb 100644 --- a/libc/src/math/generic/cospif16.cpp +++ b/libc/src/math/generic/cospif16.cpp @@ -17,7 +17,6 @@ #include "src/__support/macros/optimization.h" namespace LIBC_NAMESPACE_DECL { - LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { using FPBits = typename fputil::FPBits; FPBits xbits(x); @@ -40,12 +39,8 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { // Once k and y are computed, we then deduce the answer by the sine of sum // formula: // cos(x * pi) = cos((k + y) * pi/32) - // = cos(k * pi/32) * cos(y * pi/32) - // + sin(y * pi/32) * sin(k * pi/32) - // The values of sin(k * pi/32) and cos (k * pi/32) for k = 0...63 are - // precomputed and stored using a vector of 64 single precision floats. sin(y - // * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev - // polynomials generated by Sollya. + // = cos(k * pi/32) * cos(y * pi/32) + + // sin(y * pi/32) * sin(k * pi/32) // For signed zeros if (LIBC_UNLIKELY(x_abs == 0U)) diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h index 39121f58f2024..6cad95285f3e6 100644 --- a/libc/src/math/generic/sincosf16_utils.h +++ b/libc/src/math/generic/sincosf16_utils.h @@ -56,11 +56,11 @@ LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, // Recall; // sin(x * pi/32) = sin((k + y) * pi/32) - // = sin(y * pi/32) * cos(k * pi/32) + cos(y * pi/32) * sin(k * - // pi/32) Recall, after range reduction, -0.5 <= y <= 0.5. For very small - // values of y, calculating sin(y * p/32) can be inaccurate. Generating a - // polynomial for sin(y * p/32)/y instead significantly reduces the relative - // errors. + // = sin(y * pi/32) * cos(k * pi/32) + + // cos(y * pi/32) * sin(k * pi/32) + // Recall, after range reduction, -0.5 <= y <= 0.5. For very small values of + // y, calculating sin(y * p/32) can be inaccurate. Generating a polynomial for + // sin(y * p/32)/y instead significantly reduces the relative errors. float ysq = y * y; // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya @@ -74,7 +74,7 @@ LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, // = sin_y * cos_k + cos_y * sin_k // = cos_k * sin_y + sin_k * (1 + cos_y - 1) // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with: - // > P = fpminimax(cos(y*pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]); + // > P = fpminimax(cos(y * pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]); cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, 0x1.a6f7a2p-29f); } diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 2c52d4cd42688..e663cc0774b32 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -6,12 +6,15 @@ // //===----------------------------------------------------------------------===// + #include "src/math/sinpif16.h" -#include "sincosf16_utils.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/cast.h" #include "src/__support/FPUtil/multiply_add.h" +#include "sincosf16_utils.h" namespace LIBC_NAMESPACE_DECL { LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { @@ -36,13 +39,9 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // Once k and y are computed, we then deduce the answer by the sine of sum // formula: // sin(x * pi) = sin((k + y) * pi/32) - // = sin(k * pi/32) * cos(y * pi/32) - // + sin (y * pi/32) * cos (k * pi/32) - // The values of sin(k * pi/32) and cos (k * pi/32) for k = 0...63 are - // precomputed and stored using a vector of 64 single precision floats. sin(y - // * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev - // polynomials generated by Sollya. - + // = sin(k * pi/32) * cos(y * pi/32) + + // sin(y * pi/32) * cos(k * pi/32) + // For signed zeros if (LIBC_UNLIKELY(x_abs == 0U)) return x; diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index d9cfd00670746..60e4abadb5e3c 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -260,7 +260,7 @@ class MPFRNumber { mpz_init(integer); mpfr_get_z(integer, value, mpfr_rounding); - auto d = mpz_tstbit(integer, 0); + int d = mpz_tstbit(integer, 0); mpfr_set_si(result.value, d ? -1 : 1, mpfr_rounding); return result; } From 31ea8dc35329055b61a3a29470ccc92c3b25fed5 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Tue, 29 Oct 2024 10:46:22 +0100 Subject: [PATCH 14/15] more formatting changes --- libc/src/math/generic/cospif16.cpp | 2 +- libc/src/math/generic/sincosf16_utils.h | 12 ++---------- libc/src/math/generic/sinpif16.cpp | 5 ++--- 3 files changed, 5 insertions(+), 14 deletions(-) diff --git a/libc/src/math/generic/cospif16.cpp b/libc/src/math/generic/cospif16.cpp index af9bf73291ddb..84ec7c0a6baa4 100644 --- a/libc/src/math/generic/cospif16.cpp +++ b/libc/src/math/generic/cospif16.cpp @@ -39,7 +39,7 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { // Once k and y are computed, we then deduce the answer by the sine of sum // formula: // cos(x * pi) = cos((k + y) * pi/32) - // = cos(k * pi/32) * cos(y * pi/32) + + // = cos(k * pi/32) * cos(y * pi/32) + // sin(y * pi/32) * sin(k * pi/32) // For signed zeros diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h index 6cad95285f3e6..135fa8f7cf73c 100644 --- a/libc/src/math/generic/sincosf16_utils.h +++ b/libc/src/math/generic/sincosf16_utils.h @@ -54,25 +54,17 @@ LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, sin_k = SIN_K_PI_OVER_32[k & 63]; cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; - // Recall; - // sin(x * pi/32) = sin((k + y) * pi/32) - // = sin(y * pi/32) * cos(k * pi/32) + - // cos(y * pi/32) * sin(k * pi/32) // Recall, after range reduction, -0.5 <= y <= 0.5. For very small values of // y, calculating sin(y * p/32) can be inaccurate. Generating a polynomial for // sin(y * p/32)/y instead significantly reduces the relative errors. float ysq = y * y; // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya - // with: > Q = fpminimax(sin(y*pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]); + // with: + // > Q = fpminimax(sin(y * pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]); sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f, 0x1.a03354p-21f, -0x1.ad02d2p-20f); - // Note that cosm1_y = cos(y*pi/32) - 1 = cos_y - 1 - // Derivation: - // sin(x * pi) = sin((k + y) * pi/32) - // = sin_y * cos_k + cos_y * sin_k - // = cos_k * sin_y + sin_k * (1 + cos_y - 1) // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with: // > P = fpminimax(cos(y * pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]); cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index e663cc0774b32..6642d7d578827 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -6,15 +6,14 @@ // //===----------------------------------------------------------------------===// - #include "src/math/sinpif16.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 "sincosf16_utils.h" namespace LIBC_NAMESPACE_DECL { LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { @@ -41,7 +40,7 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // sin(x * pi) = sin((k + y) * pi/32) // = sin(k * pi/32) * cos(y * pi/32) + // sin(y * pi/32) * cos(k * pi/32) - + // For signed zeros if (LIBC_UNLIKELY(x_abs == 0U)) return x; From e5dc117152ef5b22dd825635f7945df580ab17dd Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Tue, 29 Oct 2024 10:52:29 +0100 Subject: [PATCH 15/15] nit formatting changes --- libc/src/math/generic/cospif16.cpp | 2 ++ libc/src/math/generic/sincosf16_utils.h | 5 +++-- libc/src/math/generic/sinpif16.cpp | 2 ++ 3 files changed, 7 insertions(+), 2 deletions(-) diff --git a/libc/src/math/generic/cospif16.cpp b/libc/src/math/generic/cospif16.cpp index 84ec7c0a6baa4..dd8c7ab6afa3d 100644 --- a/libc/src/math/generic/cospif16.cpp +++ b/libc/src/math/generic/cospif16.cpp @@ -17,6 +17,7 @@ #include "src/__support/macros/optimization.h" namespace LIBC_NAMESPACE_DECL { + LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { using FPBits = typename fputil::FPBits; FPBits xbits(x); @@ -76,4 +77,5 @@ LLVM_LIBC_FUNCTION(float16, cospif16, (float16 x)) { return fputil::cast(fputil::multiply_add( cos_k, cosm1_y, fputil::multiply_add(-sin_k, sin_y, cos_k))); } + } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/generic/sincosf16_utils.h b/libc/src/math/generic/sincosf16_utils.h index 135fa8f7cf73c..83511755a56c4 100644 --- a/libc/src/math/generic/sincosf16_utils.h +++ b/libc/src/math/generic/sincosf16_utils.h @@ -60,12 +60,13 @@ LIBC_INLINE void sincospif16_eval(float xf, float &sin_k, float &cos_k, float ysq = y * y; // Degree-6 minimax even polynomial for sin(y*pi/32)/y generated by Sollya - // with: + // with: // > Q = fpminimax(sin(y * pi/32)/y, [|0, 2, 4, 6|], [|SG...|], [0, 0.5]); sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f, 0x1.a03354p-21f, -0x1.ad02d2p-20f); - // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with: + // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya + // with: // > P = fpminimax(cos(y * pi/32), [|0, 2, 4, 6|],[|1, SG...|], [0, 0.5]); cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, 0x1.a6f7a2p-29f); diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 6642d7d578827..51ea595653b4d 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -16,6 +16,7 @@ #include "src/__support/FPUtil/multiply_add.h" namespace LIBC_NAMESPACE_DECL { + LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { using FPBits = typename fputil::FPBits; FPBits xbits(x); @@ -71,4 +72,5 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { return fputil::cast(fputil::multiply_add( sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k))); } + } // namespace LIBC_NAMESPACE_DECL