| Message ID | DB5PR08MB10301908C58D8A9DD1C847BA83490@DB5PR08MB1030.eurprd08.prod.outlook.com |
|---|---|
| State | New |
| Headers | show |
| Series | [PATCH1/2] Improve performance of sincosf | expand |
On Tue, 26 Jun 2018, Wilco Dijkstra wrote: > +/* New sincosf. */ Comments describing something as new are never appropriate; they'll always become out of date. > +extern sincos_t sincosf_table[2] attribute_hidden; I'd expect this array, and thus all pointers to it, to be const; any writable static / global data is suspect without a clear requirement for it to be writable.
Joseph Myers wrote: > > +/* New sincosf. */ > > Comments describing something as new are never appropriate; they'll > always become out of date. That comment was from early development, the 2nd patch removes it. I've removed it from patch 1 as well. >> +extern sincos_t sincosf_table[2] attribute_hidden; > > I'd expect this array, and thus all pointers to it, to be const; any > writable static / global data is suspect without a clear requirement for > it to be writable. sincos_t itself is cost, however since that seems unclear I've added it to the uses of sincos_t. I also needed to change the patch because of recent changes to math barriers, here is the updated version: This patch is a complete rewrite of sinf, cosf and sincosf. The new version is significantly faster, as well as simple and accurate. The worst-case ULP is 0.56072, maximum relative error is 0.5303p-23 over all 4 billion inputs. In non-nearest rounding modes the error is 1ULP. The algorithm uses 3 main cases: small inputs which don't need argument reduction, small inputs which need a simple range reduction and large inputs requiring complex range reduction. The code uses approximate integer comparisons to quickly decide between these cases - on some targets this may be slow, so this can be configured to use floating point comparisons. The small range reducer uses a single reduction step to handle values up to 120.0. It is fastest on targets which support inlined round instructions. The large range reducer uses integer arithmetic for simplicity. It does a 32x96 bit multiply to compute a 64-bit modulo result. This is more than accurate enough to handle the worst-case cancellation for values close to an integer multiple of PI/4. It could be further optimized, however it is already much faster than necessary. sincosf throughput gains on Cortex-A72: * |x| < 0x1p-12 : 1.6x * |x| < M_PI_4 : 1.7x * |x| < 2 * M_PI: 1.5x * |x| < 120.0 : 1.8x * |x| < Inf : 2.3x On a benchmark with significant use of sincosf the overall speedup is >33%. ChangeLog: 2018-06-29 Wilco Dijkstra <wdijkstr@arm.com> * math/Makefile: Add s_sincosf_data.c. * sysdeps/aarch64/fpu/math_private.h (roundtoint): Use round. (converttoint): Use lround. * sysdeps/ieee754/flt-32/math_config.h (PREFER_FLOAT_COMPARISON): Define. * sysdeps/ieee754/flt-32/s_sincosf.h (abstop12): Add new function. (sincosf_poly): Likewise. (reduce_small): Likewise. (reduce_large): Likewise. * sysdeps/ieee754/flt-32/s_sincosf.c (sincosf): Rewrite. * sysdeps/ieee754/flt-32/s_sincosf_data.c: New file with sincosf data. -- diff --git a/math/Makefile b/math/Makefile index 90b3b68916e12d8563ed57772013f4420415928b..e73ceb8d4e8d520c3047f63bf6c6d05856a89186 100644 --- a/math/Makefile +++ b/math/Makefile @@ -131,7 +131,7 @@ type-double-routines := branred doasin dosincos mpa mpatan2 \ # float support type-float-suffix := f type-float-routines := k_rem_pio2f math_errf e_exp2f_data e_logf_data \ - e_log2f_data e_powf_log2_data + e_log2f_data e_powf_log2_data s_sincosf_data # _Float128 support type-float128-suffix := f128 diff --git a/sysdeps/aarch64/fpu/math_private.h b/sysdeps/aarch64/fpu/math_private.h index fcd02c06543c0841cac9949c7ca052c042b3d85b..161e4cdca265edac2537eebb85a49ec13d65c7d6 100644 --- a/sysdeps/aarch64/fpu/math_private.h +++ b/sysdeps/aarch64/fpu/math_private.h @@ -21,6 +21,8 @@ #include <fenv.h> #include <fpu_control.h> +#include <stdint.h> +#include <math.h> static __always_inline void libc_feholdexcept_aarch64 (fenv_t *envp) @@ -298,26 +300,20 @@ libc_feresetround_noex_aarch64_ctx (struct rm_ctx *ctx) #define libc_feresetround_noexf_ctx libc_feresetround_noex_aarch64_ctx #define libc_feresetround_noexl_ctx libc_feresetround_noex_aarch64_ctx -/* Hack: only include the large arm_neon.h when needed. */ -#ifdef _MATH_CONFIG_H -# include <arm_neon.h> - /* ACLE intrinsics for frintn and fcvtns instructions. */ # define TOINT_INTRINSICS 1 static inline double_t roundtoint (double_t x) { - return vget_lane_f64 (vrndn_f64 (vld1_f64 (&x)), 0); + return round (x); } static inline uint64_t converttoint (double_t x) { - return vcvtnd_s64_f64 (x); + return (uint64_t) lround (x); } #endif #include_next <math_private.h> - -#endif diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h index c4def9bbc1898d71865572d0084b219034bf3805..b68fc47acc811d9e9844faac1ca71ec3b4d7de41 100644 --- a/sysdeps/ieee754/flt-32/math_config.h +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -46,6 +46,9 @@ #ifndef TOINT_SHIFT # define TOINT_SHIFT 1 #endif +#ifndef PREFER_FLOAT_COMPARISON +#define PREFER_FLOAT_COMPARISON 0 +#endif static inline uint32_t asuint (float f) diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h index 35b5eee536ef0ca55e68e3d3a091486dca39d36a..88785aa5ef895bd39b1340ff6f85257669760c63 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.h +++ b/sysdeps/ieee754/flt-32/s_sincosf.h @@ -16,6 +16,10 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <stdint.h> +#include <math.h> +#include "math_config.h" + /* Chebyshev constants for cos, range -PI/4 - PI/4. */ static const double C0 = -0x1.ffffffffe98aep-2; static const double C1 = 0x1.55555545c50c7p-5; @@ -153,3 +157,118 @@ reduced_cos (double theta, unsigned int n) } return sign * cx; } + + +/* PI * 2^-64. */ +static const double pi64 = 0x1.921FB54442D18p-62; +/* PI / 4. */ +static const double pio4 = 0x1.921FB54442D18p-1; + +typedef struct +{ + double sign[4]; + double hpi_inv, hpi, c0, c1, c2, c3, c4, s1, s2, s3; +} sincos_t; + +extern const sincos_t sincosf_table[2] attribute_hidden; + +extern const uint32_t inv_pio4[] attribute_hidden; + +/* abstop12 assumes floating point reinterpret is fast by default. + If floating point comparisons are faster, define PREFER_FLOAT_COMPARISON. */ +#if PREFER_FLOAT_COMPARISON +static inline float +abstop12 (float x) +{ + return fabsf (x); +} +#else +static inline uint32_t +abstop12 (float x) +{ + return (asuint (x) >> 20) & 0x7ff; +} +#endif + +/* 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; +} + +/* 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, + only 2 multiplies are required and the result is accurate for |X| <= 120.0. + Use round/lround if inlined, otherwise convert to int. To avoid inaccuracies + introduced by truncating negative values, compute the quadrant * 2^24. */ +static inline double +reduce_fast (double x, const sincos_t *p, int *np) +{ + double r; +#if TOINT_INTRINSICS + r = x * p->hpi_inv; + *np = converttoint (r); + return x - roundtoint (r) * p->hpi; +#else + 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/4 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 * pi64; +} diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c index d4a5a1b22c2f9a5cf607cd2d31157dd1dbd14e9a..f8fec38f0fd7c702dfe7e3a6abcc3b19f228cedd 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.c +++ b/sysdeps/ieee754/flt-32/s_sincosf.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 2017-2018 Free Software Foundation, Inc. + Copyright (C) 2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -17,9 +17,11 @@ <http://www.gnu.org/licenses/>. */ #include <errno.h> +#include <stdint.h> #include <math.h> -#include <math_private.h> +#include <math-barriers.h> #include <libm-alias-float.h> +#include "math_config.h" #include "s_sincosf.h" #ifndef SINCOSF @@ -28,141 +30,72 @@ # define SINCOSF_FUNC SINCOSF #endif +/* Fast sincosf implementation. Worst-case ULP is 0.56072, maximum relative + error is 0.5303p-23. A single-step signed range reduction is used for + small values. Large inputs have their range reduced using fast integer + arithmetic. +*/ void -SINCOSF_FUNC (float x, float *sinx, float *cosx) +SINCOSF_FUNC (float y, float *sinp, float *cosp) { - double cx; - double theta = x; - double abstheta = fabs (theta); - /* If |x|< Pi/4. */ - if (isless (abstheta, M_PI_4)) + double x = y; + double s; + int n; + const sincos_t *p = &sincosf_table[0]; + + if (abstop12 (y) < abstop12 (pio4)) + { + double x2 = x * x; + + if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-12f))) + { + /* Force underflow for tiny y. */ + if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-126f))) + math_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)) { - if (abstheta >= 0x1p-5) /* |x| >= 2^-5. */ - { - const double theta2 = theta * theta; - /* Chebyshev polynomial of the form for sin and cos. */ - cx = C3 + theta2 * C4; - cx = C2 + theta2 * cx; - cx = C1 + theta2 * cx; - cx = C0 + theta2 * cx; - cx = 1.0 + theta2 * cx; - *cosx = cx; - cx = S3 + theta2 * S4; - cx = S2 + theta2 * cx; - cx = S1 + theta2 * cx; - cx = S0 + theta2 * cx; - cx = theta + theta * theta2 * cx; - *sinx = cx; - } - else if (abstheta >= 0x1p-27) /* |x| >= 2^-27. */ - { - /* A simpler Chebyshev approximation is close enough for this range: - for sin: x+x^3*(SS0+x^2*SS1) - for cos: 1.0+x^2*(CC0+x^3*CC1). */ - const double theta2 = theta * theta; - cx = CC0 + theta * theta2 * CC1; - cx = 1.0 + theta2 * cx; - *cosx = cx; - cx = SS0 + theta2 * SS1; - cx = theta + theta * theta2 * cx; - *sinx = cx; - } - else - { - /* Handle some special cases. */ - if (theta) - *sinx = theta - (theta * SMALL); - else - *sinx = theta; - *cosx = 1.0 - abstheta; - } + 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 /* |x| >= Pi/4. */ + else if (__glibc_likely (abstop12 (y) < abstop12 (INFINITY))) { - unsigned int signbit = isless (x, 0); - if (isless (abstheta, 9 * M_PI_4)) /* |x| < 9*Pi/4. */ - { - /* There are cases where FE_UPWARD rounding mode can - produce a result of abstheta * inv_PI_4 == 9, - where abstheta < 9pi/4, so the domain for - pio2_table must go to 5 (9 / 2 + 1). */ - unsigned int n = (abstheta * inv_PI_4) + 1; - theta = abstheta - pio2_table[n / 2]; - *sinx = reduced_sin (theta, n, signbit); - *cosx = reduced_cos (theta, n); - } - else if (isless (abstheta, INFINITY)) - { - if (abstheta < 0x1p+23) /* |x| < 2^23. */ - { - unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; - double x = n / 2; - theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; - /* Argument reduction needed. */ - *sinx = reduced_sin (theta, n, signbit); - *cosx = reduced_cos (theta, n); - } - else /* |x| >= 2^23. */ - { - x = fabsf (x); - int exponent; - GET_FLOAT_WORD (exponent, x); - exponent - = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS; - exponent += 3; - exponent /= 28; - double a = invpio4_table[exponent] * x; - double b = invpio4_table[exponent + 1] * x; - double c = invpio4_table[exponent + 2] * x; - double d = invpio4_table[exponent + 3] * x; - uint64_t l = a; - l &= ~0x7; - a -= l; - double e = a + b; - l = e; - e = a - l; - if (l & 1) - { - e -= 1.0; - e += b; - e += c; - e += d; - e *= M_PI_4; - *sinx = reduced_sin (e, l + 1, signbit); - *cosx = reduced_cos (e, l + 1); - } - else - { - e += b; - e += c; - e += d; - if (e <= 1.0) - { - e *= M_PI_4; - *sinx = reduced_sin (e, l + 1, signbit); - *cosx = reduced_cos (e, l + 1); - } - else - { - l++; - e -= 2.0; - e *= M_PI_4; - *sinx = reduced_sin (e, l + 1, signbit); - *cosx = reduced_cos (e, l + 1); - } - } - } - } - else - { - int32_t ix; - /* High word of x. */ - GET_FLOAT_WORD (ix, abstheta); - /* sin/cos(Inf or NaN) is NaN. */ - *sinx = *cosx = x - x; - if (ix == 0x7f800000) - __set_errno (EDOM); - } + 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/sysdeps/ieee754/flt-32/s_sincosf_data.c b/sysdeps/ieee754/flt-32/s_sincosf_data.c new file mode 100644 index 0000000000000000000000000000000000000000..c3dac5fe2ead99e279df26e5d2bc5a46d8d425d2 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_sincosf_data.c @@ -0,0 +1,74 @@ +/* Compute sine and cosine of argument. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <stdint.h> +#include <math.h> +#include "math_config.h" +#include "s_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 +};
On Fri, 29 Jun 2018, Wilco Dijkstra wrote: > * math/Makefile: Add s_sincosf_data.c. You then need to add a dummy s_sincosf_data.c for architectures with their own sincosf implementation, to avoid wasting space in libm.so with unused data from s_sincosf_data.c. That looks like ia64, m68k, x86_64. (i686 always builds in the generic version in case of a lack of SSE2 support, so no change is needed there, though I've no idea of the performance merits of the SSE2 version for i686 versus building the new generic one with SSE2 enabled.) > -/* Hack: only include the large arm_neon.h when needed. */ > -#ifdef _MATH_CONFIG_H > -# include <arm_neon.h> > - > /* ACLE intrinsics for frintn and fcvtns instructions. */ > # define TOINT_INTRINSICS 1 Since you're removing the #ifdef you also need to remove the space between '#' and define (update the preprocessor indentation on what was conditional code). > #endif > > #include_next <math_private.h> > - > -#endif So you seem to be moving the #include_next outside the multiple-include guards, so probably defeating the multiple-include optimization on this header. Why? > +#ifndef PREFER_FLOAT_COMPARISON > +#define PREFER_FLOAT_COMPARISON 0 > +#endif Missing preprocessor indentation, "# define". > +/* 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] = This should be __sincosf_table. I'd have expected you to get linknamespace errors from testing your second patch with the name sincosf_table, which is in the user's namespace. (They wouldn't have shown up with the first patch because sincosf isn't a standard function, but sinf and cosf are standard functions.)
Joseph Myers wrote: > You then need to add a dummy s_sincosf_data.c for architectures with their > own sincosf implementation, to avoid wasting space in libm.so with unused > data from s_sincosf_data.c. That looks like ia64, m68k, x86_64. It's about 256 bytes so it's not a huge waste, I'll do it as a separate enhancement when I get back from holiday. Here is v3 with spaces and namespace fixed: This patch is a complete rewrite of sinf, cosf and sincosf. The new version is significantly faster, as well as simple and accurate. The worst-case ULP is 0.56072, maximum relative error is 0.5303p-23 over all 4 billion inputs. In non-nearest rounding modes the error is 1ULP. The algorithm uses 3 main cases: small inputs which don't need argument reduction, small inputs which need a simple range reduction and large inputs requiring complex range reduction. The code uses approximate integer comparisons to quickly decide between these cases - on some targets this may be slow, so this can be configured to use floating point comparisons. The small range reducer uses a single reduction step to handle values up to 120.0. It is fastest on targets which support inlined round instructions. The large range reducer uses integer arithmetic for simplicity. It does a 32x96 bit multiply to compute a 64-bit modulo result. This is more than accurate enough to handle the worst-case cancellation for values close to an integer multiple of PI/4. It could be further optimized, however it is already much faster than necessary. sincosf throughput gains on Cortex-A72: * |x| < 0x1p-12 : 1.6x * |x| < M_PI_4 : 1.7x * |x| < 2 * M_PI: 1.5x * |x| < 120.0 : 1.8x * |x| < Inf : 2.3x On a benchmark with significant use of sincosf the overall speedup is >33%. ChangeLog: 2018-06-30 Wilco Dijkstra <wdijkstr@arm.com> * math/Makefile: Add s_sincosf_data.c. * sysdeps/aarch64/fpu/math_private.h (roundtoint): Use round. (converttoint): Use lround. * sysdeps/ieee754/flt-32/math_config.h (PREFER_FLOAT_COMPARISON): Define. * sysdeps/ieee754/flt-32/s_sincosf.h (abstop12): Add new function. (sincosf_poly): Likewise. (reduce_small): Likewise. (reduce_large): Likewise. * sysdeps/ieee754/flt-32/s_sincosf.c (sincosf): Rewrite. * sysdeps/ieee754/flt-32/s_sincosf_data.c: New file with sincosf data. -- diff --git a/math/Makefile b/math/Makefile index 90b3b68916e12d8563ed57772013f4420415928b..e73ceb8d4e8d520c3047f63bf6c6d05856a89186 100644 --- a/math/Makefile +++ b/math/Makefile @@ -131,7 +131,7 @@ type-double-routines := branred doasin dosincos mpa mpatan2 \ # float support type-float-suffix := f type-float-routines := k_rem_pio2f math_errf e_exp2f_data e_logf_data \ - e_log2f_data e_powf_log2_data + e_log2f_data e_powf_log2_data s_sincosf_data # _Float128 support type-float128-suffix := f128 diff --git a/sysdeps/aarch64/fpu/math_private.h b/sysdeps/aarch64/fpu/math_private.h index fcd02c06543c0841cac9949c7ca052c042b3d85b..657f654176760125b96348004129af9b544b7206 100644 --- a/sysdeps/aarch64/fpu/math_private.h +++ b/sysdeps/aarch64/fpu/math_private.h @@ -21,6 +21,8 @@ #include <fenv.h> #include <fpu_control.h> +#include <stdint.h> +#include <math.h> static __always_inline void libc_feholdexcept_aarch64 (fenv_t *envp) @@ -298,25 +300,20 @@ libc_feresetround_noex_aarch64_ctx (struct rm_ctx *ctx) #define libc_feresetround_noexf_ctx libc_feresetround_noex_aarch64_ctx #define libc_feresetround_noexl_ctx libc_feresetround_noex_aarch64_ctx -/* Hack: only include the large arm_neon.h when needed. */ -#ifdef _MATH_CONFIG_H -# include <arm_neon.h> - -/* ACLE intrinsics for frintn and fcvtns instructions. */ -# define TOINT_INTRINSICS 1 +/* Use inline round and lround instructions. */ +#define TOINT_INTRINSICS 1 static inline double_t roundtoint (double_t x) { - return vget_lane_f64 (vrndn_f64 (vld1_f64 (&x)), 0); + return round (x); } static inline uint64_t converttoint (double_t x) { - return vcvtnd_s64_f64 (x); + return (uint64_t) lround (x); } -#endif #include_next <math_private.h> diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h index c4def9bbc1898d71865572d0084b219034bf3805..6729ea06c9d687452e9891874b2c022348d274e7 100644 --- a/sysdeps/ieee754/flt-32/math_config.h +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -46,6 +46,9 @@ #ifndef TOINT_SHIFT # define TOINT_SHIFT 1 #endif +#ifndef PREFER_FLOAT_COMPARISON +# define PREFER_FLOAT_COMPARISON 0 +#endif static inline uint32_t asuint (float f) diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h index 35b5eee536ef0ca55e68e3d3a091486dca39d36a..f1fcf4ce5fa35459808c32d7461f83036716b096 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.h +++ b/sysdeps/ieee754/flt-32/s_sincosf.h @@ -16,6 +16,10 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <stdint.h> +#include <math.h> +#include "math_config.h" + /* Chebyshev constants for cos, range -PI/4 - PI/4. */ static const double C0 = -0x1.ffffffffe98aep-2; static const double C1 = 0x1.55555545c50c7p-5; @@ -153,3 +157,118 @@ reduced_cos (double theta, unsigned int n) } return sign * cx; } + + +/* PI * 2^-64. */ +static const double pi64 = 0x1.921FB54442D18p-62; +/* PI / 4. */ +static const double pio4 = 0x1.921FB54442D18p-1; + +typedef struct +{ + double sign[4]; + double hpi_inv, hpi, c0, c1, c2, c3, c4, s1, s2, s3; +} sincos_t; + +extern const sincos_t __sincosf_table[2] attribute_hidden; + +extern const uint32_t __inv_pio4[] attribute_hidden; + +/* abstop12 assumes floating point reinterpret is fast by default. + If floating point comparisons are faster, define PREFER_FLOAT_COMPARISON. */ +#if PREFER_FLOAT_COMPARISON +static inline float +abstop12 (float x) +{ + return fabsf (x); +} +#else +static inline uint32_t +abstop12 (float x) +{ + return (asuint (x) >> 20) & 0x7ff; +} +#endif + +/* 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; +} + +/* 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, + only 2 multiplies are required and the result is accurate for |X| <= 120.0. + Use round/lround if inlined, otherwise convert to int. To avoid inaccuracies + introduced by truncating negative values, compute the quadrant * 2^24. */ +static inline double +reduce_fast (double x, const sincos_t *p, int *np) +{ + double r; +#if TOINT_INTRINSICS + r = x * p->hpi_inv; + *np = converttoint (r); + return x - roundtoint (r) * p->hpi; +#else + 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/4 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 * pi64; +} diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c index d4a5a1b22c2f9a5cf607cd2d31157dd1dbd14e9a..cb3c379f882d6b29b46f32684e57426350770d9a 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.c +++ b/sysdeps/ieee754/flt-32/s_sincosf.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 2017-2018 Free Software Foundation, Inc. + Copyright (C) 2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -17,9 +17,11 @@ <http://www.gnu.org/licenses/>. */ #include <errno.h> +#include <stdint.h> #include <math.h> -#include <math_private.h> +#include <math-barriers.h> #include <libm-alias-float.h> +#include "math_config.h" #include "s_sincosf.h" #ifndef SINCOSF @@ -28,141 +30,72 @@ # define SINCOSF_FUNC SINCOSF #endif +/* Fast sincosf implementation. Worst-case ULP is 0.56072, maximum relative + error is 0.5303p-23. A single-step signed range reduction is used for + small values. Large inputs have their range reduced using fast integer + arithmetic. +*/ void -SINCOSF_FUNC (float x, float *sinx, float *cosx) +SINCOSF_FUNC (float y, float *sinp, float *cosp) { - double cx; - double theta = x; - double abstheta = fabs (theta); - /* If |x|< Pi/4. */ - if (isless (abstheta, M_PI_4)) + double x = y; + double s; + int n; + const sincos_t *p = &__sincosf_table[0]; + + if (abstop12 (y) < abstop12 (pio4)) + { + double x2 = x * x; + + if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-12f))) + { + /* Force underflow for tiny y. */ + if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-126f))) + math_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)) { - if (abstheta >= 0x1p-5) /* |x| >= 2^-5. */ - { - const double theta2 = theta * theta; - /* Chebyshev polynomial of the form for sin and cos. */ - cx = C3 + theta2 * C4; - cx = C2 + theta2 * cx; - cx = C1 + theta2 * cx; - cx = C0 + theta2 * cx; - cx = 1.0 + theta2 * cx; - *cosx = cx; - cx = S3 + theta2 * S4; - cx = S2 + theta2 * cx; - cx = S1 + theta2 * cx; - cx = S0 + theta2 * cx; - cx = theta + theta * theta2 * cx; - *sinx = cx; - } - else if (abstheta >= 0x1p-27) /* |x| >= 2^-27. */ - { - /* A simpler Chebyshev approximation is close enough for this range: - for sin: x+x^3*(SS0+x^2*SS1) - for cos: 1.0+x^2*(CC0+x^3*CC1). */ - const double theta2 = theta * theta; - cx = CC0 + theta * theta2 * CC1; - cx = 1.0 + theta2 * cx; - *cosx = cx; - cx = SS0 + theta2 * SS1; - cx = theta + theta * theta2 * cx; - *sinx = cx; - } - else - { - /* Handle some special cases. */ - if (theta) - *sinx = theta - (theta * SMALL); - else - *sinx = theta; - *cosx = 1.0 - abstheta; - } + 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 /* |x| >= Pi/4. */ + else if (__glibc_likely (abstop12 (y) < abstop12 (INFINITY))) { - unsigned int signbit = isless (x, 0); - if (isless (abstheta, 9 * M_PI_4)) /* |x| < 9*Pi/4. */ - { - /* There are cases where FE_UPWARD rounding mode can - produce a result of abstheta * inv_PI_4 == 9, - where abstheta < 9pi/4, so the domain for - pio2_table must go to 5 (9 / 2 + 1). */ - unsigned int n = (abstheta * inv_PI_4) + 1; - theta = abstheta - pio2_table[n / 2]; - *sinx = reduced_sin (theta, n, signbit); - *cosx = reduced_cos (theta, n); - } - else if (isless (abstheta, INFINITY)) - { - if (abstheta < 0x1p+23) /* |x| < 2^23. */ - { - unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; - double x = n / 2; - theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; - /* Argument reduction needed. */ - *sinx = reduced_sin (theta, n, signbit); - *cosx = reduced_cos (theta, n); - } - else /* |x| >= 2^23. */ - { - x = fabsf (x); - int exponent; - GET_FLOAT_WORD (exponent, x); - exponent - = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS; - exponent += 3; - exponent /= 28; - double a = invpio4_table[exponent] * x; - double b = invpio4_table[exponent + 1] * x; - double c = invpio4_table[exponent + 2] * x; - double d = invpio4_table[exponent + 3] * x; - uint64_t l = a; - l &= ~0x7; - a -= l; - double e = a + b; - l = e; - e = a - l; - if (l & 1) - { - e -= 1.0; - e += b; - e += c; - e += d; - e *= M_PI_4; - *sinx = reduced_sin (e, l + 1, signbit); - *cosx = reduced_cos (e, l + 1); - } - else - { - e += b; - e += c; - e += d; - if (e <= 1.0) - { - e *= M_PI_4; - *sinx = reduced_sin (e, l + 1, signbit); - *cosx = reduced_cos (e, l + 1); - } - else - { - l++; - e -= 2.0; - e *= M_PI_4; - *sinx = reduced_sin (e, l + 1, signbit); - *cosx = reduced_cos (e, l + 1); - } - } - } - } - else - { - int32_t ix; - /* High word of x. */ - GET_FLOAT_WORD (ix, abstheta); - /* sin/cos(Inf or NaN) is NaN. */ - *sinx = *cosx = x - x; - if (ix == 0x7f800000) - __set_errno (EDOM); - } + 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/sysdeps/ieee754/flt-32/s_sincosf_data.c b/sysdeps/ieee754/flt-32/s_sincosf_data.c new file mode 100644 index 0000000000000000000000000000000000000000..21fc2b60f9ead7714deedb78c05529a8c58fe096 --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_sincosf_data.c @@ -0,0 +1,74 @@ +/* Compute sine and cosine of argument. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <stdint.h> +#include <math.h> +#include "math_config.h" +#include "s_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 +};
On Sat, 30 Jun 2018, Wilco Dijkstra wrote: > Joseph Myers wrote: > > > You then need to add a dummy s_sincosf_data.c for architectures with their > > own sincosf implementation, to avoid wasting space in libm.so with unused > > data from s_sincosf_data.c. That looks like ia64, m68k, x86_64. > > It's about 256 bytes so it's not a huge waste, I'll do it as a separate > enhancement when I get back from holiday. I think the patch that adds a file only relevant to some targets should be the same patch that also avoids it building any code for other targets.
diff --git a/math/Makefile b/math/Makefile index 7cc0d16b0c83562589f1ae29fe01e45a377a2c97..69180e3c84b14c5815c9321ea1eac46d1515f0a1 100644 --- a/math/Makefile +++ b/math/Makefile @@ -130,7 +130,7 @@ type-double-routines := branred doasin dosincos mpa mpatan2 \ # float support type-float-suffix := f type-float-routines := k_rem_pio2f math_errf e_exp2f_data e_logf_data \ - e_log2f_data e_powf_log2_data + e_log2f_data e_powf_log2_data s_sincosf_data # _Float128 support type-float128-suffix := f128 diff --git a/sysdeps/aarch64/fpu/math_private.h b/sysdeps/aarch64/fpu/math_private.h index d9c2d710a95ebe76b25c9dc86ed0bea729a3a6cd..1a8bb66e55124416fa8f554a2e71f0a438abd40e 100644 --- a/sysdeps/aarch64/fpu/math_private.h +++ b/sysdeps/aarch64/fpu/math_private.h @@ -21,6 +21,8 @@ #include <fenv.h> #include <fpu_control.h> +#include <stdint.h> +#include <math.h> #define math_opt_barrier(x) \ ({ __typeof (x) __x = (x); __asm ("" : "+w" (__x)); __x; }) @@ -303,26 +305,20 @@ libc_feresetround_noex_aarch64_ctx (struct rm_ctx *ctx) #define libc_feresetround_noexf_ctx libc_feresetround_noex_aarch64_ctx #define libc_feresetround_noexl_ctx libc_feresetround_noex_aarch64_ctx -/* Hack: only include the large arm_neon.h when needed. */ -#ifdef _MATH_CONFIG_H -# include <arm_neon.h> - /* ACLE intrinsics for frintn and fcvtns instructions. */ # define TOINT_INTRINSICS 1 static inline double_t roundtoint (double_t x) { - return vget_lane_f64 (vrndn_f64 (vld1_f64 (&x)), 0); + return round (x); } static inline uint64_t converttoint (double_t x) { - return vcvtnd_s64_f64 (x); + return (uint64_t) lround (x); } #endif #include_next <math_private.h> - -#endif diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h index c4def9bbc1898d71865572d0084b219034bf3805..b68fc47acc811d9e9844faac1ca71ec3b4d7de41 100644 --- a/sysdeps/ieee754/flt-32/math_config.h +++ b/sysdeps/ieee754/flt-32/math_config.h @@ -46,6 +46,9 @@ #ifndef TOINT_SHIFT # define TOINT_SHIFT 1 #endif +#ifndef PREFER_FLOAT_COMPARISON +#define PREFER_FLOAT_COMPARISON 0 +#endif static inline uint32_t asuint (float f) diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h index 35b5eee536ef0ca55e68e3d3a091486dca39d36a..9da17dc7fe3e1d494164edbf26001303343e7c34 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.h +++ b/sysdeps/ieee754/flt-32/s_sincosf.h @@ -16,6 +16,10 @@ License along with the GNU C Library; if not, see <http://www.gnu.org/licenses/>. */ +#include <stdint.h> +#include <math.h> +#include "math_config.h" + /* Chebyshev constants for cos, range -PI/4 - PI/4. */ static const double C0 = -0x1.ffffffffe98aep-2; static const double C1 = 0x1.55555545c50c7p-5; @@ -153,3 +157,118 @@ reduced_cos (double theta, unsigned int n) } return sign * cx; } + +/* New sincosf. */ + +/* PI * 2^-64. */ +static const double pi64 = 0x1.921FB54442D18p-62; +/* PI / 4. */ +static const double pio4 = 0x1.921FB54442D18p-1; + +typedef const struct +{ + double sign[4]; + double hpi_inv, hpi, c0, c1, c2, c3, c4, s1, s2, s3; +} sincos_t; + +extern sincos_t sincosf_table[2] attribute_hidden; + +extern const uint32_t inv_pio4[] attribute_hidden; + +/* abstop12 assumes floating point reinterpret is fast by default. + If floating point comparisons are faster, define PREFER_FLOAT_COMPARISON. */ +#if PREFER_FLOAT_COMPARISON +static inline float +abstop12 (float x) +{ + return fabsf (x); +} +#else +static inline uint32_t +abstop12 (float x) +{ + return (asuint (x) >> 20) & 0x7ff; +} +#endif + +/* 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, 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; +} + +/* 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, + only 2 multiplies are required and the result is accurate for |X| <= 120.0. + Use round/lround if inlined, otherwise convert to int. To avoid inaccuracies + introduced by truncating negative values, compute the quadrant * 2^24. */ +static inline double +reduce_fast (double x, sincos_t *p, int *np) +{ + double r; +#if TOINT_INTRINSICS + r = x * p->hpi_inv; + *np = converttoint (r); + return x - roundtoint (r) * p->hpi; +#else + 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/4 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 * pi64; +} diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c index d4a5a1b22c2f9a5cf607cd2d31157dd1dbd14e9a..629673cbd37ec863c96aaebf2393e0722c2aad9c 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.c +++ b/sysdeps/ieee754/flt-32/s_sincosf.c @@ -1,5 +1,5 @@ /* Compute sine and cosine of argument. - Copyright (C) 2017-2018 Free Software Foundation, Inc. + Copyright (C) 2018 Free Software Foundation, Inc. This file is part of the GNU C Library. The GNU C Library is free software; you can redistribute it and/or @@ -17,9 +17,11 @@ <http://www.gnu.org/licenses/>. */ #include <errno.h> +#include <stdint.h> #include <math.h> #include <math_private.h> #include <libm-alias-float.h> +#include "math_config.h" #include "s_sincosf.h" #ifndef SINCOSF @@ -28,141 +30,72 @@ # define SINCOSF_FUNC SINCOSF #endif +/* Fast sincosf implementation. Worst-case ULP is 0.56072, maximum relative + error is 0.5303p-23. A single-step signed range reduction is used for + small values. Large inputs have their range reduced using fast integer + arithmetic. +*/ void -SINCOSF_FUNC (float x, float *sinx, float *cosx) +SINCOSF_FUNC (float y, float *sinp, float *cosp) { - double cx; - double theta = x; - double abstheta = fabs (theta); - /* If |x|< Pi/4. */ - if (isless (abstheta, M_PI_4)) + double x = y; + double s; + int n; + sincos_t *p = &sincosf_table[0]; + + if (abstop12 (y) < abstop12 (pio4)) + { + double x2 = x * x; + + if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-12f))) + { + /* Force underflow for tiny y. */ + if (__glibc_unlikely (abstop12 (y) < abstop12 (0x1p-126f))) + math_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)) { - if (abstheta >= 0x1p-5) /* |x| >= 2^-5. */ - { - const double theta2 = theta * theta; - /* Chebyshev polynomial of the form for sin and cos. */ - cx = C3 + theta2 * C4; - cx = C2 + theta2 * cx; - cx = C1 + theta2 * cx; - cx = C0 + theta2 * cx; - cx = 1.0 + theta2 * cx; - *cosx = cx; - cx = S3 + theta2 * S4; - cx = S2 + theta2 * cx; - cx = S1 + theta2 * cx; - cx = S0 + theta2 * cx; - cx = theta + theta * theta2 * cx; - *sinx = cx; - } - else if (abstheta >= 0x1p-27) /* |x| >= 2^-27. */ - { - /* A simpler Chebyshev approximation is close enough for this range: - for sin: x+x^3*(SS0+x^2*SS1) - for cos: 1.0+x^2*(CC0+x^3*CC1). */ - const double theta2 = theta * theta; - cx = CC0 + theta * theta2 * CC1; - cx = 1.0 + theta2 * cx; - *cosx = cx; - cx = SS0 + theta2 * SS1; - cx = theta + theta * theta2 * cx; - *sinx = cx; - } - else - { - /* Handle some special cases. */ - if (theta) - *sinx = theta - (theta * SMALL); - else - *sinx = theta; - *cosx = 1.0 - abstheta; - } + 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 /* |x| >= Pi/4. */ + else if (__glibc_likely (abstop12 (y) < abstop12 (INFINITY))) { - unsigned int signbit = isless (x, 0); - if (isless (abstheta, 9 * M_PI_4)) /* |x| < 9*Pi/4. */ - { - /* There are cases where FE_UPWARD rounding mode can - produce a result of abstheta * inv_PI_4 == 9, - where abstheta < 9pi/4, so the domain for - pio2_table must go to 5 (9 / 2 + 1). */ - unsigned int n = (abstheta * inv_PI_4) + 1; - theta = abstheta - pio2_table[n / 2]; - *sinx = reduced_sin (theta, n, signbit); - *cosx = reduced_cos (theta, n); - } - else if (isless (abstheta, INFINITY)) - { - if (abstheta < 0x1p+23) /* |x| < 2^23. */ - { - unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1; - double x = n / 2; - theta = (abstheta - x * PI_2_hi) - x * PI_2_lo; - /* Argument reduction needed. */ - *sinx = reduced_sin (theta, n, signbit); - *cosx = reduced_cos (theta, n); - } - else /* |x| >= 2^23. */ - { - x = fabsf (x); - int exponent; - GET_FLOAT_WORD (exponent, x); - exponent - = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS; - exponent += 3; - exponent /= 28; - double a = invpio4_table[exponent] * x; - double b = invpio4_table[exponent + 1] * x; - double c = invpio4_table[exponent + 2] * x; - double d = invpio4_table[exponent + 3] * x; - uint64_t l = a; - l &= ~0x7; - a -= l; - double e = a + b; - l = e; - e = a - l; - if (l & 1) - { - e -= 1.0; - e += b; - e += c; - e += d; - e *= M_PI_4; - *sinx = reduced_sin (e, l + 1, signbit); - *cosx = reduced_cos (e, l + 1); - } - else - { - e += b; - e += c; - e += d; - if (e <= 1.0) - { - e *= M_PI_4; - *sinx = reduced_sin (e, l + 1, signbit); - *cosx = reduced_cos (e, l + 1); - } - else - { - l++; - e -= 2.0; - e *= M_PI_4; - *sinx = reduced_sin (e, l + 1, signbit); - *cosx = reduced_cos (e, l + 1); - } - } - } - } - else - { - int32_t ix; - /* High word of x. */ - GET_FLOAT_WORD (ix, abstheta); - /* sin/cos(Inf or NaN) is NaN. */ - *sinx = *cosx = x - x; - if (ix == 0x7f800000) - __set_errno (EDOM); - } + 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/sysdeps/ieee754/flt-32/s_sincosf_data.c b/sysdeps/ieee754/flt-32/s_sincosf_data.c new file mode 100644 index 0000000000000000000000000000000000000000..02eacf29ddf55ea93aba250a45e0569303239b3a --- /dev/null +++ b/sysdeps/ieee754/flt-32/s_sincosf_data.c @@ -0,0 +1,74 @@ +/* Compute sine and cosine of argument. + Copyright (C) 2018 Free Software Foundation, Inc. + This file is part of the GNU C Library. + + The GNU C Library is free software; you can redistribute it and/or + modify it under the terms of the GNU Lesser General Public + License as published by the Free Software Foundation; either + version 2.1 of the License, or (at your option) any later version. + + The GNU C Library is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU + Lesser General Public License for more details. + + You should have received a copy of the GNU Lesser General Public + License along with the GNU C Library; if not, see + <http://www.gnu.org/licenses/>. */ + +#include <stdint.h> +#include <math.h> +#include "math_config.h" +#include "s_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. */ +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 +};