diff options
| -rw-r--r-- | ChangeLog | 11 | ||||
| -rw-r--r-- | NEWS | 4 | ||||
| -rw-r--r-- | sysdeps/x86_64/fpu/multiarch/Makefile | 5 | ||||
| -rw-r--r-- | sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c | 240 | ||||
| -rw-r--r-- | sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S | 2 | ||||
| -rw-r--r-- | sysdeps/x86_64/fpu/multiarch/s_sincosf.c | 28 | ||||
| -rw-r--r-- | sysdeps/x86_64/fpu/s_sincosf.S | 2 |
7 files changed, 288 insertions, 4 deletions
diff --git a/ChangeLog b/ChangeLog index 09ea55ea2b..ad0641e232 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,14 @@ +2018-01-08 H.J. Lu <hongjiu.lu@intel.com> + + * sysdeps/x86_64/fpu/multiarch/Makefile (libm-sysdep_routines): + Add s_sincosf-sse2 and s_sincosf-fma. + (CFLAGS-s_sincosf-fma.c): New. + * sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: New file. + * sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S: Likewise. + * sysdeps/x86_64/fpu/multiarch/s_sincosf.c: Likewise. + * sysdeps/x86_64/fpu/s_sincosf.S: Don't add alias if + __sincosf is defined. + 2018-01-08 Florian Weimer <fweimer@redhat.com> * nptl/tst-thread-exit-clobber.cc: New file. diff --git a/NEWS b/NEWS index 0d6f826dcc..75bf46776a 100644 --- a/NEWS +++ b/NEWS @@ -19,8 +19,8 @@ Major new features: supported on i386, x86_64, x32 and aarch64. * Optimized x86-64 asin, atan2, exp, expf, log, pow, atan, sin, cosf, - sinf and tan with FMA, contributed by Arjan van de Ven and H.J. Lu - from Intel. + sinf, sincosf and tan with FMA, contributed by Arjan van de Ven and + H.J. Lu from Intel. * Optimized x86-64 trunc and truncf for processors with SSE4.1. diff --git a/sysdeps/x86_64/fpu/multiarch/Makefile b/sysdeps/x86_64/fpu/multiarch/Makefile index 0825340c0c..9a89bfc286 100644 --- a/sysdeps/x86_64/fpu/multiarch/Makefile +++ b/sysdeps/x86_64/fpu/multiarch/Makefile @@ -37,10 +37,10 @@ CFLAGS-slowpow-fma.c = -mfma -mavx2 CFLAGS-s_sin-fma.c = -mfma -mavx2 CFLAGS-s_tan-fma.c = -mfma -mavx2 -libm-sysdep_routines += s_sinf-sse2 s_cosf-sse2 +libm-sysdep_routines += s_sinf-sse2 s_cosf-sse2 s_sincosf-sse2 libm-sysdep_routines += e_exp2f-fma e_expf-fma e_log2f-fma e_logf-fma \ - e_powf-fma s_sinf-fma s_cosf-fma + e_powf-fma s_sinf-fma s_cosf-fma s_sincosf-fma CFLAGS-e_exp2f-fma.c = -mfma -mavx2 CFLAGS-e_expf-fma.c = -mfma -mavx2 @@ -49,6 +49,7 @@ CFLAGS-e_logf-fma.c = -mfma -mavx2 CFLAGS-e_powf-fma.c = -mfma -mavx2 CFLAGS-s_sinf-fma.c = -mfma -mavx2 CFLAGS-s_cosf-fma.c = -mfma -mavx2 +CFLAGS-s_sincosf-fma.c = -mfma -mavx2 libm-sysdep_routines += e_exp-fma4 e_log-fma4 e_pow-fma4 s_atan-fma4 \ e_asin-fma4 e_atan2-fma4 s_sin-fma4 s_tan-fma4 \ diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c b/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c new file mode 100644 index 0000000000..64abe7abca --- /dev/null +++ b/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c @@ -0,0 +1,240 @@ +/* Compute sine and cosine of argument optimized with vector. + Copyright (C) 2017 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 <errno.h> +#include <math.h> +#include <math_private.h> +#include <x86intrin.h> +#include <libm-alias-float.h> +#include "s_sincosf.h" + +#define SINCOSF __sincosf_fma + +#ifndef SINCOSF +# define SINCOSF_FUNC __sincosf +#else +# define SINCOSF_FUNC SINCOSF +#endif + +/* Chebyshev constants for sin and cos, range -PI/4 - PI/4. */ +static const __v2df V0 = { -0x1.5555555551cd9p-3, -0x1.ffffffffe98aep-2}; +static const __v2df V1 = { 0x1.1111110c2688bp-7, 0x1.55555545c50c7p-5 }; +static const __v2df V2 = { -0x1.a019f8b4bd1f9p-13, -0x1.6c16b348b6874p-10 }; +static const __v2df V3 = { 0x1.71d7264e6b5b4p-19, 0x1.a00eb9ac43ccp-16 }; +static const __v2df V4 = { -0x1.a947e1674b58ap-26, -0x1.23c97dd8844d7p-22 }; + +/* Chebyshev constants for sin and cos, range 2^-27 - 2^-5. */ +static const __v2df VC0 = { -0x1.555555543d49dp-3, -0x1.fffffff5cc6fdp-2 }; +static const __v2df VC1 = { 0x1.110f475cec8c5p-7, 0x1.55514b178dac5p-5 }; + +static const __v2df v2ones = { 1.0, 1.0 }; + +/* Compute the sine and cosine values using Chebyshev polynomials where + THETA is the range reduced absolute value of the input + and it is less than Pi/4, + N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide + whether a sine or cosine approximation is more accurate and + SIGNBIT is used to add the correct sign after the Chebyshev + polynomial is computed. */ +static void +reduced_sincos (const double theta, const unsigned int n, + const unsigned int signbit, float *sinx, float *cosx) +{ + __v2df v2x, v2sx, v2cx; + const __v2df v2theta = { theta, theta }; + const __v2df v2theta2 = v2theta * v2theta; + /* Here sinf() and cosf() are calculated using sin Chebyshev polynomial: + x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))). */ + v2x = V3 + v2theta2 * V4; /* S3+x^2*S4. */ + v2x = V2 + v2theta2 * v2x; /* S2+x^2*(S3+x^2*S4). */ + v2x = V1 + v2theta2 * v2x; /* S1+x^2*(S2+x^2*(S3+x^2*S4)). */ + v2x = V0 + v2theta2 * v2x; /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))). */ + v2x = v2theta2 * v2x; + v2cx = v2ones + v2x; + v2sx = v2theta + v2theta * v2x; + /* We are operating on |x|, so we need to add back the original + signbit for sinf. */ + /* Determine positive or negative primary interval. */ + /* Are we in the primary interval of sin or cos? */ + if ((n & 2) == 0) + { + const __v2df v2sign = + { + ones[((n >> 2) & 1) ^ signbit], + ones[((n + 2) >> 2) & 1] + }; + v2cx[0] = v2sx[0]; + v2cx *= v2sign; + __v4sf v4sx = _mm_cvtpd_ps (v2cx); + *sinx = v4sx[0]; + *cosx = v4sx[1]; + } + else + { + const __v2df v2sign = + { + ones[((n + 2) >> 2) & 1], + ones[((n >> 2) & 1) ^ signbit] + }; + v2cx[0] = v2sx[0]; + v2cx *= v2sign; + __v4sf v4sx = _mm_cvtpd_ps (v2cx); + *sinx = v4sx[1]; + *cosx = v4sx[0]; + } +} + +void +SINCOSF_FUNC (float x, float *sinx, float *cosx) +{ + double theta = x; + double abstheta = fabs (theta); + uint32_t ix, xi; + GET_FLOAT_WORD (xi, x); + /* |x| */ + ix = xi & 0x7fffffff; + /* If |x|< Pi/4. */ + if (ix < 0x3f490fdb) + { + if (ix >= 0x3d000000) /* |x| >= 2^-5. */ + { + __v2df v2x, v2sx, v2cx; + const __v2df v2theta = { theta, theta }; + const __v2df v2theta2 = v2theta * v2theta; + /* Chebyshev polynomial of the form for sin and cos. */ + v2x = V3 + v2theta2 * V4; + v2x = V2 + v2theta2 * v2x; + v2x = V1 + v2theta2 * v2x; + v2x = V0 + v2theta2 * v2x; + v2x = v2theta2 * v2x; + v2cx = v2ones + v2x; + v2sx = v2theta + v2theta * v2x; + v2cx[0] = v2sx[0]; + __v4sf v4sx = _mm_cvtpd_ps (v2cx); + *sinx = v4sx[0]; + *cosx = v4sx[1]; + } + else if (ix >= 0x32000000) /* |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). */ + __v2df v2x, v2sx, v2cx; + const __v2df v2theta = { theta, theta }; + const __v2df v2theta2 = v2theta * v2theta; + v2x = VC0 + v2theta * v2theta2 * VC1; + v2x = v2theta2 * v2x; + v2cx = v2ones + v2x; + v2sx = v2theta + v2theta * v2x; + v2cx[0] = v2sx[0]; + __v4sf v4sx = _mm_cvtpd_ps (v2cx); + *sinx = v4sx[0]; + *cosx = v4sx[1]; + } + else + { + /* Handle some special cases. */ + if (ix) + *sinx = theta - (theta * SMALL); + else + *sinx = theta; + *cosx = 1.0 - abstheta; + } + } + else /* |x| >= Pi/4. */ + { + unsigned int signbit = xi >> 31; + if (ix < 0x40e231d6) /* |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]; + reduced_sincos (theta, n, signbit, sinx, cosx); + } + else if (ix < 0x7f800000) + { + if (ix < 0x4b000000) /* |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. */ + reduced_sincos (theta, n, signbit, sinx, cosx); + } + else /* |x| >= 2^23. */ + { + x = fabsf (x); + int exponent + = (ix >> 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; + reduced_sincos (e, l + 1, signbit, sinx, cosx); + } + else + { + e += b; + e += c; + e += d; + if (e <= 1.0) + { + e *= M_PI_4; + reduced_sincos (e, l + 1, signbit, sinx, cosx); + } + else + { + l++; + e -= 2.0; + e *= M_PI_4; + reduced_sincos (e, l + 1, signbit, sinx, cosx); + } + } + } + } + else + { + if (ix == 0x7f800000) + __set_errno (EDOM); + /* sin/cos(Inf or NaN) is NaN. */ + *sinx = *cosx = x - x; + } + } +} + +#ifndef SINCOSF +libm_alias_float (__sincos, sincos) +#endif diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S b/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S new file mode 100644 index 0000000000..51d012bb12 --- /dev/null +++ b/sysdeps/x86_64/fpu/multiarch/s_sincosf-sse2.S @@ -0,0 +1,2 @@ +#define __sincosf __sincosf_sse2 +#include <sysdeps/x86_64/fpu/s_sincosf.S> diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf.c b/sysdeps/x86_64/fpu/multiarch/s_sincosf.c new file mode 100644 index 0000000000..6cb4295558 --- /dev/null +++ b/sysdeps/x86_64/fpu/multiarch/s_sincosf.c @@ -0,0 +1,28 @@ +/* Multiple versions of sincosf. + Copyright (C) 2017 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 <libm-alias-float.h> + +extern void __redirect_sincosf (float, float *, float *); + +#define SYMBOL_NAME sincosf +#include "ifunc-fma.h" + +libc_ifunc_redirected (__redirect_sincosf, __sincosf, IFUNC_SELECTOR ()); + +libm_alias_float (__sincos, sincos) diff --git a/sysdeps/x86_64/fpu/s_sincosf.S b/sysdeps/x86_64/fpu/s_sincosf.S index f00cd0eb3c..2086e8ca5c 100644 --- a/sysdeps/x86_64/fpu/s_sincosf.S +++ b/sysdeps/x86_64/fpu/s_sincosf.S @@ -561,4 +561,6 @@ L(SP_ONE): .type L(SP_ONE), @object ASM_SIZE_DIRECTIVE(L(SP_ONE)) +#ifndef __sincosf libm_alias_float (__sincos, sincos) +#endif |