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/AOR_v20.02/math/test/testcases/directed/cosf.tst b/libc/AOR_v20.02/math/test/testcases/directed/cosf.tst deleted file mode 100644 --- a/libc/AOR_v20.02/math/test/testcases/directed/cosf.tst +++ /dev/null @@ -1,26 +0,0 @@ -; cosf.tst - Directed test cases for SP cosine -; -; 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 - -func=cosf op1=7fc00001 result=7fc00001 errno=0 -func=cosf op1=ffc00001 result=7fc00001 errno=0 -func=cosf op1=7f800001 result=7fc00001 errno=0 status=i -func=cosf op1=ff800001 result=7fc00001 errno=0 status=i -func=cosf op1=7f800000 result=7fc00001 errno=EDOM status=i -func=cosf op1=ff800000 result=7fc00001 errno=EDOM status=i -func=cosf op1=00000000 result=3f800000 errno=0 -func=cosf op1=80000000 result=3f800000 errno=0 -; SDCOMP-26094: check cosf in the cases for which the range reducer -; returns values furthest beyond its nominal upper bound of pi/4. -func=cosf op1=46427f1b result=3f34dc5c.565 error=0 -func=cosf op1=4647e568 result=3f34dc33.c1f error=0 -func=cosf op1=46428bac result=bf34dbf2.8e3 error=0 -func=cosf op1=4647f1f9 result=bf34dbc9.f9b error=0 -func=cosf op1=4647fe8a result=3f34db60.313 error=0 -func=cosf op1=45d8d7f1 result=bf35006a.7fd error=0 -func=cosf op1=45d371a4 result=3f350056.39b error=0 -func=cosf op1=45ce0b57 result=bf350041.f38 error=0 -func=cosf op1=45d35882 result=bf34ffec.868 error=0 -func=cosf op1=45cdf235 result=3f34ffd8.404 error=0 diff --git a/libc/AOR_v20.02/math/test/testcases/directed/sincosf.tst b/libc/AOR_v20.02/math/test/testcases/directed/sincosf.tst deleted file mode 100644 --- a/libc/AOR_v20.02/math/test/testcases/directed/sincosf.tst +++ /dev/null @@ -1,52 +0,0 @@ -; Directed test cases for SP sincos -; -; 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 - - -func=sincosf_sinf op1=7fc00001 result=7fc00001 errno=0 -func=sincosf_sinf op1=ffc00001 result=7fc00001 errno=0 -func=sincosf_sinf op1=7f800001 result=7fc00001 errno=0 status=i -func=sincosf_sinf op1=ff800001 result=7fc00001 errno=0 status=i -func=sincosf_sinf op1=7f800000 result=7fc00001 errno=EDOM status=i -func=sincosf_sinf op1=ff800000 result=7fc00001 errno=EDOM status=i -func=sincosf_sinf op1=00000000 result=00000000 errno=0 -func=sincosf_sinf op1=80000000 result=80000000 errno=0 -func=sincosf_sinf op1=c70d39a1 result=be37fad5.7ed errno=0 -func=sincosf_sinf op1=46427f1b result=3f352d80.f9b error=0 -func=sincosf_sinf op1=4647e568 result=3f352da9.7be error=0 -func=sincosf_sinf op1=46428bac result=bf352dea.924 error=0 -func=sincosf_sinf op1=4647f1f9 result=bf352e13.146 error=0 -func=sincosf_sinf op1=4647fe8a result=3f352e7c.ac9 error=0 -func=sincosf_sinf op1=45d8d7f1 result=3f35097b.cb0 error=0 -func=sincosf_sinf op1=45d371a4 result=bf350990.102 error=0 -func=sincosf_sinf op1=45ce0b57 result=3f3509a4.554 error=0 -func=sincosf_sinf op1=45d35882 result=3f3509f9.bdb error=0 -func=sincosf_sinf op1=45cdf235 result=bf350a0e.02c error=0 - -func=sincosf_cosf op1=7fc00001 result=7fc00001 errno=0 -func=sincosf_cosf op1=ffc00001 result=7fc00001 errno=0 -func=sincosf_cosf op1=7f800001 result=7fc00001 errno=0 status=i -func=sincosf_cosf op1=ff800001 result=7fc00001 errno=0 status=i -func=sincosf_cosf op1=7f800000 result=7fc00001 errno=EDOM status=i -func=sincosf_cosf op1=ff800000 result=7fc00001 errno=EDOM status=i -func=sincosf_cosf op1=00000000 result=3f800000 errno=0 -func=sincosf_cosf op1=80000000 result=3f800000 errno=0 -func=sincosf_cosf op1=46427f1b result=3f34dc5c.565 error=0 -func=sincosf_cosf op1=4647e568 result=3f34dc33.c1f error=0 -func=sincosf_cosf op1=46428bac result=bf34dbf2.8e3 error=0 -func=sincosf_cosf op1=4647f1f9 result=bf34dbc9.f9b error=0 -func=sincosf_cosf op1=4647fe8a result=3f34db60.313 error=0 -func=sincosf_cosf op1=45d8d7f1 result=bf35006a.7fd error=0 -func=sincosf_cosf op1=45d371a4 result=3f350056.39b error=0 -func=sincosf_cosf op1=45ce0b57 result=bf350041.f38 error=0 -func=sincosf_cosf op1=45d35882 result=bf34ffec.868 error=0 -func=sincosf_cosf op1=45cdf235 result=3f34ffd8.404 error=0 - -; no underflow -func=sincosf_sinf op1=17800000 result=17800000.000 -func=sincosf_cosf op1=17800000 result=3f800000.000 -; underflow -func=sincosf_sinf op1=00400000 result=00400000.000 status=ux -func=sincosf_cosf op1=00400000 result=3f800000.000 status=ux diff --git a/libc/AOR_v20.02/math/test/testcases/directed/sinf.tst b/libc/AOR_v20.02/math/test/testcases/directed/sinf.tst deleted file mode 100644 --- a/libc/AOR_v20.02/math/test/testcases/directed/sinf.tst +++ /dev/null @@ -1,29 +0,0 @@ -; sinf.tst - Directed test cases for SP sine -; -; 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 - - -func=sinf op1=7fc00001 result=7fc00001 errno=0 -func=sinf op1=ffc00001 result=7fc00001 errno=0 -func=sinf op1=7f800001 result=7fc00001 errno=0 status=i -func=sinf op1=ff800001 result=7fc00001 errno=0 status=i -func=sinf op1=7f800000 result=7fc00001 errno=EDOM status=i -func=sinf op1=ff800000 result=7fc00001 errno=EDOM status=i -func=sinf op1=00000000 result=00000000 errno=0 -func=sinf op1=80000000 result=80000000 errno=0 -; Directed test for a failure I found while developing mathbench -func=sinf op1=c70d39a1 result=be37fad5.7ed errno=0 -; SDCOMP-26094: check sinf in the cases for which the range reducer -; returns values furthest beyond its nominal upper bound of pi/4. -func=sinf op1=46427f1b result=3f352d80.f9b error=0 -func=sinf op1=4647e568 result=3f352da9.7be error=0 -func=sinf op1=46428bac result=bf352dea.924 error=0 -func=sinf op1=4647f1f9 result=bf352e13.146 error=0 -func=sinf op1=4647fe8a result=3f352e7c.ac9 error=0 -func=sinf op1=45d8d7f1 result=3f35097b.cb0 error=0 -func=sinf op1=45d371a4 result=bf350990.102 error=0 -func=sinf op1=45ce0b57 result=3f3509a4.554 error=0 -func=sinf op1=45d35882 result=3f3509f9.bdb error=0 -func=sinf op1=45cdf235 result=bf350a0e.02c error=0 diff --git a/libc/AOR_v20.02/math/test/testcases/random/float.tst b/libc/AOR_v20.02/math/test/testcases/random/float.tst --- a/libc/AOR_v20.02/math/test/testcases/random/float.tst +++ b/libc/AOR_v20.02/math/test/testcases/random/float.tst @@ -4,10 +4,6 @@ !! See https://llvm.org/LICENSE.txt for license information. !! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception -test sinf 10000 -test cosf 10000 -test sincosf_sinf 5000 -test sincosf_cosf 5000 test tanf 10000 test expf 10000 test exp2f 10000 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" @@ -123,7 +124,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 + libc.src.math.cosf libc.src.math.round + libc.src.math.sincosf + libc.src.math.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,23 @@ +add_header_library( + math_utils + HDRS + math_utils.h + DEPENDS + libc.include.errno + libc.include.math + libc.src.errno.__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 +32,39 @@ SRC round_redirector.cpp ) + +add_entrypoint_object( + cosf + SRCS + cosf.cpp + HDRS + cosf.h + DEPENDS + .sincosf_utils + libc.include.math + libc.src.errno.__errno_location +) + +add_entrypoint_object( + sinf + SRCS + sinf.cpp + HDRS + sinf.h + DEPENDS + .sincosf_utils + libc.include.math + libc.src.errno.__errno_location +) + +add_entrypoint_object( + sincosf + SRCS + sincosf.cpp + HDRS + sincosf.h + DEPENDS + .sincosf_utils + libc.include.math + libc.src.errno.__errno_location +) 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 --------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC_MATH_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(as_float(0x39800000)))) + 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 = as_uint32_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); + } + + 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,49 @@ +//===-- Collection of utils for implementing math functions -----*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC_MATH_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 { + +static inline float with_errnof(float x, int err) { + if (math_errhandling & MATH_ERRNO) + llvmlibc_errno = err; + return x; +} + +static inline uint32_t as_uint32_bits(float x) { + return *reinterpret_cast(&x); +} + +static inline float as_float(uint32_t x) { + return *reinterpret_cast(&x); +} + +static inline double as_double(uint64_t x) { + return *reinterpret_cast(&x); +} + +static inline constexpr 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 sincosf -----------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC_MATH_SINCOSF_H +#define LLVM_LIBC_SRC_MATH_SINCOSF_H + +namespace __llvm_libc { + +void sincosf(float x, float *sinx, float *cosx); + +} // namespace __llvm_libc + +#endif // LLVM_LIBC_SRC_MATH_SINCOSF_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/math.h" +#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(as_float(0x39800000)))) { + if (unlikely(abstop12(y) < abstop12(as_float(0x800000)))) + // 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 = as_uint32_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,51 @@ +//===-- 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}, + as_double(0x41645f306dc9c883), + as_double(0x3ff921fb54442d18), + as_double(0x3ff0000000000000), + as_double(0xbfdffffffd0c621c), + as_double(0x3fa55553e1068f19), + as_double(0xbf56c087e89a359d), + as_double(0x3ef99343027bf8c3), + as_double(0xbfc555545995a603), + as_double(0x3f81107605230bc4), + as_double(0xbf2994eb3774cf24)}, + {{1.0, -1.0, -1.0, 1.0}, + as_double(0x41645f306dc9c883), + as_double(0x3ff921fb54442d18), + as_double(0xbff0000000000000), + as_double(0x3fdffffffd0c621c), + as_double(0xbfa55553e1068f19), + as_double(0x3f56c087e89a359d), + as_double(0xbef99343027bf8c3), + as_double(0xbfc555545995a603), + as_double(0x3f81107605230bc4), + as_double(0xbf2994eb3774cf24)}, +}; + +// 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 ---------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H +#define LLVM_LIBC_SRC_MATH_SINCOSF_UTILS_H + +#include "math_utils.h" + +#include + +namespace __llvm_libc { + +// 2PI * 2^-64. +static const double pi63 = as_double(0x3c1921fb54442d18); +// PI / 4. +static const double pio4 = as_double(0x3fe921fb54442d18); + +// 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 (as_uint32_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_UTILS_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 --------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_SRC_MATH_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,68 @@ +//===-- 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(as_float(0x39800000)))) { + if (unlikely(abstop12(y) < abstop12(as_float(0x800000)))) + // 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 = as_uint32_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); + } + + return invalidf(y); +} + +} // namespace __llvm_libc diff --git a/libc/test/src/CMakeLists.txt b/libc/test/src/CMakeLists.txt --- a/libc/test/src/CMakeLists.txt +++ b/libc/test/src/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(assert) add_subdirectory(errno) +add_subdirectory(math) add_subdirectory(signal) add_subdirectory(stdio) add_subdirectory(stdlib) diff --git a/libc/test/src/math/CMakeLists.txt b/libc/test/src/math/CMakeLists.txt new file mode 100644 --- /dev/null +++ b/libc/test/src/math/CMakeLists.txt @@ -0,0 +1,80 @@ +add_libc_testsuite(libc_math_unittests) + +function(add_math_unittest name) + cmake_parse_arguments( + "MATH_UNITTEST" + "NEED_MPFR" # No optional arguments + "" # Single value arguments + "" # Multi-value arguments + ${ARGN} + ) + + if(MATH_UNITTEST_NEED_MPFR) + if(NOT LIBC_TESTS_CAN_USE_MPFR) + message("WARNING: Math test ${name} will be skipped as MPFR library is not available.") + return() + endif() + endif() + + add_libc_unittest(${name} ${MATH_UNITTEST_UNPARSED_ARGUMENTS}) + if(MATH_UNITTEST_NEED_MPFR) + get_fq_target_name(${name} fq_target_name) + target_link_libraries(${fq_target_name} PRIVATE libcMPFRWrapper -lmpfr -lgmp) + endif() +endfunction(add_math_unittest) + +add_header_library( + float_utils + HDRS + float.h +) + +# TODO(sivachandra): Remove the dependency on __errno_location as the tested +# entry points depend on the already. +add_math_unittest( + cosf_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + cosf_test.cpp + HDRS + sdcomp26094.h + DEPENDS + .float_utils + libc.src.errno.__errno_location + libc.src.math.cosf + libc.utils.CPP.standalone_cpp +) + +add_math_unittest( + sinf_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + sinf_test.cpp + HDRS + sdcomp26094.h + DEPENDS + .float_utils + libc.src.errno.__errno_location + libc.src.math.sinf + libc.utils.CPP.standalone_cpp +) + +add_math_unittest( + sincosf_test + NEED_MPFR + SUITE + libc_math_unittests + SRCS + sincosf_test.cpp + HDRS + sdcomp26094.h + DEPENDS + .float_utils + libc.src.errno.__errno_location + libc.src.math.sincosf + libc.utils.CPP.standalone_cpp +) diff --git a/libc/test/src/math/cosf_test.cpp b/libc/test/src/math/cosf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/cosf_test.cpp @@ -0,0 +1,103 @@ +//===-- Unittests 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 +// +//===----------------------------------------------------------------------===// + +#include "include/math.h" +#include "src/errno/llvmlibc_errno.h" +#include "src/math/cosf.h" +#include "src/math/math_utils.h" +#include "test/src/math/float.h" +#include "test/src/math/sdcomp26094.h" +#include "utils/CPP/Array.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/Test.h" + +#include + +using __llvm_libc::as_float; +using __llvm_libc::as_uint32_bits; + +using __llvm_libc::testing::FloatBits; +using __llvm_libc::testing::sdcomp26094Values; + +namespace mpfr = __llvm_libc::testing::mpfr; + +// 12 additional bits of precision over the base precision of a |float| +// value. +static constexpr mpfr::Tolerance tolerance{mpfr::Tolerance::floatPrecision, 12, + 3 * 0x1000 / 4}; + +TEST(CosfTest, SpecialNumbers) { + llvmlibc_errno = 0; + + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::QNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::NegQNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::SNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::NegSNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::One, + as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::Zero)))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::One, + as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::NegZero)))); + EXPECT_EQ(llvmlibc_errno, 0); + + llvmlibc_errno = 0; + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::Inf))))); + EXPECT_EQ(llvmlibc_errno, EDOM); + + llvmlibc_errno = 0; + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::cosf(as_float(FloatBits::NegInf))))); + EXPECT_EQ(llvmlibc_errno, EDOM); +} + +TEST(CosfTest, InFloatRange) { + constexpr uint32_t count = 1000000; + constexpr uint32_t step = UINT32_MAX / count; + for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) { + float x = as_float(v); + if (isnan(x) || isinf(x)) + continue; + EXPECT_TRUE(mpfr::equalsCos(x, __llvm_libc::cosf(x), tolerance)); + } +} + +// For small values, cos(x) is 1. +TEST(CosfTest, SmallValues) { + float x = as_float(0x17800000); + float result = __llvm_libc::cosf(x); + EXPECT_TRUE(mpfr::equalsCos(x, result, tolerance)); + EXPECT_EQ(FloatBits::One, as_uint32_bits(result)); + + x = as_float(0x00400000); + result = __llvm_libc::cosf(x); + EXPECT_TRUE(mpfr::equalsCos(x, result, tolerance)); + EXPECT_EQ(FloatBits::One, as_uint32_bits(result)); +} + +// SDCOMP-26094: check cosf in the cases for which the range reducer +// returns values furthest beyond its nominal upper bound of pi/4. +TEST(CosfTest, SDCOMP_26094) { + for (uint32_t v : sdcomp26094Values) { + float x = as_float(v); + EXPECT_TRUE(mpfr::equalsCos(x, __llvm_libc::cosf(x), tolerance)); + } +} diff --git a/libc/test/src/math/float.h b/libc/test/src/math/float.h new file mode 100644 --- /dev/null +++ b/libc/test/src/math/float.h @@ -0,0 +1,49 @@ +//===-- Single precision floating point test utils --------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_TEST_SRC_MATH_FLOAT_H +#define LLVM_LIBC_TEST_SRC_MATH_FLOAT_H + +#include "src/math/math_utils.h" + +namespace __llvm_libc { +namespace testing { + +struct FloatBits { + // The various NaN bit patterns here are just one of the many possible + // patterns. The functions isQNan and isNegQNan can help understand why. + + static const uint32_t QNan = 0x7fc00000; + static const uint32_t NegQNan = 0xffc00000; + + static const uint32_t SNan = 0x7f800001; + static const uint32_t NegSNan = 0xff800001; + + static bool isQNan(float f) { + uint32_t bits = as_uint32_bits(f); + return ((0x7fc00000 & bits) != 0) && ((0x80000000 & bits) == 0); + } + + static bool isNegQNan(float f) { + uint32_t bits = as_uint32_bits(f); + return 0xffc00000 & bits; + } + + static constexpr uint32_t Zero = 0x0; + static constexpr uint32_t NegZero = 0x80000000; + + static constexpr uint32_t Inf = 0x7f800000; + static constexpr uint32_t NegInf = 0xff800000; + + static constexpr uint32_t One = 0x3f800000; +}; + +} // namespace testing +} // namespace __llvm_libc + +#endif // LLVM_LIBC_TEST_SRC_MATH_FLOAT_H diff --git a/libc/test/src/math/sdcomp26094.h b/libc/test/src/math/sdcomp26094.h new file mode 100644 --- /dev/null +++ b/libc/test/src/math/sdcomp26094.h @@ -0,0 +1,25 @@ +//===-- SDCOMP-26094 specific items -----------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_TEST_SRC_MATH_SDCOMP26094_H +#define LLVM_LIBC_TEST_SRC_MATH_SDCOMP26094_H + +#include "utils/CPP/Array.h" + +namespace __llvm_libc { +namespace testing { + +static constexpr __llvm_libc::cpp::Array sdcomp26094Values{ + 0x46427f1b, 0x4647e568, 0x46428bac, 0x4647f1f9, 0x4647fe8a, + 0x45d8d7f1, 0x45d371a4, 0x45ce0b57, 0x45d35882, 0x45cdf235, +}; + +} // namespace testing +} // namespace __llvm_libc + +#endif // LLVM_LIBC_TEST_SRC_MATH_SDCOMP26094_H diff --git a/libc/test/src/math/sincosf_test.cpp b/libc/test/src/math/sincosf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/sincosf_test.cpp @@ -0,0 +1,125 @@ +//===-- Unittests for 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/math.h" +#include "src/errno/llvmlibc_errno.h" +#include "src/math/math_utils.h" +#include "src/math/sincosf.h" +#include "test/src/math/float.h" +#include "test/src/math/sdcomp26094.h" +#include "utils/CPP/Array.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/Test.h" + +#include + +using __llvm_libc::as_float; +using __llvm_libc::as_uint32_bits; + +using __llvm_libc::testing::FloatBits; +using __llvm_libc::testing::sdcomp26094Values; + +namespace mpfr = __llvm_libc::testing::mpfr; +// 12 additional bits of precision over the base precision of a |float| +// value. +static constexpr mpfr::Tolerance tolerance{mpfr::Tolerance::floatPrecision, 12, + 3 * 0x1000 / 4}; + +TEST(SinCosfTest, SpecialNumbers) { + llvmlibc_errno = 0; + float sin, cos; + + __llvm_libc::sincosf(as_float(FloatBits::QNan), &sin, &cos); + EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(cos))); + EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(sin))); + EXPECT_EQ(llvmlibc_errno, 0); + + __llvm_libc::sincosf(as_float(FloatBits::NegQNan), &sin, &cos); + EXPECT_TRUE(FloatBits::isNegQNan(as_uint32_bits(cos))); + EXPECT_TRUE(FloatBits::isNegQNan(as_uint32_bits(sin))); + EXPECT_EQ(llvmlibc_errno, 0); + + __llvm_libc::sincosf(as_float(FloatBits::SNan), &sin, &cos); + EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(cos))); + EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(sin))); + EXPECT_EQ(llvmlibc_errno, 0); + + __llvm_libc::sincosf(as_float(FloatBits::NegSNan), &sin, &cos); + EXPECT_TRUE(FloatBits::isNegQNan(as_uint32_bits(cos))); + EXPECT_TRUE(FloatBits::isNegQNan(as_uint32_bits(sin))); + EXPECT_EQ(llvmlibc_errno, 0); + + __llvm_libc::sincosf(as_float(FloatBits::Zero), &sin, &cos); + EXPECT_EQ(FloatBits::One, as_uint32_bits(cos)); + EXPECT_EQ(FloatBits::Zero, as_uint32_bits(sin)); + EXPECT_EQ(llvmlibc_errno, 0); + + __llvm_libc::sincosf(as_float(FloatBits::NegZero), &sin, &cos); + EXPECT_EQ(FloatBits::One, as_uint32_bits(cos)); + EXPECT_EQ(FloatBits::NegZero, as_uint32_bits(sin)); + EXPECT_EQ(llvmlibc_errno, 0); + + llvmlibc_errno = 0; + __llvm_libc::sincosf(as_float(FloatBits::Inf), &sin, &cos); + EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(cos))); + EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(sin))); + EXPECT_EQ(llvmlibc_errno, EDOM); + + llvmlibc_errno = 0; + __llvm_libc::sincosf(as_float(FloatBits::NegInf), &sin, &cos); + EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(cos))); + EXPECT_TRUE(FloatBits::isQNan(as_uint32_bits(sin))); + EXPECT_EQ(llvmlibc_errno, EDOM); +} + +TEST(SinCosfTest, InFloatRange) { + constexpr uint32_t count = 1000000; + constexpr uint32_t step = UINT32_MAX / count; + for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) { + float x = as_float(v); + if (isnan(x) || isinf(x)) + continue; + + float sin, cos; + __llvm_libc::sincosf(x, &sin, &cos); + EXPECT_TRUE(mpfr::equalsCos(x, cos, tolerance)); + EXPECT_TRUE(mpfr::equalsSin(x, sin, tolerance)); + } +} + +// For small values, cos(x) is 1 and sin(x) is x. +TEST(SinCosfTest, SmallValues) { + uint32_t bits = 0x17800000; + float x = as_float(bits); + float result_cos, result_sin; + __llvm_libc::sincosf(x, &result_sin, &result_cos); + EXPECT_TRUE(mpfr::equalsCos(x, result_cos, tolerance)); + EXPECT_TRUE(mpfr::equalsSin(x, result_sin, tolerance)); + EXPECT_EQ(FloatBits::One, as_uint32_bits(result_cos)); + EXPECT_EQ(bits, as_uint32_bits(result_sin)); + + bits = 0x00400000; + x = as_float(bits); + __llvm_libc::sincosf(x, &result_sin, &result_cos); + EXPECT_TRUE(mpfr::equalsCos(x, result_cos, tolerance)); + EXPECT_TRUE(mpfr::equalsSin(x, result_sin, tolerance)); + EXPECT_EQ(FloatBits::One, as_uint32_bits(result_cos)); + EXPECT_EQ(bits, as_uint32_bits(result_sin)); +} + +// SDCOMP-26094: check sinf in the cases for which the range reducer +// returns values furthest beyond its nominal upper bound of pi/4. +TEST(SinCosfTest, SDCOMP_26094) { + for (uint32_t v : sdcomp26094Values) { + float x = as_float(v); + float sin, cos; + __llvm_libc::sincosf(x, &sin, &cos); + EXPECT_TRUE(mpfr::equalsCos(x, cos, tolerance)); + EXPECT_TRUE(mpfr::equalsSin(x, sin, tolerance)); + } +} diff --git a/libc/test/src/math/sinf_test.cpp b/libc/test/src/math/sinf_test.cpp new file mode 100644 --- /dev/null +++ b/libc/test/src/math/sinf_test.cpp @@ -0,0 +1,110 @@ +//===-- Unittests 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 +// +//===----------------------------------------------------------------------===// + +#include "include/math.h" +#include "src/errno/llvmlibc_errno.h" +#include "src/math/math_utils.h" +#include "src/math/sinf.h" +#include "test/src/math/float.h" +#include "test/src/math/sdcomp26094.h" +#include "utils/CPP/Array.h" +#include "utils/MPFRWrapper/MPFRUtils.h" +#include "utils/UnitTest/Test.h" + +#include + +using __llvm_libc::as_float; +using __llvm_libc::as_uint32_bits; + +using __llvm_libc::testing::FloatBits; +using __llvm_libc::testing::sdcomp26094Values; + +namespace mpfr = __llvm_libc::testing::mpfr; + +// 12 additional bits of precision over the base precision of a |float| +// value. +static constexpr mpfr::Tolerance tolerance{mpfr::Tolerance::floatPrecision, 12, + 3 * 0x1000 / 4}; + +TEST(SinfTest, SpecialNumbers) { + llvmlibc_errno = 0; + + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::QNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::NegQNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::SNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::NegSNan))))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::Zero, + as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::Zero)))); + EXPECT_EQ(llvmlibc_errno, 0); + + EXPECT_EQ(FloatBits::NegZero, + as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::NegZero)))); + EXPECT_EQ(llvmlibc_errno, 0); + + llvmlibc_errno = 0; + EXPECT_TRUE(FloatBits::isQNan( + as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::Inf))))); + EXPECT_EQ(llvmlibc_errno, EDOM); + + llvmlibc_errno = 0; + EXPECT_TRUE(FloatBits::isNegQNan( + as_uint32_bits(__llvm_libc::sinf(as_float(FloatBits::NegInf))))); + EXPECT_EQ(llvmlibc_errno, EDOM); +} + +TEST(SinfTest, InFloatRange) { + constexpr uint32_t count = 1000000; + constexpr uint32_t step = UINT32_MAX / count; + for (uint32_t i = 0, v = 0; i <= count; ++i, v += step) { + float x = as_float(v); + if (isnan(x) || isinf(x)) + continue; + EXPECT_TRUE(mpfr::equalsSin(x, __llvm_libc::sinf(x), tolerance)); + } +} + +TEST(SinfTest, SpecificBitPatterns) { + float x = as_float(0xc70d39a1); + EXPECT_TRUE(mpfr::equalsSin(x, __llvm_libc::sinf(x), tolerance)); +} + +// For small values, sin(x) is x. +TEST(SinfTest, SmallValues) { + uint32_t bits = 0x17800000; + float x = as_float(bits); + float result = __llvm_libc::sinf(x); + EXPECT_TRUE(mpfr::equalsSin(x, result, tolerance)); + EXPECT_EQ(bits, as_uint32_bits(result)); + + bits = 0x00400000; + x = as_float(bits); + result = __llvm_libc::sinf(x); + EXPECT_TRUE(mpfr::equalsSin(x, result, tolerance)); + EXPECT_EQ(bits, as_uint32_bits(result)); +} + +// SDCOMP-26094: check sinf in the cases for which the range reducer +// returns values furthest beyond its nominal upper bound of pi/4. +TEST(SinfTest, SDCOMP_26094) { + for (uint32_t v : sdcomp26094Values) { + float x = as_float(v); + EXPECT_TRUE(mpfr::equalsSin(x, __llvm_libc::sinf(x), tolerance)); + } +} diff --git a/libc/utils/CMakeLists.txt b/libc/utils/CMakeLists.txt --- a/libc/utils/CMakeLists.txt +++ b/libc/utils/CMakeLists.txt @@ -1,5 +1,6 @@ add_subdirectory(CPP) add_subdirectory(HdrGen) +add_subdirectory(MPFRWrapper) add_subdirectory(testutils) add_subdirectory(UnitTest) add_subdirectory(benchmarks) diff --git a/libc/utils/MPFRWrapper/CMakeLists.txt b/libc/utils/MPFRWrapper/CMakeLists.txt new file mode 100644 --- /dev/null +++ b/libc/utils/MPFRWrapper/CMakeLists.txt @@ -0,0 +1,17 @@ +try_compile( + LIBC_TESTS_CAN_USE_MPFR + ${CMAKE_CURRENT_BINARY_DIR} + SOURCES + ${CMAKE_CURRENT_SOURCE_DIR}/check_mpfr.cpp + LINK_LIBRARIES + -lmpfr -lgmp +) + +if(LIBC_TESTS_CAN_USE_MPFR) + add_library(libcMPFRWrapper + MPFRUtils.cpp + MPFRUtils.h + ) +else() + message(WARNING "Math tests using MPFR will be skipped.") +endif() diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h new file mode 100644 --- /dev/null +++ b/libc/utils/MPFRWrapper/MPFRUtils.h @@ -0,0 +1,51 @@ +//===-- MPFRUtils.h ---------------------------------------------*- C++ -*-===// +// +// Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions. +// See https://llvm.org/LICENSE.txt for license information. +// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception +// +//===----------------------------------------------------------------------===// + +#ifndef LLVM_LIBC_UTILS_TESTUTILS_MPFRUTILS_H +#define LLVM_LIBC_UTILS_TESTUTILS_MPFRUTILS_H + +#include + +namespace __llvm_libc { +namespace testing { +namespace mpfr { + +struct Tolerance { + // Number of bits used to represent the fractional + // part of a value of type 'float'. + static constexpr unsigned int floatPrecision = 23; + + // Number of bits used to represent the fractional + // part of a value of type 'double'. + static constexpr unsigned int doublePrecision = 52; + + // The base precision of the number. For example, for values of + // type float, the base precision is the value |floatPrecision|. + unsigned int basePrecision; + + unsigned int width; // Number of valid LSB bits in |value|. + + // The bits in the tolerance value. The tolerance value will be + // sum(bits[width - i] * 2 ^ (- basePrecision - i)) for |i| in + // range [1, width]. + uint32_t bits; +}; + +// Return true if |libcOutput| is within the tolerance |t| of the cos(x) +// value as evaluated by MPFR. +bool equalsCos(float x, float libcOutput, const Tolerance &t); + +// Return true if |libcOutput| is within the tolerance |t| of the sin(x) +// value as evaluated by MPFR. +bool equalsSin(float x, float libcOutput, const Tolerance &t); + +} // namespace mpfr +} // namespace testing +} // namespace __llvm_libc + +#endif // LLVM_LIBC_UTILS_TESTUTILS_MPFRUTILS_H diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp new file mode 100644 --- /dev/null +++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp @@ -0,0 +1,97 @@ +//===-- Utils which wrap MPFR ---------------------------------------------===// +// +// 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 "MPFRUtils.h" + +#include +#include + +namespace __llvm_libc { +namespace testing { +namespace mpfr { + +class MPFRNumber { + // A precision value which allows sufficiently large additional + // precision even compared to double precision floating point values. + static constexpr unsigned int mpfrPrecision = 96; + + mpfr_t value; + +public: + MPFRNumber() { mpfr_init2(value, mpfrPrecision); } + + explicit MPFRNumber(float x) { + mpfr_init2(value, mpfrPrecision); + mpfr_set_flt(value, x, MPFR_RNDN); + } + + MPFRNumber(const MPFRNumber &other) { + mpfr_set(value, other.value, MPFR_RNDN); + } + + ~MPFRNumber() { mpfr_clear(value); } + + // Returns true if |other| is within the tolerance value |t| of this + // number. + bool isEqual(const MPFRNumber &other, const Tolerance &t) { + MPFRNumber tolerance(0.0); + uint32_t bitMask = 1 << (t.width - 1); + for (int exponent = -t.basePrecision; bitMask > 0; bitMask >>= 1) { + --exponent; + if (t.bits & bitMask) { + MPFRNumber delta; + mpfr_set_ui_2exp(delta.value, 1, exponent, MPFR_RNDN); + mpfr_add(tolerance.value, tolerance.value, delta.value, MPFR_RNDN); + } + } + + MPFRNumber difference; + if (mpfr_cmp(value, other.value) >= 0) + mpfr_sub(difference.value, value, other.value, MPFR_RNDN); + else + mpfr_sub(difference.value, other.value, value, MPFR_RNDN); + + return mpfr_lessequal_p(difference.value, tolerance.value); + } + + // These functions are useful for debugging. + float asFloat() const { return mpfr_get_flt(value, MPFR_RNDN); } + double asDouble() const { return mpfr_get_d(value, MPFR_RNDN); } + void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); } + +public: + static MPFRNumber cos(float x) { + MPFRNumber result; + MPFRNumber mpfrX(x); + mpfr_cos(result.value, mpfrX.value, MPFR_RNDN); + return result; + } + + static MPFRNumber sin(float x) { + MPFRNumber result; + MPFRNumber mpfrX(x); + mpfr_sin(result.value, mpfrX.value, MPFR_RNDN); + return result; + } +}; + +bool equalsCos(float input, float libcOutput, const Tolerance &t) { + MPFRNumber mpfrResult = MPFRNumber::cos(input); + MPFRNumber libcResult(libcOutput); + return mpfrResult.isEqual(libcResult, t); +} + +bool equalsSin(float input, float libcOutput, const Tolerance &t) { + MPFRNumber mpfrResult = MPFRNumber::sin(input); + MPFRNumber libcResult(libcOutput); + return mpfrResult.isEqual(libcResult, t); +} + +} // namespace mpfr +} // namespace testing +} // namespace __llvm_libc diff --git a/libc/utils/MPFRWrapper/check_mpfr.cpp b/libc/utils/MPFRWrapper/check_mpfr.cpp new file mode 100644 --- /dev/null +++ b/libc/utils/MPFRWrapper/check_mpfr.cpp @@ -0,0 +1,8 @@ +#include + +int main() { + mpfr_t x; + mpfr_init(x); + mpfr_clear(x); + return 0; +}