diff --git a/libc/AOR_v20.02/math/cosf.c b/libc/AOR_v20.02/math/cosf.c deleted file mode 100644 --- a/libc/AOR_v20.02/math/cosf.c +++ /dev/null @@ -1,64 +0,0 @@ -/* - * Single-precision cos 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 -#include -#include "math_config.h" -#include "sincosf.h" - -/* Fast cosf implementation. Worst-case ULP is 0.5607, maximum relative - error is 0.5303 * 2^-23. A single-step range reduction is used for - small values. Large inputs have their range reduced using fast integer - arithmetic. */ -float -cosf (float y) -{ - double x = y; - double s; - int n; - const sincos_t *p = &__sincosf_table[0]; - - if (abstop12 (y) < abstop12 (pio4)) - { - double x2 = x * x; - - if (unlikely (abstop12 (y) < abstop12 (0x1p-12f))) - return 1.0f; - - return sinf_poly (x, x2, p, 1); - } - else if (likely (abstop12 (y) < abstop12 (120.0f))) - { - x = reduce_fast (x, p, &n); - - /* Setup the signs for sin and cos. */ - s = p->sign[n & 3]; - - if (n & 2) - p = &__sincosf_table[1]; - - return sinf_poly (x * s, x * x, p, n ^ 1); - } - else if (abstop12 (y) < abstop12 (INFINITY)) - { - uint32_t xi = asuint (y); - int sign = xi >> 31; - - x = reduce_large (xi, &n); - - /* Setup signs for sin and cos - include original sign. */ - s = p->sign[(n + sign) & 3]; - - if ((n + sign) & 2) - p = &__sincosf_table[1]; - - return sinf_poly (x * s, x * x, p, n ^ 1); - } - else - return __math_invalidf (y); -} diff --git a/libc/AOR_v20.02/math/sincosf.h b/libc/AOR_v20.02/math/sincosf.h deleted file mode 100644 --- a/libc/AOR_v20.02/math/sincosf.h +++ /dev/null @@ -1,154 +0,0 @@ -/* - * Header for sinf, cosf and sincosf. - * - * 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 -#include -#include "math_config.h" - -/* 2PI * 2^-64. */ -static const double pi63 = 0x1.921FB54442D18p-62; -/* PI / 4. */ -static const double pio4 = 0x1.921FB54442D18p-1; - -/* The constants and polynomials for sine and cosine. */ -typedef struct -{ - double sign[4]; /* Sign of sine in quadrants 0..3. */ - double hpi_inv; /* 2 / PI ( * 2^24 if !TOINT_INTRINSICS). */ - double hpi; /* PI / 2. */ - double c0, c1, c2, c3, c4; /* Cosine polynomial. */ - double s1, s2, s3; /* Sine polynomial. */ -} sincos_t; - -/* Polynomial data (the cosine polynomial is negated in the 2nd entry). */ -extern const sincos_t __sincosf_table[2] HIDDEN; - -/* Table with 4/PI to 192 bit precision. */ -extern const uint32_t __inv_pio4[] HIDDEN; - -/* Top 12 bits of the float representation with the sign bit cleared. */ -static inline uint32_t -abstop12 (float x) -{ - return (asuint (x) >> 20) & 0x7ff; -} - -/* Compute the sine and cosine of inputs X and X2 (X squared), using the - polynomial P and store the results in SINP and COSP. N is the quadrant, - if odd the cosine and sine polynomials are swapped. */ -static inline void -sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp, - float *cosp) -{ - double x3, x4, x5, x6, s, c, c1, c2, s1; - - x4 = x2 * x2; - x3 = x2 * x; - c2 = p->c3 + x2 * p->c4; - s1 = p->s2 + x2 * p->s3; - - /* Swap sin/cos result based on quadrant. */ - float *tmp = (n & 1 ? cosp : sinp); - cosp = (n & 1 ? sinp : cosp); - sinp = tmp; - - c1 = p->c0 + x2 * p->c1; - x5 = x3 * x2; - x6 = x4 * x2; - - s = x + x3 * p->s1; - c = c1 + x4 * p->c2; - - *sinp = s + x5 * s1; - *cosp = c + x6 * c2; -} - -/* Return the sine of inputs X and X2 (X squared) using the polynomial P. - N is the quadrant, and if odd the cosine polynomial is used. */ -static inline float -sinf_poly (double x, double x2, const sincos_t *p, int n) -{ - double x3, x4, x6, x7, s, c, c1, c2, s1; - - if ((n & 1) == 0) - { - x3 = x * x2; - s1 = p->s2 + x2 * p->s3; - - x7 = x3 * x2; - s = x + x3 * p->s1; - - return s + x7 * s1; - } - else - { - x4 = x2 * x2; - c2 = p->c3 + x2 * p->c4; - c1 = p->c0 + x2 * p->c1; - - x6 = x4 * x2; - c = c1 + x4 * p->c2; - - return c + x6 * c2; - } -} - -/* Fast range reduction using single multiply-subtract. Return the modulo of - X as a value between -PI/4 and PI/4 and store the quadrant in NP. - The values for PI/2 and 2/PI are accessed via P. Since PI/2 as a double - is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4, - the result is accurate for |X| <= 120.0. */ -static inline double -reduce_fast (double x, const sincos_t *p, int *np) -{ - double r; -#if TOINT_INTRINSICS - /* Use fast round and lround instructions when available. */ - r = x * p->hpi_inv; - *np = converttoint (r); - return x - roundtoint (r) * p->hpi; -#else - /* Use scaled float to int conversion with explicit rounding. - hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31. - This avoids inaccuracies introduced by truncating negative values. */ - r = x * p->hpi_inv; - int n = ((int32_t)r + 0x800000) >> 24; - *np = n; - return x - n * p->hpi; -#endif -} - -/* Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic. - XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored). - Return the modulo between -PI/4 and PI/4 and store the quadrant in NP. - Reduction uses a table of 4/PI with 192 bits of precision. A 32x96->128 bit - multiply computes the exact 2.62-bit fixed-point modulo. Since the result - can have at most 29 leading zeros after the binary point, the double - precision result is accurate to 33 bits. */ -static inline double -reduce_large (uint32_t xi, int *np) -{ - const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15]; - int shift = (xi >> 23) & 7; - uint64_t n, res0, res1, res2; - - xi = (xi & 0xffffff) | 0x800000; - xi <<= shift; - - res0 = xi * arr[0]; - res1 = (uint64_t)xi * arr[4]; - res2 = (uint64_t)xi * arr[8]; - res0 = (res2 >> 32) | (res0 << 32); - res0 += res1; - - n = (res0 + (1ULL << 61)) >> 62; - res0 -= n << 62; - double x = (int64_t)res0; - *np = n; - return x * pi63; -} diff --git a/libc/AOR_v20.02/math/sincosf.c b/libc/AOR_v20.02/math/sincosf.c deleted file mode 100644 --- a/libc/AOR_v20.02/math/sincosf.c +++ /dev/null @@ -1,80 +0,0 @@ -/* - * Single-precision sin/cos 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 -#include -#include "math_config.h" -#include "sincosf.h" - -/* Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative - error is 0.5303 * 2^-23. A single-step range reduction is used for - small values. Large inputs have their range reduced using fast integer - arithmetic. */ -void -sincosf (float y, float *sinp, float *cosp) -{ - double x = y; - double s; - int n; - const sincos_t *p = &__sincosf_table[0]; - - if (abstop12 (y) < abstop12 (pio4)) - { - double x2 = x * x; - - if (unlikely (abstop12 (y) < abstop12 (0x1p-12f))) - { - if (unlikely (abstop12 (y) < abstop12 (0x1p-126f))) - /* Force underflow for tiny y. */ - force_eval_float (x2); - *sinp = y; - *cosp = 1.0f; - return; - } - - sincosf_poly (x, x2, p, 0, sinp, cosp); - } - else if (abstop12 (y) < abstop12 (120.0f)) - { - x = reduce_fast (x, p, &n); - - /* Setup the signs for sin and cos. */ - s = p->sign[n & 3]; - - if (n & 2) - p = &__sincosf_table[1]; - - sincosf_poly (x * s, x * x, p, n, sinp, cosp); - } - else if (likely (abstop12 (y) < abstop12 (INFINITY))) - { - uint32_t xi = asuint (y); - int sign = xi >> 31; - - x = reduce_large (xi, &n); - - /* Setup signs for sin and cos - include original sign. */ - s = p->sign[(n + sign) & 3]; - - if ((n + sign) & 2) - p = &__sincosf_table[1]; - - sincosf_poly (x * s, x * x, p, n, sinp, cosp); - } - else - { - /* Return NaN if Inf or NaN for both sin and cos. */ - *sinp = *cosp = y - y; -#if WANT_ERRNO - /* Needed to set errno for +-Inf, the add is a hack to work - around a gcc register allocation issue: just passing y - affects code generation in the fast path. */ - __math_invalidf (y + y); -#endif - } -} diff --git a/libc/AOR_v20.02/math/sincosf_data.c b/libc/AOR_v20.02/math/sincosf_data.c deleted file mode 100644 --- a/libc/AOR_v20.02/math/sincosf_data.c +++ /dev/null @@ -1,64 +0,0 @@ -/* - * Data definition for sinf, cosf and sincosf. - * - * 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 -#include -#include "math_config.h" -#include "sincosf.h" - -/* The constants and polynomials for sine and cosine. The 2nd entry - computes -cos (x) rather than cos (x) to get negation for free. */ -const sincos_t __sincosf_table[2] = -{ - { - { 1.0, -1.0, -1.0, 1.0 }, -#if TOINT_INTRINSICS - 0x1.45F306DC9C883p-1, -#else - 0x1.45F306DC9C883p+23, -#endif - 0x1.921FB54442D18p0, - 0x1p0, - -0x1.ffffffd0c621cp-2, - 0x1.55553e1068f19p-5, - -0x1.6c087e89a359dp-10, - 0x1.99343027bf8c3p-16, - -0x1.555545995a603p-3, - 0x1.1107605230bc4p-7, - -0x1.994eb3774cf24p-13 - }, - { - { 1.0, -1.0, -1.0, 1.0 }, -#if TOINT_INTRINSICS - 0x1.45F306DC9C883p-1, -#else - 0x1.45F306DC9C883p+23, -#endif - 0x1.921FB54442D18p0, - -0x1p0, - 0x1.ffffffd0c621cp-2, - -0x1.55553e1068f19p-5, - 0x1.6c087e89a359dp-10, - -0x1.99343027bf8c3p-16, - -0x1.555545995a603p-3, - 0x1.1107605230bc4p-7, - -0x1.994eb3774cf24p-13 - } -}; - -/* Table with 4/PI to 192 bit precision. To avoid unaligned accesses - only 8 new bits are added per entry, making the table 4 times larger. */ -const uint32_t __inv_pio4[24] = -{ - 0xa2, 0xa2f9, 0xa2f983, 0xa2f9836e, - 0xf9836e4e, 0x836e4e44, 0x6e4e4415, 0x4e441529, - 0x441529fc, 0x1529fc27, 0x29fc2757, 0xfc2757d1, - 0x2757d1f5, 0x57d1f534, 0xd1f534dd, 0xf534ddc0, - 0x34ddc0db, 0xddc0db62, 0xc0db6295, 0xdb629599, - 0x6295993c, 0x95993c43, 0x993c4390, 0x3c439041 -}; diff --git a/libc/AOR_v20.02/math/sinf.c b/libc/AOR_v20.02/math/sinf.c deleted file mode 100644 --- a/libc/AOR_v20.02/math/sinf.c +++ /dev/null @@ -1,68 +0,0 @@ -/* - * Single-precision sin function. - * - * Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. - * See https://llvm.org/LICENSE.txt for license information. - * SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception - */ - -#include -#include "math_config.h" -#include "sincosf.h" - -/* Fast sinf implementation. Worst-case ULP is 0.5607, maximum relative - error is 0.5303 * 2^-23. A single-step range reduction is used for - small values. Large inputs have their range reduced using fast integer - arithmetic. */ -float -sinf (float y) -{ - double x = y; - double s; - int n; - const sincos_t *p = &__sincosf_table[0]; - - if (abstop12 (y) < abstop12 (pio4)) - { - s = x * x; - - if (unlikely (abstop12 (y) < abstop12 (0x1p-12f))) - { - if (unlikely (abstop12 (y) < abstop12 (0x1p-126f))) - /* Force underflow for tiny y. */ - force_eval_float (s); - return y; - } - - return sinf_poly (x, s, p, 0); - } - else if (likely (abstop12 (y) < abstop12 (120.0f))) - { - x = reduce_fast (x, p, &n); - - /* Setup the signs for sin and cos. */ - s = p->sign[n & 3]; - - if (n & 2) - p = &__sincosf_table[1]; - - return sinf_poly (x * s, x * x, p, n); - } - else if (abstop12 (y) < abstop12 (INFINITY)) - { - uint32_t xi = asuint (y); - int sign = xi >> 31; - - x = reduce_large (xi, &n); - - /* Setup signs for sin and cos - include original sign. */ - s = p->sign[(n + sign) & 3]; - - if ((n + sign) & 2) - p = &__sincosf_table[1]; - - return sinf_poly (x * s, x * x, p, n); - } - else - return __math_invalidf (y); -} diff --git a/libc/config/linux/api.td b/libc/config/linux/api.td --- a/libc/config/linux/api.td +++ b/libc/config/linux/api.td @@ -1,5 +1,6 @@ include "config/public_api.td" +include "spec/gnu_ext.td" include "spec/linux.td" include "spec/posix.td" include "spec/stdc.td" @@ -117,7 +118,10 @@ IsNanMacro, ]; let Functions = [ + "cosf", "round", + "sincosf", + "sinf", ]; } diff --git a/libc/lib/CMakeLists.txt b/libc/lib/CMakeLists.txt --- a/libc/lib/CMakeLists.txt +++ b/libc/lib/CMakeLists.txt @@ -41,7 +41,10 @@ llvmlibm DEPENDS # math.h entrypoints + cosf round + sincosf + sinf ) add_redirector_library( diff --git a/libc/src/__support/common.h.def b/libc/src/__support/common.h.def --- a/libc/src/__support/common.h.def +++ b/libc/src/__support/common.h.def @@ -11,6 +11,10 @@ #define LIBC_INLINE_ASM __asm__ __volatile__ +#define likely(x) __builtin_expect (!!(x), 1) +#define unlikely(x) __builtin_expect (x, 0) +#define UNUSED __attribute__((unused)) + Include the platform specific definitions at build time. For example, that of entrypoint macro. %%include_file(${platform_defs}) diff --git a/libc/src/math/CMakeLists.txt b/libc/src/math/CMakeLists.txt --- a/libc/src/math/CMakeLists.txt +++ b/libc/src/math/CMakeLists.txt @@ -1,3 +1,24 @@ +add_header_library( + math_utils + HDRS + math_utils.h + DEPENDS + errno_h + math_h + DEPENDS + __errno_location +) + +add_object_library( + sincosf_utils + HDRS + sincosf_utils.h + SRCS + sincosf_data.cpp + DEPENDS + math_utils +) + add_entrypoint_object( round REDIRECTED @@ -12,3 +33,39 @@ SRC round_redirector.cpp ) + +add_entrypoint_object( + cosf + SRCS + cosf.cpp + HDRS + cosf.h + DEPENDS + math_h + math_utils + sincosf_utils +) + +add_entrypoint_object( + sinf + SRCS + sinf.cpp + HDRS + sinf.h + DEPENDS + math_h + math_utils + sincosf_utils +) + +add_entrypoint_object( + sincosf + SRCS + sincosf.cpp + HDRS + sincosf.h + DEPENDS + math_h + math_utils + sincosf_utils +) diff --git a/libc/src/math/cosf.h b/libc/src/math/cosf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/cosf.h @@ -0,0 +1,18 @@ +//===------------------- Implementation header for cosf ------------------===// +// +// 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_COSF_H +#define LLVM_LIBC_SRC_MATH_COSF_H + +namespace __llvm_libc { + +float cosf(float x); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_COSF_H diff --git a/libc/src/math/cosf.cpp b/libc/src/math/cosf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/cosf.cpp @@ -0,0 +1,64 @@ +//===------------------- Single-precision cos 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 "math_utils.h" +#include "sincosf_utils.h" + +#include "include/math.h" +#include "src/__support/common.h" + +#include + +namespace __llvm_libc { + +// Fast cosf implementation. Worst-case ULP is 0.5607, maximum relative +// error is 0.5303 * 2^-23. A single-step range reduction is used for +// small values. Large inputs have their range reduced using fast integer +// arithmetic. +float LLVM_LIBC_ENTRYPOINT(cosf)(float y) { + double x = y; + double s; + int n; + const sincos_t *p = &__sincosf_table[0]; + + if (abstop12(y) < abstop12(pio4)) { + double x2 = x * x; + + if (unlikely(abstop12(y) < abstop12(0x1p-12f))) + return 1.0f; + + return sinf_poly(x, x2, p, 1); + } else if (likely(abstop12(y) < abstop12(120.0f))) { + x = reduce_fast(x, p, &n); + + // Setup the signs for sin and cos. + s = p->sign[n & 3]; + + if (n & 2) + p = &__sincosf_table[1]; + + return sinf_poly(x * s, x * x, p, n ^ 1); + } else if (abstop12(y) < abstop12(INFINITY)) { + uint32_t xi = asuint32_bits(y); + int sign = xi >> 31; + + x = reduce_large(xi, &n); + + // Setup signs for sin and cos - include original sign. + s = p->sign[(n + sign) & 3]; + + if ((n + sign) & 2) + p = &__sincosf_table[1]; + + return sinf_poly(x * s, x * x, p, n ^ 1); + } else { + return invalidf(y); + } +} + +} // namespace __llvm_libc diff --git a/libc/src/math/math_utils.h b/libc/src/math/math_utils.h new file mode 100644 --- /dev/null +++ b/libc/src/math/math_utils.h @@ -0,0 +1,35 @@ +#ifndef LLVM_LIBC_SRC_MATH_MATH_UTILS_H +#define LLVM_LIBC_SRC_MATH_MATH_UTILS_H + +#include "include/errno.h" +#include "include/math.h" + +#include "src/__support/common.h" +#include "src/errno/llvmlibc_errno.h" + +#include + +namespace __llvm_libc { + +float with_errnof(float x, int err) { + if (math_errhandling & MATH_ERRNO) + llvmlibc_errno = err; + return x; +} + +static inline uint32_t asuint32_bits(float x) { + return *reinterpret_cast(&x); +} + +static inline float invalidf(float x) { + float y = (x - x) / (x - x); + return isnan(x) ? y : with_errnof(y, EDOM); +} + +static inline void force_eval_float(float x) { + volatile float y UNUSED = x; +} + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_MATH_UTILS_H diff --git a/libc/src/math/sincosf.h b/libc/src/math/sincosf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/sincosf.h @@ -0,0 +1,18 @@ +//===------------------- Implementation header for sinf ------------------===// +// +// 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_SINF_H +#define LLVM_LIBC_SRC_MATH_SINF_H + +namespace __llvm_libc { + +void sincosf(float x, float *sinx, float *cosx); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_SINF_H diff --git a/libc/src/math/sincosf.cpp b/libc/src/math/sincosf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/sincosf.cpp @@ -0,0 +1,76 @@ +//===------------------ Single-precision sincos 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 "math_utils.h" +#include "sincosf_utils.h" + +#include +#include "src/__support/common.h" + +#include + +namespace __llvm_libc { + +// Fast sincosf implementation. Worst-case ULP is 0.5607, maximum relative +// error is 0.5303 * 2^-23. A single-step range reduction is used for +// small values. Large inputs have their range reduced using fast integer +// arithmetic. +void LLVM_LIBC_ENTRYPOINT(sincosf)(float y, float *sinp, float *cosp) { + double x = y; + double s; + int n; + const sincos_t *p = &__sincosf_table[0]; + + if (abstop12(y) < abstop12(pio4)) { + double x2 = x * x; + + if (unlikely(abstop12(y) < abstop12(0x1p-12f))) { + if (unlikely(abstop12(y) < abstop12(0x1p-126f))) + // Force underflow for tiny y. + force_eval_float(x2); + *sinp = y; + *cosp = 1.0f; + return; + } + + sincosf_poly(x, x2, p, 0, sinp, cosp); + } else if (abstop12(y) < abstop12(120.0f)) { + x = reduce_fast(x, p, &n); + + // Setup the signs for sin and cos. + s = p->sign[n & 3]; + + if (n & 2) + p = &__sincosf_table[1]; + + sincosf_poly(x * s, x * x, p, n, sinp, cosp); + } else if (likely(abstop12(y) < abstop12(INFINITY))) { + uint32_t xi = asuint32_bits(y); + int sign = xi >> 31; + + x = reduce_large(xi, &n); + + // Setup signs for sin and cos - include original sign. + s = p->sign[(n + sign) & 3]; + + if ((n + sign) & 2) + p = &__sincosf_table[1]; + + sincosf_poly(x * s, x * x, p, n, sinp, cosp); + } else { + //Return NaN if Inf or NaN for both sin and cos. + *sinp = *cosp = y - y; + + // Needed to set errno for +-Inf, the add is a hack to work + // around a gcc register allocation issue: just passing y + // affects code generation in the fast path. + invalidf(y + y); + } +} + +} // namespace __llvm_libc diff --git a/libc/src/math/sincosf_data.cpp b/libc/src/math/sincosf_data.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/sincosf_data.cpp @@ -0,0 +1,49 @@ +//===------------------------- sinf/cosf data tables ----------------------===// +// +// 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 "math_utils.h" +#include "sincosf_utils.h" + +#include + +namespace __llvm_libc { + +// The constants and polynomials for sine and cosine. The 2nd entry +// computes -cos (x) rather than cos (x) to get negation for free. +const sincos_t __sincosf_table[2] = {{{1.0, -1.0, -1.0, 1.0}, + 0x1.45F306DC9C883p+23, + 0x1.921FB54442D18p0, + 0x1p0, + -0x1.ffffffd0c621cp-2, + 0x1.55553e1068f19p-5, + -0x1.6c087e89a359dp-10, + 0x1.99343027bf8c3p-16, + -0x1.555545995a603p-3, + 0x1.1107605230bc4p-7, + -0x1.994eb3774cf24p-13}, + {{1.0, -1.0, -1.0, 1.0}, + 0x1.45F306DC9C883p+23, + 0x1.921FB54442D18p0, + -0x1p0, + 0x1.ffffffd0c621cp-2, + -0x1.55553e1068f19p-5, + 0x1.6c087e89a359dp-10, + -0x1.99343027bf8c3p-16, + -0x1.555545995a603p-3, + 0x1.1107605230bc4p-7, + -0x1.994eb3774cf24p-13}}; + +// Table with 4/PI to 192 bit precision. To avoid unaligned accesses +// only 8 new bits are added per entry, making the table 4 times larger. +const uint32_t __inv_pio4[24] = { + 0xa2, 0xa2f9, 0xa2f983, 0xa2f9836e, 0xf9836e4e, 0x836e4e44, + 0x6e4e4415, 0x4e441529, 0x441529fc, 0x1529fc27, 0x29fc2757, 0xfc2757d1, + 0x2757d1f5, 0x57d1f534, 0xd1f534dd, 0xf534ddc0, 0x34ddc0db, 0xddc0db62, + 0xc0db6295, 0xdb629599, 0x6295993c, 0x95993c43, 0x993c4390, 0x3c439041}; + +} // namespace __llvm_libc diff --git a/libc/src/math/sincosf_utils.h b/libc/src/math/sincosf_utils.h new file mode 100644 --- /dev/null +++ b/libc/src/math/sincosf_utils.h @@ -0,0 +1,142 @@ +//===------------- Collection of utils for cosf/sinf/sincosf --------------===// +// +// 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_SINCOSF_H +#define LLVM_LIBC_SRC_MATH_SINCOSF_H + +#include "math_utils.h" + +#include + +namespace __llvm_libc { + +// 2PI * 2^-64. +static const double pi63 = 0x1.921FB54442D18p-62; +// PI / 4. +static const double pio4 = 0x1.921FB54442D18p-1; + +// The constants and polynomials for sine and cosine. +typedef struct { + double sign[4]; // Sign of sine in quadrants 0..3. + double hpi_inv; // 2 / PI ( * 2^24 ). + double hpi; // PI / 2. + double c0, c1, c2, c3, c4; // Cosine polynomial. + double s1, s2, s3; // Sine polynomial. +} sincos_t; + +// Polynomial data (the cosine polynomial is negated in the 2nd entry). +extern const sincos_t __sincosf_table[2]; + +// Table with 4/PI to 192 bit precision. +extern const uint32_t __inv_pio4[]; + +// Top 12 bits of the float representation with the sign bit cleared. +static inline uint32_t abstop12(float x) { + return (asuint32_bits(x) >> 20) & 0x7ff; +} + +// Compute the sine and cosine of inputs X and X2 (X squared), using the +// polynomial P and store the results in SINP and COSP. N is the quadrant, +// if odd the cosine and sine polynomials are swapped. +static inline void sincosf_poly(double x, double x2, const sincos_t *p, int n, + float *sinp, float *cosp) { + double x3, x4, x5, x6, s, c, c1, c2, s1; + + x4 = x2 * x2; + x3 = x2 * x; + c2 = p->c3 + x2 * p->c4; + s1 = p->s2 + x2 * p->s3; + + // Swap sin/cos result based on quadrant. + float *tmp = (n & 1 ? cosp : sinp); + cosp = (n & 1 ? sinp : cosp); + sinp = tmp; + + c1 = p->c0 + x2 * p->c1; + x5 = x3 * x2; + x6 = x4 * x2; + + s = x + x3 * p->s1; + c = c1 + x4 * p->c2; + + *sinp = s + x5 * s1; + *cosp = c + x6 * c2; +} + +// Return the sine of inputs X and X2 (X squared) using the polynomial P. +// N is the quadrant, and if odd the cosine polynomial is used. +static inline float sinf_poly(double x, double x2, const sincos_t *p, int n) { + double x3, x4, x6, x7, s, c, c1, c2, s1; + + if ((n & 1) == 0) { + x3 = x * x2; + s1 = p->s2 + x2 * p->s3; + + x7 = x3 * x2; + s = x + x3 * p->s1; + + return s + x7 * s1; + } else { + x4 = x2 * x2; + c2 = p->c3 + x2 * p->c4; + c1 = p->c0 + x2 * p->c1; + + x6 = x4 * x2; + c = c1 + x4 * p->c2; + + return c + x6 * c2; + } +} + +// Fast range reduction using single multiply-subtract. Return the modulo of +// X as a value between -PI/4 and PI/4 and store the quadrant in NP. +// The values for PI/2 and 2/PI are accessed via P. Since PI/2 as a double +// is accurate to 55 bits and the worst-case cancellation happens at 6 * PI/4, +// the result is accurate for |X| <= 120.0. +static inline double reduce_fast(double x, const sincos_t *p, int *np) { + double r; + // Use scaled float to int conversion with explicit rounding. + // hpi_inv is prescaled by 2^24 so the quadrant ends up in bits 24..31. + // This avoids inaccuracies introduced by truncating negative values. + r = x * p->hpi_inv; + int n = ((int32_t)r + 0x800000) >> 24; + *np = n; + return x - n * p->hpi; +} + +// Reduce the range of XI to a multiple of PI/2 using fast integer arithmetic. +// XI is a reinterpreted float and must be >= 2.0f (the sign bit is ignored). +// Return the modulo between -PI/4 and PI/4 and store the quadrant in NP. +// Reduction uses a table of 4/PI with 192 bits of precision. A 32x96->128 bit +// multiply computes the exact 2.62-bit fixed-point modulo. Since the result +// can have at most 29 leading zeros after the binary point, the double +// precision result is accurate to 33 bits. +static inline double reduce_large(uint32_t xi, int *np) { + const uint32_t *arr = &__inv_pio4[(xi >> 26) & 15]; + int shift = (xi >> 23) & 7; + uint64_t n, res0, res1, res2; + + xi = (xi & 0xffffff) | 0x800000; + xi <<= shift; + + res0 = xi * arr[0]; + res1 = (uint64_t)xi * arr[4]; + res2 = (uint64_t)xi * arr[8]; + res0 = (res2 >> 32) | (res0 << 32); + res0 += res1; + + n = (res0 + (1ULL << 61)) >> 62; + res0 -= n << 62; + double x = (int64_t)res0; + *np = n; + return x * pi63; +} + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_SINCOSF_H diff --git a/libc/src/math/sinf.h b/libc/src/math/sinf.h new file mode 100644 --- /dev/null +++ b/libc/src/math/sinf.h @@ -0,0 +1,18 @@ +//===------------------- Implementation header for sinf ------------------===// +// +// 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_SINF_H +#define LLVM_LIBC_SRC_MATH_SINF_H + +namespace __llvm_libc { + +float sinf(float x); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_SINF_H diff --git a/libc/src/math/sinf.cpp b/libc/src/math/sinf.cpp new file mode 100644 --- /dev/null +++ b/libc/src/math/sinf.cpp @@ -0,0 +1,67 @@ +//===------------------- Single-precision sin function --------------------===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#include "math_utils.h" +#include "sincosf_utils.h" + +#include "include/math.h" +#include "src/__support/common.h" + +#include + +namespace __llvm_libc { + +// Fast sinf implementation. Worst-case ULP is 0.5607, maximum relative +// error is 0.5303 * 2^-23. A single-step range reduction is used for +// small values. Large inputs have their range reduced using fast integer +// arithmetic. +float LLVM_LIBC_ENTRYPOINT(sinf)(float y) { + double x = y; + double s; + int n; + const sincos_t *p = &__sincosf_table[0]; + + if (abstop12(y) < abstop12(pio4)) { + s = x * x; + + if (unlikely(abstop12(y) < abstop12(0x1p-12f))) { + if (unlikely(abstop12(y) < abstop12(0x1p-126f))) + // Force underflow for tiny y. + force_eval_float(s); + return y; + } + + return sinf_poly(x, s, p, 0); + } else if (likely(abstop12(y) < abstop12(120.0f))) { + x = reduce_fast(x, p, &n); + + // Setup the signs for sin and cos. + s = p->sign[n & 3]; + + if (n & 2) + p = &__sincosf_table[1]; + + return sinf_poly(x * s, x * x, p, n); + } else if (abstop12(y) < abstop12(INFINITY)) { + uint32_t xi = asuint32_bits(y); + int sign = xi >> 31; + + x = reduce_large(xi, &n); + + // Setup signs for sin and cos - include original sign. + s = p->sign[(n + sign) & 3]; + + if ((n + sign) & 2) + p = &__sincosf_table[1]; + + return sinf_poly(x * s, x * x, p, n); + } else + return invalidf(y); +} + +} // namespace __llvm_libc