From 1dd121c8f55d32c7152389e20d1cb2005d983630 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 23 Sep 2024 23:29:33 +0100 Subject: [PATCH 01/20] add sinpif16 function implementation and unittests --- libc/config/linux/aarch64/entrypoints.txt | 1 + libc/config/linux/arm/entrypoints.txt | 1 + libc/config/linux/riscv/entrypoints.txt | 1 + libc/config/linux/x86_64/entrypoints.txt | 1 + libc/config/windows/entrypoints.txt | 3 +- libc/newhdrgen/yaml/math.yaml | 6 + libc/src/math/CMakeLists.txt | 1 + libc/src/math/generic/CMakeLists.txt | 16 +++ libc/src/math/generic/sinpif16.cpp | 151 +++++++++++++++++++++ libc/src/math/sinpif16.h | 23 ++++ libc/test/src/math/smoke/CMakeLists.txt | 13 ++ libc/test/src/math/smoke/sinpif16_test.cpp | 44 ++++++ libc/utils/MPFRWrapper/MPFRUtils.h | 1 + 13 files changed, 261 insertions(+), 1 deletion(-) create mode 100644 libc/src/math/generic/sinpif16.cpp create mode 100644 libc/src/math/sinpif16.h create mode 100644 libc/test/src/math/smoke/sinpif16_test.cpp diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt index 64fbe1a250c0b..47f46c4b9160e 100644 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -564,6 +564,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sinf libc.src.math.sinhf libc.src.math.sinpif + libc.src.math.sinpif16 libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt index 1be9a872dd2f7..dbd21f6b5d438 100644 --- a/libc/config/linux/arm/entrypoints.txt +++ b/libc/config/linux/arm/entrypoints.txt @@ -390,6 +390,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sincosf libc.src.math.sinf libc.src.math.sinhf + libc.src.math.sinpif16 libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt index ff3d821c664c5..c45ba2fe5b9aa 100644 --- a/libc/config/linux/riscv/entrypoints.txt +++ b/libc/config/linux/riscv/entrypoints.txt @@ -567,6 +567,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sinf libc.src.math.sinhf libc.src.math.sinpif + libc.src.math.sinpif16 libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index dd658af3bfb67..12b86f541bc91 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -567,6 +567,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sinf libc.src.math.sinhf libc.src.math.sinpif + libc.src.math.sinpif16 libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt index 8f0b50bcc83ea..5b0d84669e924 100644 --- a/libc/config/windows/entrypoints.txt +++ b/libc/config/windows/entrypoints.txt @@ -275,7 +275,8 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sincosf libc.src.math.sincosf libc.src.math.sinf - libc.src.math.sinhf + libc.src.math.sinhfi + libc.src.math.sinpif16 libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/newhdrgen/yaml/math.yaml b/libc/newhdrgen/yaml/math.yaml index 04b6a073deace..fb79f15a31460 100644 --- a/libc/newhdrgen/yaml/math.yaml +++ b/libc/newhdrgen/yaml/math.yaml @@ -2290,6 +2290,12 @@ functions: return_type: float arguments: - type: float + - name: sinpif16 + standards: + - stdc + return_type: _Float16 + arguments: + - type: _Float16 - name: sqrt standards: - stdc diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt index 3cba34fc24932..dbfc0ce6fcaeb 100644 --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -462,6 +462,7 @@ add_math_entrypoint_object(sincosf) add_math_entrypoint_object(sin) add_math_entrypoint_object(sinf) add_math_entrypoint_object(sinpif) +add_math_entrypoint_object(sinpif16) add_math_entrypoint_object(sinh) add_math_entrypoint_object(sinhf) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 5a1ee3b8b83c7..6a2e45dbd4023 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -473,6 +473,22 @@ add_entrypoint_object( -O3 ) +add_entrypoint_object( + sinpif16 + SRCS + sinpif16.cpp + HDRS + ../sinpif16.h + DEPENDS + libc.src.__support.macros.properties.types + libc.src.__support.FPUtil.fp_bits + libc.src.__support.FPUtil.multiply_add + libc.src.__support.common + libc.src.__support.macros.config + COMPILE_OPTIONS + -O3 +) + add_entrypoint_object( tan SRCS diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp new file mode 100644 index 0000000000000..bb028d5b88fc0 --- /dev/null +++ b/libc/src/math/generic/sinpif16.cpp @@ -0,0 +1,151 @@ +//===-- Half-precision sinpif 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 +// +//===----------------------------------------------------------------------===// + +#define M_PI 3.1415925f + +#include "src/math/sinpif16.h" +#include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/multiply_add.h" +#include "src/__support/common.h" +#include "src/__support/macros/config.h" + +// TODO: Should probably create a new file; sincospif16_utils.h +// To store the following helper functions and constants. +// I'd defer to @lntue for suggestions regarding that + +// HELPER_START +namespace LIBC_NAMESPACE_DECL { + +constexpr float PI_OVER_32 = M_PI / 32; + +// In Sollya generate 10 coeffecients for a degree-9 chebyshev polynomial +// approximating the sine function in [-pi / 32, pi / 32] with the following +// commands: +// > prec=23; +// > TL = chebyshevform(sin(x), 9, [-pi / 32, pi / 32]); +// > TL[0]; +const float SIN_COEFF[10] = { + 0x1.801p-27, 0x1.000078p0, + -0x1.7e98p-14, -0x1.6bf4p-3, + 0x1.95ccp-5, 0x1.1baep2, + -0x1.030ap3, -0x1.3dap9, + 0x1.98e4p8, 0x1.d3d8p14 +}; + +// In Sollya generate 10 coefficients for a degree-9 chebyshev polynomial +// approximating the sine function in [-pi/32, pi/32] with the following +// commands: +// > prec = 23; +// > TL = chebyshevform(cos(x), 9, [-pi / 32, pi / 32]); +// > TL[0]; +const float COS_COEFF[10] = { + 0x1.00001p0, -0x1.48p-17, + -0x1.01259cp-1, -0x1.17fp-6, + 0x1.283p0, 0x1.5d1p3, + -0x1.6278p7, -0x1.c23p10, + 0x1.1444p13, 0x1.5fcp16 +}; + +// Lookup table for sin(k * pi / 32) with k = 0, ..., 63. +// Table is generated with Sollya as follows: +// > display = hexadecimmal; +// > prec = 23; +// > for k from 0 to 63 do {sin(k * pi/32);}; + +const float SIN_K_PI_OVER_32[64] = { + 0, 0x1.917a6cp-4, 0x1.8f8b84p-3, + 0x1.294064p-2, 0x1.87de2cp-2, 0x1.e2b5d4p-2, + 0x1.1c73b4p-1, 0x1.44cf34p-1, 0x1.6a09e8p-1, + 0x1.8bc808p-1, 0x1.a9b664p-1, 0x1.c38b3p-1, + 0x1.d906bcp-1, 0x1.e9f414p-1, 0x1.f6297cp-1, + 0x1.fd88dcp-1, 0x1p0, 0x1.fd88dcp-1, + 0x1.f6297cp-1, 0x1.e9f414p-1, 0x1.d906bcp-1, + 0x1.c38b3p-1, 0x1.a9b664p-1, 0x1.8bc808p-1, + 0x1.6a09e8p-1, 0x1.44cf34p-1, 0x1.1c73b4p-1, + 0x1.e2b5d4p-2, 0x1.87de2cp-2, 0x1.294064p-2, + 0x1.8f8b84p-3, 0x1.917a6cp-4, 0, + -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294064p-2, + -0x1.87de2cp-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, + -0x1.44cf34p-1, -0x1.6a09e8p-1, -0x1.8bc808p-1, + -0x1.a9b664p-1, -0x1.c38b3p-1, -0x1.d906bcp-1, + -0x1.e9f414p-1, -0x1.f6297cp-1, -0x1.fd88dcp-1, + -0x1p0, -0x1.fd88dcp-1, -0x1.f6297cp-1, + -0x1.e9f414p-1, -0x1.d906bcp-1, -0x1.c38b3p-1, + -0x1.a9b664p-1, -0x1.8bc808p-1, -0x1.6a09e8p-1, + -0x1.44cf34p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2, + -0x1.87de2cp-2, -0x1.294064p-2, -0x1.8f8b84p-3, + -0x1.917a6cp-4 +}; + +// horner's algorithm to accurately and efficiently evaluate a degree-9 +// polynomial iteratively +float horners(float x, const float COEFF[10]) { + float b8 = fputil::multiply_add(COEFF[9], x, COEFF[8]); + float b7 = fputil::multiply_add(b8, x, COEFF[7]); + float b6 = fputil::multiply_add(b7, x, COEFF[6]); + float b5 = fputil::multiply_add(b6, x, COEFF[5]); + float b4 = fputil::multiply_add(b5, x, COEFF[4]); + float b3 = fputil::multiply_add(b4, x, COEFF[3]); + float b2 = fputil::multiply_add(b3, x, COEFF[2]); + float b1 = fputil::multiply_add(b2, x, COEFF[1]); + return fputil::multiply_add(b1, x, COEFF[0]); +} + +float 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); +} +// HELPER_END + +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; + + // 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. + + 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]; + + float cos_y, sin_y; + if (y == 0) { + cos_y = 1; + sin_y = 0; + } else { + cos_y = horners(y * PI_OVER_32, COS_COEFF); + sin_y = horners(y * PI_OVER_32, SIN_COEFF); + } + + return static_cast(fputil::multiply_add(sin_k, cos_y, fputil::multiply_add(sin_y, cos_k, 0))); +} +} diff --git a/libc/src/math/sinpif16.h b/libc/src/math/sinpif16.h new file mode 100644 index 0000000000000..ce091c1737145 --- /dev/null +++ b/libc/src/math/sinpif16.h @@ -0,0 +1,23 @@ +//===-- Implementation header for sinpif16 ---------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache Licese 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_SINPIF16_H +#define LLVM_LIBC_SRC_MATH_SINPIF16_H + + +#include "include/llvm-libc-macros/float16-macros.h" +#include "src/__support/macros/config.h" +#include "src/__support/macros/properties/types.h" + +namespace LIBC_NAMESPACE_DECL { + +float16 sinpif16(float16 x); + +} // namespace LIBC_NAMESPACE_DECL + +#endif //LLVM_LIBC_SRC_MATH_SINPIF16_H diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index 47e16926f10df..33d03710abd03 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -51,6 +51,19 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + sinpif16_test + SUITE + libc-math-smoke-tests + SRCS + sinpif16_test.cpp + DEPENDS + libc.src.math.sinpif16 + libc.src.errno.errno + libc.src.__support.CPP.array + libc.src.__support.FPUtil.fp_bits +) + add_fp_unittest( sincosf_test SUITE diff --git a/libc/test/src/math/smoke/sinpif16_test.cpp b/libc/test/src/math/smoke/sinpif16_test.cpp new file mode 100644 index 0000000000000..e1356f493a103 --- /dev/null +++ b/libc/test/src/math/smoke/sinpif16_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/math/sinpif16.h" +#include "test/UnitTest/FPMatcher.h" +#include "src/errno/libc_errno.h" + +#include + +using LlvmLibcSinpif16Test = LIBC_NAMESPACE::testing::FPTest; + +TEST_F(LlvmLibcSinpif16Test, SpecialNumbers) { + LIBC_NAMESPACE::libc_errno = 0; + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif16(aNaN)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::sinpif16(0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(-0.0f, LIBC_NAMESPACE::sinpif16(-0.0f)); + EXPECT_MATH_ERRNO(0); + + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif16(inf)); + EXPECT_MATH_ERRNO(EDOM); + + EXPECT_FP_EQ(aNan, LIBC_NAMESPACE::sinpif16(neg_inf)); + EXPECT_MATH_ERRNO(EDOM); +} + +TEST_F(LlvmLibcSinpif16Test, Integers) { + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x420)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x1p+43)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x1.4p+64)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x420)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x1.cp+106)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x1.cp+21)); +} diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h index 8d51fa4e47726..00c298368e5a9 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -55,6 +55,7 @@ enum class Operation : int { RoundEven, Sin, Sinpi, + Sinpif16 Sinh, Sqrt, Tan, From d7bbbc20d1509e989a88bcb9cfd9a345dc4a2662 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 23 Sep 2024 23:30:36 +0100 Subject: [PATCH 02/20] clang format changes --- libc/src/math/generic/sinpif16.cpp | 127 +++++++++++++-------- libc/src/math/sinpif16.h | 3 +- libc/test/src/math/smoke/sinpif16_test.cpp | 2 +- libc/utils/MPFRWrapper/MPFRUtils.h | 3 +- 4 files changed, 84 insertions(+), 51 deletions(-) diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index bb028d5b88fc0..dc33188d5a5d5 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -15,7 +15,7 @@ #include "src/__support/macros/config.h" // TODO: Should probably create a new file; sincospif16_utils.h -// To store the following helper functions and constants. +// To store the following helper functions and constants. // I'd defer to @lntue for suggestions regarding that // HELPER_START @@ -23,19 +23,15 @@ namespace LIBC_NAMESPACE_DECL { constexpr float PI_OVER_32 = M_PI / 32; -// In Sollya generate 10 coeffecients for a degree-9 chebyshev polynomial -// approximating the sine function in [-pi / 32, pi / 32] with the following +// In Sollya generate 10 coeffecients for a degree-9 chebyshev polynomial +// approximating the sine function in [-pi / 32, pi / 32] with the following // commands: // > prec=23; // > TL = chebyshevform(sin(x), 9, [-pi / 32, pi / 32]); // > TL[0]; const float SIN_COEFF[10] = { - 0x1.801p-27, 0x1.000078p0, - -0x1.7e98p-14, -0x1.6bf4p-3, - 0x1.95ccp-5, 0x1.1baep2, - -0x1.030ap3, -0x1.3dap9, - 0x1.98e4p8, 0x1.d3d8p14 -}; + 0x1.801p-27, 0x1.000078p0, -0x1.7e98p-14, -0x1.6bf4p-3, 0x1.95ccp-5, + 0x1.1baep2, -0x1.030ap3, -0x1.3dap9, 0x1.98e4p8, 0x1.d3d8p14}; // In Sollya generate 10 coefficients for a degree-9 chebyshev polynomial // approximating the sine function in [-pi/32, pi/32] with the following @@ -44,12 +40,8 @@ const float SIN_COEFF[10] = { // > TL = chebyshevform(cos(x), 9, [-pi / 32, pi / 32]); // > TL[0]; const float COS_COEFF[10] = { - 0x1.00001p0, -0x1.48p-17, - -0x1.01259cp-1, -0x1.17fp-6, - 0x1.283p0, 0x1.5d1p3, - -0x1.6278p7, -0x1.c23p10, - 0x1.1444p13, 0x1.5fcp16 -}; + 0x1.00001p0, -0x1.48p-17, -0x1.01259cp-1, -0x1.17fp-6, 0x1.283p0, + 0x1.5d1p3, -0x1.6278p7, -0x1.c23p10, 0x1.1444p13, 0x1.5fcp16}; // Lookup table for sin(k * pi / 32) with k = 0, ..., 63. // Table is generated with Sollya as follows: @@ -57,30 +49,70 @@ const float COS_COEFF[10] = { // > prec = 23; // > for k from 0 to 63 do {sin(k * pi/32);}; -const float SIN_K_PI_OVER_32[64] = { - 0, 0x1.917a6cp-4, 0x1.8f8b84p-3, - 0x1.294064p-2, 0x1.87de2cp-2, 0x1.e2b5d4p-2, - 0x1.1c73b4p-1, 0x1.44cf34p-1, 0x1.6a09e8p-1, - 0x1.8bc808p-1, 0x1.a9b664p-1, 0x1.c38b3p-1, - 0x1.d906bcp-1, 0x1.e9f414p-1, 0x1.f6297cp-1, - 0x1.fd88dcp-1, 0x1p0, 0x1.fd88dcp-1, - 0x1.f6297cp-1, 0x1.e9f414p-1, 0x1.d906bcp-1, - 0x1.c38b3p-1, 0x1.a9b664p-1, 0x1.8bc808p-1, - 0x1.6a09e8p-1, 0x1.44cf34p-1, 0x1.1c73b4p-1, - 0x1.e2b5d4p-2, 0x1.87de2cp-2, 0x1.294064p-2, - 0x1.8f8b84p-3, 0x1.917a6cp-4, 0, - -0x1.917a6cp-4, -0x1.8f8b84p-3, -0x1.294064p-2, - -0x1.87de2cp-2, -0x1.e2b5d4p-2, -0x1.1c73b4p-1, - -0x1.44cf34p-1, -0x1.6a09e8p-1, -0x1.8bc808p-1, - -0x1.a9b664p-1, -0x1.c38b3p-1, -0x1.d906bcp-1, - -0x1.e9f414p-1, -0x1.f6297cp-1, -0x1.fd88dcp-1, - -0x1p0, -0x1.fd88dcp-1, -0x1.f6297cp-1, - -0x1.e9f414p-1, -0x1.d906bcp-1, -0x1.c38b3p-1, - -0x1.a9b664p-1, -0x1.8bc808p-1, -0x1.6a09e8p-1, - -0x1.44cf34p-1, -0x1.1c73b4p-1, -0x1.e2b5d4p-2, - -0x1.87de2cp-2, -0x1.294064p-2, -0x1.8f8b84p-3, - -0x1.917a6cp-4 -}; +const float SIN_K_PI_OVER_32[64] = {0, + 0x1.917a6cp-4, + 0x1.8f8b84p-3, + 0x1.294064p-2, + 0x1.87de2cp-2, + 0x1.e2b5d4p-2, + 0x1.1c73b4p-1, + 0x1.44cf34p-1, + 0x1.6a09e8p-1, + 0x1.8bc808p-1, + 0x1.a9b664p-1, + 0x1.c38b3p-1, + 0x1.d906bcp-1, + 0x1.e9f414p-1, + 0x1.f6297cp-1, + 0x1.fd88dcp-1, + 0x1p0, + 0x1.fd88dcp-1, + 0x1.f6297cp-1, + 0x1.e9f414p-1, + 0x1.d906bcp-1, + 0x1.c38b3p-1, + 0x1.a9b664p-1, + 0x1.8bc808p-1, + 0x1.6a09e8p-1, + 0x1.44cf34p-1, + 0x1.1c73b4p-1, + 0x1.e2b5d4p-2, + 0x1.87de2cp-2, + 0x1.294064p-2, + 0x1.8f8b84p-3, + 0x1.917a6cp-4, + 0, + -0x1.917a6cp-4, + -0x1.8f8b84p-3, + -0x1.294064p-2, + -0x1.87de2cp-2, + -0x1.e2b5d4p-2, + -0x1.1c73b4p-1, + -0x1.44cf34p-1, + -0x1.6a09e8p-1, + -0x1.8bc808p-1, + -0x1.a9b664p-1, + -0x1.c38b3p-1, + -0x1.d906bcp-1, + -0x1.e9f414p-1, + -0x1.f6297cp-1, + -0x1.fd88dcp-1, + -0x1p0, + -0x1.fd88dcp-1, + -0x1.f6297cp-1, + -0x1.e9f414p-1, + -0x1.d906bcp-1, + -0x1.c38b3p-1, + -0x1.a9b664p-1, + -0x1.8bc808p-1, + -0x1.6a09e8p-1, + -0x1.44cf34p-1, + -0x1.1c73b4p-1, + -0x1.e2b5d4p-2, + -0x1.87de2cp-2, + -0x1.294064p-2, + -0x1.8f8b84p-3, + -0x1.917a6cp-4}; // horner's algorithm to accurately and efficiently evaluate a degree-9 // polynomial iteratively @@ -96,7 +128,7 @@ float horners(float x, const float COEFF[10]) { return fputil::multiply_add(b1, x, COEFF[0]); } -float range_reduction(float x, float& y) { +float range_reduction(float x, float &y) { float kf = fputil::nearest_integer(x * 32); y = fputil::multiply_add(x, 32.0, -kf); @@ -125,10 +157,12 @@ 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) + // 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. float f32 = x; float y; @@ -146,6 +180,7 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { sin_y = horners(y * PI_OVER_32, SIN_COEFF); } - return static_cast(fputil::multiply_add(sin_k, cos_y, fputil::multiply_add(sin_y, cos_k, 0))); -} + return static_cast(fputil::multiply_add( + sin_k, cos_y, fputil::multiply_add(sin_y, cos_k, 0))); } +} // namespace LIBC_NAMESPACE_DECL diff --git a/libc/src/math/sinpif16.h b/libc/src/math/sinpif16.h index ce091c1737145..f7c653383ba5e 100644 --- a/libc/src/math/sinpif16.h +++ b/libc/src/math/sinpif16.h @@ -9,7 +9,6 @@ #ifndef LLVM_LIBC_SRC_MATH_SINPIF16_H #define LLVM_LIBC_SRC_MATH_SINPIF16_H - #include "include/llvm-libc-macros/float16-macros.h" #include "src/__support/macros/config.h" #include "src/__support/macros/properties/types.h" @@ -20,4 +19,4 @@ float16 sinpif16(float16 x); } // namespace LIBC_NAMESPACE_DECL -#endif //LLVM_LIBC_SRC_MATH_SINPIF16_H +#endif // LLVM_LIBC_SRC_MATH_SINPIF16_H diff --git a/libc/test/src/math/smoke/sinpif16_test.cpp b/libc/test/src/math/smoke/sinpif16_test.cpp index e1356f493a103..54eb12ae100ef 100644 --- a/libc/test/src/math/smoke/sinpif16_test.cpp +++ b/libc/test/src/math/smoke/sinpif16_test.cpp @@ -7,9 +7,9 @@ // // ===----------------------------------------------------------------------==// +#include "src/errno/libc_errno.h" #include "src/math/sinpif16.h" #include "test/UnitTest/FPMatcher.h" -#include "src/errno/libc_errno.h" #include diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h index 00c298368e5a9..43a30ad87c57a 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -55,8 +55,7 @@ enum class Operation : int { RoundEven, Sin, Sinpi, - Sinpif16 - Sinh, + Sinpif16 Sinh, Sqrt, Tan, Tanh, From 663fbd51d8992daf763772bc61b5942a5a1e8249 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Thu, 3 Oct 2024 14:42:57 +0100 Subject: [PATCH 03/20] add implementation of sinpif16 function --- libc/src/math/generic/sinpif16.cpp | 181 ++++++++++----------- libc/test/src/math/smoke/sinpif16_test.cpp | 12 +- libc/utils/MPFRWrapper/MPFRUtils.h | 2 +- 3 files changed, 91 insertions(+), 104 deletions(-) diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index dc33188d5a5d5..b878e09d54200 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -5,12 +5,12 @@ // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //===----------------------------------------------------------------------===// - -#define M_PI 3.1415925f - #include "src/math/sinpif16.h" +#include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" +#include "src/__support/FPUtil/PolyEval.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" @@ -21,114 +21,70 @@ // HELPER_START namespace LIBC_NAMESPACE_DECL { -constexpr float PI_OVER_32 = M_PI / 32; +constexpr float PI_OVER_32 = 0x1.921fb6p-4f; // In Sollya generate 10 coeffecients for a degree-9 chebyshev polynomial // approximating the sine function in [-pi / 32, pi / 32] with the following // commands: -// > prec=23; +// > prec=24; // > TL = chebyshevform(sin(x), 9, [-pi / 32, pi / 32]); // > TL[0]; const float SIN_COEFF[10] = { - 0x1.801p-27, 0x1.000078p0, -0x1.7e98p-14, -0x1.6bf4p-3, 0x1.95ccp-5, - 0x1.1baep2, -0x1.030ap3, -0x1.3dap9, 0x1.98e4p8, 0x1.d3d8p14}; - + 0x1.d333p-26, 0x1.000048p0, -0x1.a5d2p-14, -0x1.628588p-3, 0x1.c1eep-5, + 0x1.4455p1, -0x1.317a8p3, -0x1.6bb9p8, 0x1.00ef8p9, 0x1.0edcp14 +}; // In Sollya generate 10 coefficients for a degree-9 chebyshev polynomial // approximating the sine function in [-pi/32, pi/32] with the following // commands: -// > prec = 23; +// > prec = 24; // > TL = chebyshevform(cos(x), 9, [-pi / 32, pi / 32]); // > TL[0]; const float COS_COEFF[10] = { - 0x1.00001p0, -0x1.48p-17, -0x1.01259cp-1, -0x1.17fp-6, 0x1.283p0, - 0x1.5d1p3, -0x1.6278p7, -0x1.c23p10, 0x1.1444p13, 0x1.5fcp16}; - + 0x1.000006p0, 0x1.e1eap-15, -0x1.0071p-1, -0x1.3b56p-4, 0x1.f3dfp-2, + 0x1.ccbap4, -0x1.3034p6, -0x1.f817p11, 0x1.fc59p11, 0x1.7079p17 +}; // Lookup table for sin(k * pi / 32) with k = 0, ..., 63. // Table is generated with Sollya as follows: // > display = hexadecimmal; -// > prec = 23; +// > prec = 24; // > for k from 0 to 63 do {sin(k * pi/32);}; -const float SIN_K_PI_OVER_32[64] = {0, - 0x1.917a6cp-4, - 0x1.8f8b84p-3, - 0x1.294064p-2, - 0x1.87de2cp-2, - 0x1.e2b5d4p-2, - 0x1.1c73b4p-1, - 0x1.44cf34p-1, - 0x1.6a09e8p-1, - 0x1.8bc808p-1, - 0x1.a9b664p-1, - 0x1.c38b3p-1, - 0x1.d906bcp-1, - 0x1.e9f414p-1, - 0x1.f6297cp-1, - 0x1.fd88dcp-1, - 0x1p0, - 0x1.fd88dcp-1, - 0x1.f6297cp-1, - 0x1.e9f414p-1, - 0x1.d906bcp-1, - 0x1.c38b3p-1, - 0x1.a9b664p-1, - 0x1.8bc808p-1, - 0x1.6a09e8p-1, - 0x1.44cf34p-1, - 0x1.1c73b4p-1, - 0x1.e2b5d4p-2, - 0x1.87de2cp-2, - 0x1.294064p-2, - 0x1.8f8b84p-3, - 0x1.917a6cp-4, - 0, - -0x1.917a6cp-4, - -0x1.8f8b84p-3, - -0x1.294064p-2, - -0x1.87de2cp-2, - -0x1.e2b5d4p-2, - -0x1.1c73b4p-1, - -0x1.44cf34p-1, - -0x1.6a09e8p-1, - -0x1.8bc808p-1, - -0x1.a9b664p-1, - -0x1.c38b3p-1, - -0x1.d906bcp-1, - -0x1.e9f414p-1, - -0x1.f6297cp-1, - -0x1.fd88dcp-1, - -0x1p0, - -0x1.fd88dcp-1, - -0x1.f6297cp-1, - -0x1.e9f414p-1, - -0x1.d906bcp-1, - -0x1.c38b3p-1, - -0x1.a9b664p-1, - -0x1.8bc808p-1, - -0x1.6a09e8p-1, - -0x1.44cf34p-1, - -0x1.1c73b4p-1, - -0x1.e2b5d4p-2, - -0x1.87de2cp-2, - -0x1.294064p-2, - -0x1.8f8b84p-3, - -0x1.917a6cp-4}; - -// horner's algorithm to accurately and efficiently evaluate a degree-9 -// polynomial iteratively -float horners(float x, const float COEFF[10]) { - float b8 = fputil::multiply_add(COEFF[9], x, COEFF[8]); - float b7 = fputil::multiply_add(b8, x, COEFF[7]); - float b6 = fputil::multiply_add(b7, x, COEFF[6]); - float b5 = fputil::multiply_add(b6, x, COEFF[5]); - float b4 = fputil::multiply_add(b5, x, COEFF[4]); - float b3 = fputil::multiply_add(b4, x, COEFF[3]); - float b2 = fputil::multiply_add(b3, x, COEFF[2]); - float b1 = fputil::multiply_add(b2, x, COEFF[1]); - return fputil::multiply_add(b1, x, COEFF[0]); -} - -float range_reduction(float x, float &y) { +const float SIN_K_PI_OVER_32[64] = { + 0, 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, + 0, -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 +}; + +int32_t range_reduction(float x, float &y) { float kf = fputil::nearest_integer(x * 32); y = fputil::multiply_add(x, 32.0, -kf); @@ -164,7 +120,28 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev // polynomials generated by Sollya. - float f32 = x; + if (LIBC_UNLIKELY(x_abs == 0U)) { + // For signed zeros + return x; + } + + // Numbers greater or equal to 2^10 are integers or NaN + if (LIBC_UNLIKELY(x_abs >= 0x6400)) { + // Check for NaN or infinity 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); + } + + // If value is NaN + return x + FPBits::quiet_nan().get_val(); + } + return FPBits::zero(xbits.sign()).get_val(); + } + + float f32 = static_cast(x); float y; int32_t k = range_reduction(f32, y); @@ -176,11 +153,21 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { cos_y = 1; sin_y = 0; } else { - cos_y = horners(y * PI_OVER_32, COS_COEFF); - sin_y = horners(y * PI_OVER_32, SIN_COEFF); + cos_y = fputil::polyeval(y * PI_OVER_32, + COS_COEFF[0], COS_COEFF[1], + COS_COEFF[2], COS_COEFF[3], + COS_COEFF[4], COS_COEFF[5], + COS_COEFF[6], COS_COEFF[7], + COS_COEFF[8], COS_COEFF[9]); + sin_y = fputil::polyeval(y * PI_OVER_32, + SIN_COEFF[0], SIN_COEFF[1], + SIN_COEFF[2], SIN_COEFF[3], + SIN_COEFF[4], SIN_COEFF[5], + SIN_COEFF[6], SIN_COEFF[7], + SIN_COEFF[8], SIN_COEFF[9]); } return static_cast(fputil::multiply_add( - sin_k, cos_y, fputil::multiply_add(sin_y, cos_k, 0))); + sin_k, cos_y, fputil::multiply_add(sin_y, cos_k, 0.0f))); } } // namespace LIBC_NAMESPACE_DECL diff --git a/libc/test/src/math/smoke/sinpif16_test.cpp b/libc/test/src/math/smoke/sinpif16_test.cpp index 54eb12ae100ef..bf55a4854bfb3 100644 --- a/libc/test/src/math/smoke/sinpif16_test.cpp +++ b/libc/test/src/math/smoke/sinpif16_test.cpp @@ -1,4 +1,4 @@ -//===-- Unittests for sinpif16 ------------------------------------===// +//===-- 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. @@ -30,15 +30,15 @@ TEST_F(LlvmLibcSinpif16Test, SpecialNumbers) { EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif16(inf)); EXPECT_MATH_ERRNO(EDOM); - EXPECT_FP_EQ(aNan, LIBC_NAMESPACE::sinpif16(neg_inf)); + EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif16(neg_inf)); EXPECT_MATH_ERRNO(EDOM); } TEST_F(LlvmLibcSinpif16Test, Integers) { EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x420)); - EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x1p+43)); - EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x1.4p+64)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x1p+10)); + EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x1.4p+14)); EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x420)); - EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x1.cp+106)); - EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x1.cp+21)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x1.cp+15)); + EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x1.cp+7)); } diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h index 43a30ad87c57a..8d51fa4e47726 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.h +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -55,7 +55,7 @@ enum class Operation : int { RoundEven, Sin, Sinpi, - Sinpif16 Sinh, + Sinh, Sqrt, Tan, Tanh, From 0f038fed86a108336958adb5b17bf00b04f36d23 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Thu, 3 Oct 2024 15:01:10 +0100 Subject: [PATCH 04/20] clang format --- libc/src/math/generic/sinpif16.cpp | 94 ++++++++++++++++++++++++++++++ 1 file changed, 94 insertions(+) diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index b878e09d54200..7dd423468fc80 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -30,9 +30,14 @@ constexpr float PI_OVER_32 = 0x1.921fb6p-4f; // > TL = chebyshevform(sin(x), 9, [-pi / 32, pi / 32]); // > TL[0]; const float SIN_COEFF[10] = { +<<<<<<< HEAD 0x1.d333p-26, 0x1.000048p0, -0x1.a5d2p-14, -0x1.628588p-3, 0x1.c1eep-5, 0x1.4455p1, -0x1.317a8p3, -0x1.6bb9p8, 0x1.00ef8p9, 0x1.0edcp14 }; +======= + 0x1.d333p-26, 0x1.000048p0, -0x1.a5d2p-14, -0x1.628588p-3, 0x1.c1eep-5, + 0x1.4455p1, -0x1.317a8p3, -0x1.6bb9p8, 0x1.00ef8p9, 0x1.0edcp14}; +>>>>>>> 1e0a1f667 (add implementation of sinpif16 function) // In Sollya generate 10 coefficients for a degree-9 chebyshev polynomial // approximating the sine function in [-pi/32, pi/32] with the following // commands: @@ -40,15 +45,21 @@ const float SIN_COEFF[10] = { // > TL = chebyshevform(cos(x), 9, [-pi / 32, pi / 32]); // > TL[0]; const float COS_COEFF[10] = { +<<<<<<< HEAD 0x1.000006p0, 0x1.e1eap-15, -0x1.0071p-1, -0x1.3b56p-4, 0x1.f3dfp-2, 0x1.ccbap4, -0x1.3034p6, -0x1.f817p11, 0x1.fc59p11, 0x1.7079p17 }; +======= + 0x1.000006p0, 0x1.e1eap-15, -0x1.0071p-1, -0x1.3b56p-4, 0x1.f3dfp-2, + 0x1.ccbap4, -0x1.3034p6, -0x1.f817p11, 0x1.fc59p11, 0x1.7079p17}; +>>>>>>> 1e0a1f667 (add implementation of sinpif16 function) // Lookup table for sin(k * pi / 32) with k = 0, ..., 63. // Table is generated with Sollya as follows: // > display = hexadecimmal; // > prec = 24; // > for k from 0 to 63 do {sin(k * pi/32);}; +<<<<<<< HEAD const float SIN_K_PI_OVER_32[64] = { 0, 0x1.917a6cp-4, 0x1.8f8b84p-3, 0x1.294062p-2, @@ -83,6 +94,72 @@ const float SIN_K_PI_OVER_32[64] = { -0x1.87de2ap-2, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4 }; +======= +const float SIN_K_PI_OVER_32[64] = {0, + 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, + 0, + -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}; +>>>>>>> 1e0a1f667 (add implementation of sinpif16 function) int32_t range_reduction(float x, float &y) { float kf = fputil::nearest_integer(x * 32); @@ -132,9 +209,15 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // If value is equal to infinity if (x_abs == 0x7c00) { fputil::set_errno_if_required(EDOM); +<<<<<<< HEAD fputil::raise_except_if_required(FE_INVALID); } +======= + fputil::raise_except_if_required(FE_INVALID); + } + +>>>>>>> 1e0a1f667 (add implementation of sinpif16 function) // If value is NaN return x + FPBits::quiet_nan().get_val(); } @@ -153,6 +236,7 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { cos_y = 1; sin_y = 0; } else { +<<<<<<< HEAD cos_y = fputil::polyeval(y * PI_OVER_32, COS_COEFF[0], COS_COEFF[1], COS_COEFF[2], COS_COEFF[3], @@ -165,6 +249,16 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { SIN_COEFF[4], SIN_COEFF[5], SIN_COEFF[6], SIN_COEFF[7], SIN_COEFF[8], SIN_COEFF[9]); +======= + cos_y = fputil::polyeval(y * PI_OVER_32, COS_COEFF[0], COS_COEFF[1], + COS_COEFF[2], COS_COEFF[3], COS_COEFF[4], + COS_COEFF[5], COS_COEFF[6], COS_COEFF[7], + COS_COEFF[8], COS_COEFF[9]); + sin_y = fputil::polyeval(y * PI_OVER_32, SIN_COEFF[0], SIN_COEFF[1], + SIN_COEFF[2], SIN_COEFF[3], SIN_COEFF[4], + SIN_COEFF[5], SIN_COEFF[6], SIN_COEFF[7], + SIN_COEFF[8], SIN_COEFF[9]); +>>>>>>> 1e0a1f667 (add implementation of sinpif16 function) } return static_cast(fputil::multiply_add( From c54defca1bcbd6146a9971651b6b78db0d667122 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Thu, 3 Oct 2024 15:03:29 +0100 Subject: [PATCH 05/20] clang format --- libc/src/math/generic/sinpif16.cpp | 70 ------------------------------ 1 file changed, 70 deletions(-) diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 7dd423468fc80..eff3871a3bfaa 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -30,14 +30,8 @@ constexpr float PI_OVER_32 = 0x1.921fb6p-4f; // > TL = chebyshevform(sin(x), 9, [-pi / 32, pi / 32]); // > TL[0]; const float SIN_COEFF[10] = { -<<<<<<< HEAD - 0x1.d333p-26, 0x1.000048p0, -0x1.a5d2p-14, -0x1.628588p-3, 0x1.c1eep-5, - 0x1.4455p1, -0x1.317a8p3, -0x1.6bb9p8, 0x1.00ef8p9, 0x1.0edcp14 -}; -======= 0x1.d333p-26, 0x1.000048p0, -0x1.a5d2p-14, -0x1.628588p-3, 0x1.c1eep-5, 0x1.4455p1, -0x1.317a8p3, -0x1.6bb9p8, 0x1.00ef8p9, 0x1.0edcp14}; ->>>>>>> 1e0a1f667 (add implementation of sinpif16 function) // In Sollya generate 10 coefficients for a degree-9 chebyshev polynomial // approximating the sine function in [-pi/32, pi/32] with the following // commands: @@ -45,56 +39,14 @@ const float SIN_COEFF[10] = { // > TL = chebyshevform(cos(x), 9, [-pi / 32, pi / 32]); // > TL[0]; const float COS_COEFF[10] = { -<<<<<<< HEAD - 0x1.000006p0, 0x1.e1eap-15, -0x1.0071p-1, -0x1.3b56p-4, 0x1.f3dfp-2, - 0x1.ccbap4, -0x1.3034p6, -0x1.f817p11, 0x1.fc59p11, 0x1.7079p17 -}; -======= 0x1.000006p0, 0x1.e1eap-15, -0x1.0071p-1, -0x1.3b56p-4, 0x1.f3dfp-2, 0x1.ccbap4, -0x1.3034p6, -0x1.f817p11, 0x1.fc59p11, 0x1.7079p17}; ->>>>>>> 1e0a1f667 (add implementation of sinpif16 function) // Lookup table for sin(k * pi / 32) with k = 0, ..., 63. // Table is generated with Sollya as follows: // > display = hexadecimmal; // > prec = 24; // > for k from 0 to 63 do {sin(k * pi/32);}; -<<<<<<< HEAD -const float SIN_K_PI_OVER_32[64] = { - 0, 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, - 0, -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 -}; -======= const float SIN_K_PI_OVER_32[64] = {0, 0x1.917a6cp-4, 0x1.8f8b84p-3, @@ -159,7 +111,6 @@ const float SIN_K_PI_OVER_32[64] = {0, -0x1.294062p-2, -0x1.8f8b84p-3, -0x1.917a6cp-4}; ->>>>>>> 1e0a1f667 (add implementation of sinpif16 function) int32_t range_reduction(float x, float &y) { float kf = fputil::nearest_integer(x * 32); @@ -209,15 +160,9 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // If value is equal to infinity if (x_abs == 0x7c00) { fputil::set_errno_if_required(EDOM); -<<<<<<< HEAD - fputil::raise_except_if_required(FE_INVALID); - } - -======= fputil::raise_except_if_required(FE_INVALID); } ->>>>>>> 1e0a1f667 (add implementation of sinpif16 function) // If value is NaN return x + FPBits::quiet_nan().get_val(); } @@ -236,20 +181,6 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { cos_y = 1; sin_y = 0; } else { -<<<<<<< HEAD - cos_y = fputil::polyeval(y * PI_OVER_32, - COS_COEFF[0], COS_COEFF[1], - COS_COEFF[2], COS_COEFF[3], - COS_COEFF[4], COS_COEFF[5], - COS_COEFF[6], COS_COEFF[7], - COS_COEFF[8], COS_COEFF[9]); - sin_y = fputil::polyeval(y * PI_OVER_32, - SIN_COEFF[0], SIN_COEFF[1], - SIN_COEFF[2], SIN_COEFF[3], - SIN_COEFF[4], SIN_COEFF[5], - SIN_COEFF[6], SIN_COEFF[7], - SIN_COEFF[8], SIN_COEFF[9]); -======= cos_y = fputil::polyeval(y * PI_OVER_32, COS_COEFF[0], COS_COEFF[1], COS_COEFF[2], COS_COEFF[3], COS_COEFF[4], COS_COEFF[5], COS_COEFF[6], COS_COEFF[7], @@ -258,7 +189,6 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { SIN_COEFF[2], SIN_COEFF[3], SIN_COEFF[4], SIN_COEFF[5], SIN_COEFF[6], SIN_COEFF[7], SIN_COEFF[8], SIN_COEFF[9]); ->>>>>>> 1e0a1f667 (add implementation of sinpif16 function) } return static_cast(fputil::multiply_add( From a8f7484896251b9c94accb5ff9aad06f4c7b5b6a Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 7 Oct 2024 13:38:19 +0100 Subject: [PATCH 06/20] nit changes and review suggestions --- libc/config/linux/aarch64/entrypoints.txt | 2 +- libc/config/linux/arm/entrypoints.txt | 1 - libc/config/linux/riscv/entrypoints.txt | 1 - libc/config/windows/entrypoints.txt | 3 +-- libc/newhdrgen/yaml/math.yaml | 1 + libc/src/math/generic/CMakeLists.txt | 10 ++++++---- libc/src/math/generic/sinpif16.cpp | 19 +++++++------------ libc/src/math/sinpif16.h | 1 - libc/test/src/math/smoke/CMakeLists.txt | 4 +--- libc/test/src/math/smoke/sinpif16_test.cpp | 3 +-- 10 files changed, 18 insertions(+), 27 deletions(-) diff --git a/libc/config/linux/aarch64/entrypoints.txt b/libc/config/linux/aarch64/entrypoints.txt index 85ae7be80b461..c439efcd3211f 100644 --- a/libc/config/linux/aarch64/entrypoints.txt +++ b/libc/config/linux/aarch64/entrypoints.txt @@ -570,7 +570,6 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sinf libc.src.math.sinhf libc.src.math.sinpif - libc.src.math.sinpif16 libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl @@ -672,6 +671,7 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.scalbnf16 libc.src.math.setpayloadf16 libc.src.math.setpayloadsigf16 + libc.src.math.sinpif16 libc.src.math.totalorderf16 libc.src.math.totalordermagf16 libc.src.math.truncf16 diff --git a/libc/config/linux/arm/entrypoints.txt b/libc/config/linux/arm/entrypoints.txt index dbd21f6b5d438..1be9a872dd2f7 100644 --- a/libc/config/linux/arm/entrypoints.txt +++ b/libc/config/linux/arm/entrypoints.txt @@ -390,7 +390,6 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sincosf libc.src.math.sinf libc.src.math.sinhf - libc.src.math.sinpif16 libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/config/linux/riscv/entrypoints.txt b/libc/config/linux/riscv/entrypoints.txt index 8b71f44cf9354..8312b2c453f23 100644 --- a/libc/config/linux/riscv/entrypoints.txt +++ b/libc/config/linux/riscv/entrypoints.txt @@ -573,7 +573,6 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sinf libc.src.math.sinhf libc.src.math.sinpif - libc.src.math.sinpif16 libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/config/windows/entrypoints.txt b/libc/config/windows/entrypoints.txt index 5b0d84669e924..8f0b50bcc83ea 100644 --- a/libc/config/windows/entrypoints.txt +++ b/libc/config/windows/entrypoints.txt @@ -275,8 +275,7 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sincosf libc.src.math.sincosf libc.src.math.sinf - libc.src.math.sinhfi - libc.src.math.sinpif16 + libc.src.math.sinhf libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl diff --git a/libc/newhdrgen/yaml/math.yaml b/libc/newhdrgen/yaml/math.yaml index fb79f15a31460..2c6656c2a7756 100644 --- a/libc/newhdrgen/yaml/math.yaml +++ b/libc/newhdrgen/yaml/math.yaml @@ -2296,6 +2296,7 @@ functions: return_type: _Float16 arguments: - type: _Float16 + guard: LIBC_TYPES_HAS_FLOAT16 - name: sqrt standards: - stdc diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index a87f041b55844..ce5ba59139a36 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -535,11 +535,13 @@ add_entrypoint_object( HDRS ../sinpif16.h DEPENDS - libc.src.__support.macros.properties.types - libc.src.__support.FPUtil.fp_bits - libc.src.__support.FPUtil.multiply_add libc.src.__support.common - libc.src.__support.macros.config + 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.poly_eval + libc.src.__support.macros.properties.types COMPILE_OPTIONS -O3 ) diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index eff3871a3bfaa..5ac541b97ec5d 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -5,6 +5,7 @@ // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // //===----------------------------------------------------------------------===// + #include "src/math/sinpif16.h" #include "src/__support/FPUtil/FEnvImpl.h" #include "src/__support/FPUtil/FPBits.h" @@ -14,11 +15,6 @@ #include "src/__support/common.h" #include "src/__support/macros/config.h" -// TODO: Should probably create a new file; sincospif16_utils.h -// To store the following helper functions and constants. -// I'd defer to @lntue for suggestions regarding that - -// HELPER_START namespace LIBC_NAMESPACE_DECL { constexpr float PI_OVER_32 = 0x1.921fb6p-4f; @@ -29,24 +25,25 @@ constexpr float PI_OVER_32 = 0x1.921fb6p-4f; // > prec=24; // > TL = chebyshevform(sin(x), 9, [-pi / 32, pi / 32]); // > TL[0]; -const float SIN_COEFF[10] = { +static constexpr float SIN_COEFF[10] = { 0x1.d333p-26, 0x1.000048p0, -0x1.a5d2p-14, -0x1.628588p-3, 0x1.c1eep-5, 0x1.4455p1, -0x1.317a8p3, -0x1.6bb9p8, 0x1.00ef8p9, 0x1.0edcp14}; + // In Sollya generate 10 coefficients for a degree-9 chebyshev polynomial // approximating the sine function in [-pi/32, pi/32] with the following // commands: // > prec = 24; // > TL = chebyshevform(cos(x), 9, [-pi / 32, pi / 32]); // > TL[0]; -const float COS_COEFF[10] = { +static constexpr float COS_COEFF[10] = { 0x1.000006p0, 0x1.e1eap-15, -0x1.0071p-1, -0x1.3b56p-4, 0x1.f3dfp-2, 0x1.ccbap4, -0x1.3034p6, -0x1.f817p11, 0x1.fc59p11, 0x1.7079p17}; + // Lookup table for sin(k * pi / 32) with k = 0, ..., 63. // Table is generated with Sollya as follows: // > display = hexadecimmal; // > prec = 24; // > for k from 0 to 63 do {sin(k * pi/32);}; - const float SIN_K_PI_OVER_32[64] = {0, 0x1.917a6cp-4, 0x1.8f8b84p-3, @@ -118,7 +115,6 @@ int32_t range_reduction(float x, float &y) { return static_cast(kf); } -// HELPER_END LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { using FPBits = typename fputil::FPBits; @@ -153,7 +149,7 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { return x; } - // Numbers greater or equal to 2^10 are integers or NaN + // Numbers greater or equal to 2^10 are integers, or infinity, or NaN if (LIBC_UNLIKELY(x_abs >= 0x6400)) { // Check for NaN or infinity values if (LIBC_UNLIKELY(x_abs >= 0x7c00)) { @@ -163,13 +159,12 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { fputil::raise_except_if_required(FE_INVALID); } - // If value is NaN return x + FPBits::quiet_nan().get_val(); } return FPBits::zero(xbits.sign()).get_val(); } - float f32 = static_cast(x); + float f32 = x; float y; int32_t k = range_reduction(f32, y); diff --git a/libc/src/math/sinpif16.h b/libc/src/math/sinpif16.h index f7c653383ba5e..33a0ae2658401 100644 --- a/libc/src/math/sinpif16.h +++ b/libc/src/math/sinpif16.h @@ -9,7 +9,6 @@ #ifndef LLVM_LIBC_SRC_MATH_SINPIF16_H #define LLVM_LIBC_SRC_MATH_SINPIF16_H -#include "include/llvm-libc-macros/float16-macros.h" #include "src/__support/macros/config.h" #include "src/__support/macros/properties/types.h" diff --git a/libc/test/src/math/smoke/CMakeLists.txt b/libc/test/src/math/smoke/CMakeLists.txt index a72ad128e21d1..67466c896f30c 100644 --- a/libc/test/src/math/smoke/CMakeLists.txt +++ b/libc/test/src/math/smoke/CMakeLists.txt @@ -58,10 +58,8 @@ add_fp_unittest( SRCS sinpif16_test.cpp DEPENDS - libc.src.math.sinpif16 libc.src.errno.errno - libc.src.__support.CPP.array - libc.src.__support.FPUtil.fp_bits + libc.src.math.sinpif16 ) add_fp_unittest( diff --git a/libc/test/src/math/smoke/sinpif16_test.cpp b/libc/test/src/math/smoke/sinpif16_test.cpp index bf55a4854bfb3..0cf6f7c6a86e2 100644 --- a/libc/test/src/math/smoke/sinpif16_test.cpp +++ b/libc/test/src/math/smoke/sinpif16_test.cpp @@ -10,8 +10,7 @@ #include "src/errno/libc_errno.h" #include "src/math/sinpif16.h" #include "test/UnitTest/FPMatcher.h" - -#include +#include "test/UnitTest/Test.h" using LlvmLibcSinpif16Test = LIBC_NAMESPACE::testing::FPTest; From 3f2282c0b869bf9a704fcf4d7729a19409bbfed4 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 7 Oct 2024 14:02:11 +0100 Subject: [PATCH 07/20] fix dependency misspell --- libc/src/math/generic/CMakeLists.txt | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index ce5ba59139a36..48078a158bfb9 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -536,11 +536,11 @@ add_entrypoint_object( ../sinpif16.h DEPENDS libc.src.__support.common - libc.src.__support.FPutil.fenv_impl + 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.poly_eval + libc.src.__support.FPUtil.polyeval libc.src.__support.macros.properties.types COMPILE_OPTIONS -O3 From 6617bf9a7cfceb725f8c94f0c3dc093649ae24e9 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Wed, 9 Oct 2024 09:18:45 +0100 Subject: [PATCH 08/20] updated function implementation to better handle extremely small non-zero fractional parts of the input x --- libc/src/math/generic/sinpif16.cpp | 68 +++++++++++++----------------- 1 file changed, 30 insertions(+), 38 deletions(-) diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 5ac541b97ec5d..0c96e6e173105 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -17,28 +17,6 @@ namespace LIBC_NAMESPACE_DECL { -constexpr float PI_OVER_32 = 0x1.921fb6p-4f; - -// In Sollya generate 10 coeffecients for a degree-9 chebyshev polynomial -// approximating the sine function in [-pi / 32, pi / 32] with the following -// commands: -// > prec=24; -// > TL = chebyshevform(sin(x), 9, [-pi / 32, pi / 32]); -// > TL[0]; -static constexpr float SIN_COEFF[10] = { - 0x1.d333p-26, 0x1.000048p0, -0x1.a5d2p-14, -0x1.628588p-3, 0x1.c1eep-5, - 0x1.4455p1, -0x1.317a8p3, -0x1.6bb9p8, 0x1.00ef8p9, 0x1.0edcp14}; - -// In Sollya generate 10 coefficients for a degree-9 chebyshev polynomial -// approximating the sine function in [-pi/32, pi/32] with the following -// commands: -// > prec = 24; -// > TL = chebyshevform(cos(x), 9, [-pi / 32, pi / 32]); -// > TL[0]; -static constexpr float COS_COEFF[10] = { - 0x1.000006p0, 0x1.e1eap-15, -0x1.0071p-1, -0x1.3b56p-4, 0x1.f3dfp-2, - 0x1.ccbap4, -0x1.3034p6, -0x1.f817p11, 0x1.fc59p11, 0x1.7079p17}; - // Lookup table for sin(k * pi / 32) with k = 0, ..., 63. // Table is generated with Sollya as follows: // > display = hexadecimmal; @@ -171,22 +149,36 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { float sin_k = SIN_K_PI_OVER_32[k & 63]; float cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; - float cos_y, sin_y; - if (y == 0) { - cos_y = 1; - sin_y = 0; - } else { - cos_y = fputil::polyeval(y * PI_OVER_32, COS_COEFF[0], COS_COEFF[1], - COS_COEFF[2], COS_COEFF[3], COS_COEFF[4], - COS_COEFF[5], COS_COEFF[6], COS_COEFF[7], - COS_COEFF[8], COS_COEFF[9]); - sin_y = fputil::polyeval(y * PI_OVER_32, SIN_COEFF[0], SIN_COEFF[1], - SIN_COEFF[2], SIN_COEFF[3], SIN_COEFF[4], - SIN_COEFF[5], SIN_COEFF[6], SIN_COEFF[7], - SIN_COEFF[8], SIN_COEFF[9]); - } - + float cosm1_y, sin_y; + + // 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); + + // Degree-6 minimax even polynomial for cos(y*pi/32) generated by Sollya with: + // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, SG...|], [0, 0.5]); + // 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) + cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, + 0x1.a6f7a2p-29f); + + // Since, cosm1_y = cos_y - 1, therefore: + // sin(x * pi) = cos_k * sin_y + sin_k + (cosm1_y * sin_k) return static_cast(fputil::multiply_add( - sin_k, cos_y, fputil::multiply_add(sin_y, cos_k, 0.0f))); + sin_y, cos_k, fputil::multiply_add(cosm1_y, sin_k, sin_k))); } } // namespace LIBC_NAMESPACE_DECL From 348eea1de6a56e3b92a102e9cb77d6783ea31141 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Wed, 9 Oct 2024 18:02:51 +0100 Subject: [PATCH 09/20] added mpfr exhaustive test --- libc/test/src/math/CMakeLists.txt | 11 ++++++++ libc/test/src/math/sinpif16_test.cpp | 41 ++++++++++++++++++++++++++++ libc/utils/MPFRWrapper/MPFRUtils.cpp | 6 ++++ 3 files changed, 58 insertions(+) create mode 100644 libc/test/src/math/sinpif16_test.cpp diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt index 0c4118c369454..8464c909721ec 100644 --- a/libc/test/src/math/CMakeLists.txt +++ b/libc/test/src/math/CMakeLists.txt @@ -90,6 +90,17 @@ add_fp_unittest( libc.src.__support.FPUtil.fp_bits ) +add_fp_unittest( + sinpif16_test + NEED_MPFR + SUITE + libc-math-unittests + SRCS + sinpif16_test.cpp + DEPENDS + libc.src.math.sinpif16 +) + add_fp_unittest( sin_test NEED_MPFR diff --git a/libc/test/src/math/sinpif16_test.cpp b/libc/test/src/math/sinpif16_test.cpp new file mode 100644 index 0000000000000..7821c6b223554 --- /dev/null +++ b/libc/test/src/math/sinpif16_test.cpp @@ -0,0 +1,41 @@ +//===-- Exhaustive test 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/math/sinpif16.h" +#include "test/UnitTest/FPMatcher.h" +#include "test/UnitTest/Test.h" +#include "utils/MPFRWrapper/MPFRUtils.h" + +using LlvmLibcSinpif16Test = 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(LlvmLibcSinpif16Test, PositiveRange) { + for (uint16_t v = POS_START; v <= POS_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::sinpif16, x, + LIBC_NAMESPACE::sinpif16(x), 0.5); + } +} + +TEST_F(LLvmLibcSinpif16Test, NegativeRange) { + for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) { + float16 x = FPBits(v).get_val(); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::sinpif16, x, + LIBC_NAMESPACE::sinpif16(x), 0.5); + } +} diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 27ff1f7190ef9..b1f44ba5eec71 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -498,6 +498,12 @@ class MPFRNumber { return result; } + MPFRNumber sinpif16() const { + MPFRNumber result(*this); + mpfr_sinpif16(result.value, value, mpfr_rounding); + return result; + } + MPFRNumber sinh() const { MPFRNumber result(*this); mpfr_sinh(result.value, value, mpfr_rounding); From 0ad7d39f2c0e6d18bae5c20937031235a0bde803 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Fri, 11 Oct 2024 14:31:14 +0100 Subject: [PATCH 10/20] minor fixes --- libc/src/math/generic/sinpif16.cpp | 19 +++++++++++-------- libc/test/src/math/sinpif16_test.cpp | 12 ++++++------ libc/utils/MPFRWrapper/MPFRUtils.cpp | 6 ------ 3 files changed, 17 insertions(+), 20 deletions(-) diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 0c96e6e173105..0b15432437b9d 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -20,8 +20,7 @@ namespace LIBC_NAMESPACE_DECL { // Lookup table for sin(k * pi / 32) with k = 0, ..., 63. // Table is generated with Sollya as follows: // > display = hexadecimmal; -// > prec = 24; -// > for k from 0 to 63 do {sin(k * pi/32);}; +// > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); }; const float SIN_K_PI_OVER_32[64] = {0, 0x1.917a6cp-4, 0x1.8f8b84p-3, @@ -166,16 +165,20 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { 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: - // > P = fpminimax(cos(x*pi/32), [|0, 2, 4, 6|], [|1, SG...|], [0, 0.5]); // 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) + // 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); + if (LIBC_UNLIKELY(sin_y == 0 && sin_k == 0)) { + return FPBits::zero(xbits.sign()).get_val(); + } + // Since, cosm1_y = cos_y - 1, therefore: // sin(x * pi) = cos_k * sin_y + sin_k + (cosm1_y * sin_k) return static_cast(fputil::multiply_add( diff --git a/libc/test/src/math/sinpif16_test.cpp b/libc/test/src/math/sinpif16_test.cpp index 7821c6b223554..fac6d4447ee8f 100644 --- a/libc/test/src/math/sinpif16_test.cpp +++ b/libc/test/src/math/sinpif16_test.cpp @@ -14,7 +14,7 @@ using LlvmLibcSinpif16Test = LIBC_NAMESPACE::testing::FPTest; -namespace mpfr = LIBC_NAMESPACE : testing::mpfr; +namespace mpfr = LIBC_NAMESPACE::testing::mpfr; // Range: [0, Inf] static constexpr uint16_t POS_START = 0x0000U; @@ -27,15 +27,15 @@ static constexpr uint16_t NEG_STOP = 0xfc00U; TEST_F(LlvmLibcSinpif16Test, PositiveRange) { for (uint16_t v = POS_START; v <= POS_STOP; ++v) { float16 x = FPBits(v).get_val(); - EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::sinpif16, x, - LIBC_NAMESPACE::sinpif16(x), 0.5); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinpi, x, + LIBC_NAMESPACE::sinpif16(x), 1); } } -TEST_F(LLvmLibcSinpif16Test, NegativeRange) { +TEST_F(LlvmLibcSinpif16Test, NegativeRange) { for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) { float16 x = FPBits(v).get_val(); - EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::sinpif16, x, - LIBC_NAMESPACE::sinpif16(x), 0.5); + EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinpi, x, + LIBC_NAMESPACE::sinpif16(x), 1); } } diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index b1f44ba5eec71..27ff1f7190ef9 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -498,12 +498,6 @@ class MPFRNumber { return result; } - MPFRNumber sinpif16() const { - MPFRNumber result(*this); - mpfr_sinpif16(result.value, value, mpfr_rounding); - return result; - } - MPFRNumber sinh() const { MPFRNumber result(*this); mpfr_sinh(result.value, value, mpfr_rounding); From 0f4af0776336de1a0b7389171b09720770b55cd0 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Fri, 11 Oct 2024 19:57:33 +0100 Subject: [PATCH 11/20] changed 0.0 to 0x0.0p0, fixed MPFRUtils computation of mpfr_sinpi(x) for mpfr versions < 4.2 --- libc/src/math/generic/sinpif16.cpp | 85 +++++++--------------------- libc/utils/MPFRWrapper/MPFRUtils.cpp | 14 ++++- 2 files changed, 32 insertions(+), 67 deletions(-) diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 0b15432437b9d..5138e0887cfc8 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -21,72 +21,25 @@ namespace LIBC_NAMESPACE_DECL { // Table is generated with Sollya as follows: // > display = hexadecimmal; // > for k from 0 to 63 do { round(sin(k * pi/32), SG, RN); }; -const float SIN_K_PI_OVER_32[64] = {0, - 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, - 0, - -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}; - -int32_t range_reduction(float x, float &y) { +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); diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 27ff1f7190ef9..22d5713c45d97 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -492,7 +492,19 @@ class MPFRNumber { MPFRNumber value_pi(0.0, 1280); mpfr_const_pi(value_pi.value, MPFR_RNDN); mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN); - mpfr_sin(result.value, value_pi.value, mpfr_rounding); + + int ret = mpfr_integer_p(value); + MPFRNumber value_mul_two(*this); + mpfr_mul_si(value_mul_two.value, value, 2, MPFR_RNDN); + if (mpfr_integer_p(value_mul_two.value) != 0) { + if (ret != 0) { + mpfr_set_si(result.value, 0, mpfr_rounding); + } else { + mpfr_sin(result.value, value_pi.value, MPFR_RNDN); + } + } else { + mpfr_sin(result.value, value_pi.value, mpfr_rounding); + } #endif return result; From 24f1173d4d34e16ac26406a92e5dbb8906818413 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Fri, 11 Oct 2024 21:07:15 +0100 Subject: [PATCH 12/20] update mpfr_sinpi implementation --- libc/utils/MPFRWrapper/MPFRUtils.cpp | 26 ++++++++++++++------------ 1 file changed, 14 insertions(+), 12 deletions(-) diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 22d5713c45d97..cb2bd6803c2d7 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -489,25 +489,27 @@ class MPFRNumber { mpfr_sinpi(result.value, value, mpfr_rounding); #else - MPFRNumber value_pi(0.0, 1280); - mpfr_const_pi(value_pi.value, MPFR_RNDN); - mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN); - int ret = mpfr_integer_p(value); + if (ret != 0) { + mpfr_set_si(result.value, 0, mpfr_rounding); + return result; + } + MPFRNumber value_mul_two(*this); mpfr_mul_si(value_mul_two.value, value, 2, MPFR_RNDN); + if (mpfr_integer_p(value_mul_two.value) != 0) { - if (ret != 0) { - mpfr_set_si(result.value, 0, mpfr_rounding); - } else { - mpfr_sin(result.value, value_pi.value, MPFR_RNDN); - } - } else { - mpfr_sin(result.value, value_pi.value, mpfr_rounding); + auto d = mpfr_get_si(value, MPFR_RNDD); + mpfr_set_si(result.value, (d & 1) ? -1 : 1, mpfr_rounding); + return result; } -#endif + MPFRNumber value_pi(0.0, 1280); + mpfr_const_pi(value_pi.value, MPFR_RNDN); + mpfr_mul(value_pi.value, value_pi.value, value, MPFR_RNDN); + mpfr_sin(result.value, value_pi.value, mpfr_rounding); return result; +#endif } MPFRNumber sinh() const { From cacbf19509556d7db89ee0f0b0fc572027a02328 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Fri, 11 Oct 2024 21:11:08 +0100 Subject: [PATCH 13/20] minor fix --- libc/utils/MPFRWrapper/MPFRUtils.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index cb2bd6803c2d7..77727cd11bcca 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -489,8 +489,7 @@ class MPFRNumber { mpfr_sinpi(result.value, value, mpfr_rounding); #else - int ret = mpfr_integer_p(value); - if (ret != 0) { + if (mpfr_integer_p(value)) { mpfr_set_si(result.value, 0, mpfr_rounding); return result; } From a543283d4f696028ff963345cf7f8c633cc3a798 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Fri, 11 Oct 2024 22:52:34 +0100 Subject: [PATCH 14/20] minor fixes, documentation update --- libc/config/linux/x86_64/entrypoints.txt | 4 +-- libc/docs/math/index.rst | 2 +- libc/src/math/generic/sinpif16.cpp | 33 +++++++++++------------- libc/test/src/math/sinpif16_test.cpp | 3 +-- libc/utils/MPFRWrapper/MPFRUtils.cpp | 1 + 5 files changed, 20 insertions(+), 23 deletions(-) diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index b7d0406fe6f00..c92c0fe75e139 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -573,7 +573,6 @@ set(TARGET_LIBM_ENTRYPOINTS libc.src.math.sinf libc.src.math.sinhf libc.src.math.sinpif - libc.src.math.sinpif16 libc.src.math.sqrt libc.src.math.sqrtf libc.src.math.sqrtl @@ -668,7 +667,8 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.scalblnf16 libc.src.math.scalbnf16 libc.src.math.setpayloadf16 - libc.src.math.setpayloadsigf16 + libc.src.math.setpayloadsigf16 + libc.src.math.sinpif16 libc.src.math.totalorderf16 libc.src.math.totalordermagf16 libc.src.math.truncf16 diff --git a/libc/docs/math/index.rst b/libc/docs/math/index.rst index 0e757e7978bfa..997cce8cef8a2 100644 --- a/libc/docs/math/index.rst +++ b/libc/docs/math/index.rst @@ -342,7 +342,7 @@ Higher Math Functions +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | sinh | |check| | | | | | 7.12.5.5 | F.10.2.5 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ -| sinpi | |check| | | | | | 7.12.4.13 | F.10.1.13 | +| sinpi | |check| | | | |check| | | 7.12.4.13 | F.10.1.13 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ | sqrt | |check| | |check| | |check| | | |check| | 7.12.7.10 | F.10.4.10 | +-----------+------------------+-----------------+------------------------+----------------------+------------------------+------------------------+----------------------------+ diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 5138e0887cfc8..4ab18986ee8b8 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -74,10 +74,9 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // * pi/32) and cos(y * pi/32) are computed using degree-9 chebyshev // polynomials generated by Sollya. - if (LIBC_UNLIKELY(x_abs == 0U)) { - // For signed zeros + // For signed zeros + if (LIBC_UNLIKELY(x_abs == 0U)) return x; - } // Numbers greater or equal to 2^10 are integers, or infinity, or NaN if (LIBC_UNLIKELY(x_abs >= 0x6400)) { @@ -101,8 +100,6 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { float sin_k = SIN_K_PI_OVER_32[k & 63]; float cos_k = SIN_K_PI_OVER_32[(k + 16) & 63]; - float cosm1_y, sin_y; - // Recall; // sin(x * pi/32) = sin((k + y) * pi/32) // = sin(y * pi/32) * cos(k * pi/32) + cos(y * pi/32) * sin(k * @@ -115,22 +112,22 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // 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); + 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]); - cosm1_y = ysq * fputil::polyeval(ysq, -0x1.3bd3ccp-8f, 0x1.03a61ap-18f, - 0x1.a6f7a2p-29f); - - if (LIBC_UNLIKELY(sin_y == 0 && sin_k == 0)) { + // 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); + + if (LIBC_UNLIKELY(sin_y == 0 && sin_k == 0)) return FPBits::zero(xbits.sign()).get_val(); - } // Since, cosm1_y = cos_y - 1, therefore: // sin(x * pi) = cos_k * sin_y + sin_k + (cosm1_y * sin_k) diff --git a/libc/test/src/math/sinpif16_test.cpp b/libc/test/src/math/sinpif16_test.cpp index fac6d4447ee8f..11b6c2b3c074d 100644 --- a/libc/test/src/math/sinpif16_test.cpp +++ b/libc/test/src/math/sinpif16_test.cpp @@ -1,5 +1,4 @@ -//===-- Exhaustive test for sinpif16 -//---------------------------------------===// +//===-- Exhaustive test 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. diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index 77727cd11bcca..fd4817e73ef7c 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -488,6 +488,7 @@ class MPFRNumber { (MPFR_VERSION_MAJOR == 4 && MPFR_VERSION_MINOR >= 2) mpfr_sinpi(result.value, value, mpfr_rounding); + return result; #else if (mpfr_integer_p(value)) { mpfr_set_si(result.value, 0, mpfr_rounding); From 743625d3eb1ad77bfe8a36bae63321ec27b8436f Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 14 Oct 2024 00:25:46 +0100 Subject: [PATCH 15/20] update pr, test branch against bot --- libc/config/linux/x86_64/entrypoints.txt | 2 +- libc/src/math/generic/sinpif16.cpp | 6 ++---- libc/test/src/math/sinpif16_test.cpp | 7 ++++--- libc/test/src/math/smoke/sinpif16_test.cpp | 17 ++++++++--------- libc/utils/MPFRWrapper/MPFRUtils.cpp | 2 +- 5 files changed, 16 insertions(+), 18 deletions(-) diff --git a/libc/config/linux/x86_64/entrypoints.txt b/libc/config/linux/x86_64/entrypoints.txt index c92c0fe75e139..d6417bfb6c3ef 100644 --- a/libc/config/linux/x86_64/entrypoints.txt +++ b/libc/config/linux/x86_64/entrypoints.txt @@ -667,7 +667,7 @@ if(LIBC_TYPES_HAS_FLOAT16) libc.src.math.scalblnf16 libc.src.math.scalbnf16 libc.src.math.setpayloadf16 - libc.src.math.setpayloadsigf16 + libc.src.math.setpayloadsigf16 libc.src.math.sinpif16 libc.src.math.totalorderf16 libc.src.math.totalordermagf16 diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 4ab18986ee8b8..054d44a3f11f0 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -110,8 +110,7 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { 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]); float sin_y = y * fputil::polyeval(ysq, 0x1.921fb6p-4f, -0x1.4aeabcp-13f, 0x1.a03354p-21f, -0x1.ad02d2p-20f); @@ -120,8 +119,7 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // 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: + // 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); diff --git a/libc/test/src/math/sinpif16_test.cpp b/libc/test/src/math/sinpif16_test.cpp index 11b6c2b3c074d..7156acd6a6785 100644 --- a/libc/test/src/math/sinpif16_test.cpp +++ b/libc/test/src/math/sinpif16_test.cpp @@ -1,4 +1,5 @@ -//===-- Exhaustive test for sinpif16---------------------------------------===// +//===-- Exhaustive test 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. @@ -27,7 +28,7 @@ TEST_F(LlvmLibcSinpif16Test, PositiveRange) { for (uint16_t v = POS_START; v <= POS_STOP; ++v) { float16 x = FPBits(v).get_val(); EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinpi, x, - LIBC_NAMESPACE::sinpif16(x), 1); + LIBC_NAMESPACE::sinpif16(x), 0.5); } } @@ -35,6 +36,6 @@ TEST_F(LlvmLibcSinpif16Test, NegativeRange) { for (uint16_t v = NEG_START; v <= NEG_STOP; ++v) { float16 x = FPBits(v).get_val(); EXPECT_MPFR_MATCH_ALL_ROUNDING(mpfr::Operation::Sinpi, x, - LIBC_NAMESPACE::sinpif16(x), 1); + LIBC_NAMESPACE::sinpif16(x), 0.5); } } diff --git a/libc/test/src/math/smoke/sinpif16_test.cpp b/libc/test/src/math/smoke/sinpif16_test.cpp index 0cf6f7c6a86e2..82623e3b3b14a 100644 --- a/libc/test/src/math/smoke/sinpif16_test.cpp +++ b/libc/test/src/math/smoke/sinpif16_test.cpp @@ -4,7 +4,6 @@ // See https://llvm.org/LICENSE.txt for license information. // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception // -// // ===----------------------------------------------------------------------==// #include "src/errno/libc_errno.h" @@ -20,10 +19,10 @@ TEST_F(LlvmLibcSinpif16Test, SpecialNumbers) { EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif16(aNaN)); EXPECT_MATH_ERRNO(0); - EXPECT_FP_EQ(0.0f, LIBC_NAMESPACE::sinpif16(0.0f)); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpif16(0.0f)); EXPECT_MATH_ERRNO(0); - EXPECT_FP_EQ(-0.0f, LIBC_NAMESPACE::sinpif16(-0.0f)); + EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpif16(-0.0f)); EXPECT_MATH_ERRNO(0); EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif16(inf)); @@ -34,10 +33,10 @@ TEST_F(LlvmLibcSinpif16Test, SpecialNumbers) { } TEST_F(LlvmLibcSinpif16Test, Integers) { - EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x420)); - EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x1p+10)); - EXPECT_FP_EQ(-0.0, LIBC_NAMESPACE::sinpif16(-0x1.4p+14)); - EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x420)); - EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x1.cp+15)); - EXPECT_FP_EQ(0.0, LIBC_NAMESPACE::sinpif16(0x1.cp+7)); + EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpif16(-0x420)); + EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpif16(-0x1p+10)); + EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpif16(-0x1.4p+14)); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpif16(0x420)); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpif16(0x1.cp+15)); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpif16(0x1.cp+7)); } diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp index fd4817e73ef7c..eecffc782c1a7 100644 --- a/libc/utils/MPFRWrapper/MPFRUtils.cpp +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -498,7 +498,7 @@ class MPFRNumber { MPFRNumber value_mul_two(*this); mpfr_mul_si(value_mul_two.value, value, 2, MPFR_RNDN); - if (mpfr_integer_p(value_mul_two.value) != 0) { + if (mpfr_integer_p(value_mul_two.value)) { auto d = mpfr_get_si(value, MPFR_RNDD); mpfr_set_si(result.value, (d & 1) ? -1 : 1, mpfr_rounding); return result; From 9f85a852f3afb3e58c4ec37c4d4b9e865127a6ce Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 14 Oct 2024 10:54:54 +0100 Subject: [PATCH 16/20] formatting --- libc/test/src/math/sinpif16_test.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/libc/test/src/math/sinpif16_test.cpp b/libc/test/src/math/sinpif16_test.cpp index 7156acd6a6785..f47d671f1b507 100644 --- a/libc/test/src/math/sinpif16_test.cpp +++ b/libc/test/src/math/sinpif16_test.cpp @@ -1,5 +1,4 @@ -//===-- Exhaustive test for sinpif16 -//---------------------------------------===// +//===-- Exhaustive test 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. From ffcbd23dbcb02ccccc0f5d0ffc9baf9eb2a58cde Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 14 Oct 2024 11:17:09 +0100 Subject: [PATCH 17/20] clang-format fix --- libc/test/src/math/sinpif16_test.cpp | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/libc/test/src/math/sinpif16_test.cpp b/libc/test/src/math/sinpif16_test.cpp index f47d671f1b507..7156acd6a6785 100644 --- a/libc/test/src/math/sinpif16_test.cpp +++ b/libc/test/src/math/sinpif16_test.cpp @@ -1,4 +1,5 @@ -//===-- Exhaustive test for sinpif16 ---------------------------------------===// +//===-- Exhaustive test 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. From d6fad4c8e1f76b7d63e68f6194bce8fa6497ac90 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 14 Oct 2024 14:59:37 +0100 Subject: [PATCH 18/20] clang-format --- libc/test/src/math/sinpif16_test.cpp | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/libc/test/src/math/sinpif16_test.cpp b/libc/test/src/math/sinpif16_test.cpp index 7156acd6a6785..cd1710ce3b445 100644 --- a/libc/test/src/math/sinpif16_test.cpp +++ b/libc/test/src/math/sinpif16_test.cpp @@ -1,11 +1,10 @@ -//===-- Exhaustive test for sinpif16 -//---------------------------------------===// +//===-- Exhaustive test 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/math/sinpif16.h" #include "test/UnitTest/FPMatcher.h" From 06c9b4b54f64cf67c998ecd1b657d602fc9bf0c0 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 14 Oct 2024 22:08:53 +0100 Subject: [PATCH 19/20] fixed float16 cast, nit formatting --- libc/src/math/generic/CMakeLists.txt | 1 + libc/src/math/generic/sinpif16.cpp | 3 ++- libc/test/src/math/sinpif16_test.cpp | 2 +- libc/test/src/math/smoke/sinpif16_test.cpp | 6 +++--- 4 files changed, 7 insertions(+), 5 deletions(-) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index 48078a158bfb9..d0611ee4253ee 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -541,6 +541,7 @@ add_entrypoint_object( libc.src.__support.FPUtil.multiply_add libc.src.__support.FPUtil.nearest_integer libc.src.__support.FPUtil.polyeval + libc.src.__support.FPUtil.cast libc.src.__support.macros.properties.types COMPILE_OPTIONS -O3 diff --git a/libc/src/math/generic/sinpif16.cpp b/libc/src/math/generic/sinpif16.cpp index 054d44a3f11f0..17cca583e0c0e 100644 --- a/libc/src/math/generic/sinpif16.cpp +++ b/libc/src/math/generic/sinpif16.cpp @@ -10,6 +10,7 @@ #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" @@ -129,7 +130,7 @@ LLVM_LIBC_FUNCTION(float16, sinpif16, (float16 x)) { // Since, cosm1_y = cos_y - 1, therefore: // sin(x * pi) = cos_k * sin_y + sin_k + (cosm1_y * sin_k) - return static_cast(fputil::multiply_add( + 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/test/src/math/sinpif16_test.cpp b/libc/test/src/math/sinpif16_test.cpp index cd1710ce3b445..8477124b2e6e3 100644 --- a/libc/test/src/math/sinpif16_test.cpp +++ b/libc/test/src/math/sinpif16_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/sinpif16.h" #include "test/UnitTest/FPMatcher.h" diff --git a/libc/test/src/math/smoke/sinpif16_test.cpp b/libc/test/src/math/smoke/sinpif16_test.cpp index 82623e3b3b14a..0bcd38a60d849 100644 --- a/libc/test/src/math/smoke/sinpif16_test.cpp +++ b/libc/test/src/math/smoke/sinpif16_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/errno/libc_errno.h" #include "src/math/sinpif16.h" @@ -19,10 +19,10 @@ TEST_F(LlvmLibcSinpif16Test, SpecialNumbers) { EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif16(aNaN)); EXPECT_MATH_ERRNO(0); - EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpif16(0.0f)); + EXPECT_FP_EQ(zero, LIBC_NAMESPACE::sinpif16(zero)); EXPECT_MATH_ERRNO(0); - EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpif16(-0.0f)); + EXPECT_FP_EQ(neg_zero, LIBC_NAMESPACE::sinpif16(neg_zero)); EXPECT_MATH_ERRNO(0); EXPECT_FP_EQ(aNaN, LIBC_NAMESPACE::sinpif16(inf)); From 87faa10933723b901df15f143c617dd51d893e62 Mon Sep 17 00:00:00 2001 From: wldfngrs Date: Mon, 14 Oct 2024 22:39:57 +0100 Subject: [PATCH 20/20] sort dependencies --- libc/src/math/generic/CMakeLists.txt | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/libc/src/math/generic/CMakeLists.txt b/libc/src/math/generic/CMakeLists.txt index d0611ee4253ee..4915b23bb383e 100644 --- a/libc/src/math/generic/CMakeLists.txt +++ b/libc/src/math/generic/CMakeLists.txt @@ -536,12 +536,12 @@ add_entrypoint_object( ../sinpif16.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.FPUtil.cast libc.src.__support.macros.properties.types COMPILE_OPTIONS -O3