Message ID | 20181121181102.27119-2-hjl.tools@gmail.com |
---|---|
State | New |
Headers | show |
Series | x86-64: Update sincosf | expand |
On 21/11/2018 16:11, H.J. Lu wrote: > Include <s_sincosf.h> in s_sincosf.c, instead of "s_sincosf.h", to allow > x86-64 s_sincosf.h with vectorized sincosf_poly. Update __sincosf_table > to allow vectorized load in vectorized sincosf_poly. On Broadwell, > bench-sincosf shows: > > Before After Improvement > max 160.273 114.198 40% > min 6.25 5.625 11% > mean 13.0325 10.6462 22% > > Vectorized sincosf_poly shows > > Before After Improvement > max 138.653 114.198 21% > min 5.004 5.625 -11% > mean 11.5934 10.6462 9% > > * sysdeps/ieee754/flt-32/s_sincosf.c: Include <s_sincosf.h> > instead of "s_sincosf.h". > * sysdeps/ieee754/flt-32/s_sincosf.h (sincos_t): Rearranged to > support vectorized load. > (sincosf_poly): Don't define if HAVE_SINCOSF_POLY is defined. > Updated for vectorized load. > (sinf_poly): Updated for vectorized load. > * sysdeps/ieee754/flt-32/s_sincosf_data.c (__sincosf_table): > Rearranged to allow vectorized load. > * sysdeps/x86_64/fpu/s_sincosf.h: New file. > * sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: Just include > <sysdeps/ieee754/flt-32/s_sincosf.c>. > --- > sysdeps/ieee754/flt-32/s_sincosf.c | 2 +- > sysdeps/ieee754/flt-32/s_sincosf.h | 38 ++- > sysdeps/ieee754/flt-32/s_sincosf_data.c | 18 +- > sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c | 271 +------------------ > sysdeps/x86_64/fpu/s_sincosf.h | 57 ++++ > 5 files changed, 93 insertions(+), 293 deletions(-) > create mode 100644 sysdeps/x86_64/fpu/s_sincosf.h > > diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c > index f7e3245097..28dd7530c5 100644 > --- a/sysdeps/ieee754/flt-32/s_sincosf.c > +++ b/sysdeps/ieee754/flt-32/s_sincosf.c > @@ -22,7 +22,7 @@ > #include <math-barriers.h> > #include <libm-alias-float.h> > #include "math_config.h" > -#include "s_sincosf.h" > +#include <s_sincosf.h> > > #ifndef SINCOSF > # define SINCOSF_FUNC __sincosf > diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h > index 1dcb04f235..a1d1639c17 100644 > --- a/sysdeps/ieee754/flt-32/s_sincosf.h > +++ b/sysdeps/ieee754/flt-32/s_sincosf.h > @@ -31,8 +31,24 @@ 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. */ > + /* Cosine polynomial: c0, c1, c2, c3, c4. > + Sine polynomial: s1, s2, s3. */ > + double c0, c1; > + struct > + { > + double s1; > + double c2; > + } s1c2; > + struct > + { > + double s2; > + double c3; > + } s2c3; > + struct > + { > + double s3; > + double c4; > + } s3c4; > } sincos_t; > I don't think this should be a problem for other architectures, do you see any possible issue about changing the layout? > /* Polynomial data (the cosine polynomial is negated in the 2nd entry). */ > @@ -48,6 +64,7 @@ abstop12 (float x) > return (asuint (x) >> 20) & 0x7ff; > } > > +#ifndef HAVE_SINCOSF_POLY > /* 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. */ > @@ -59,8 +76,8 @@ sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp, > > x4 = x2 * x2; > x3 = x2 * x; > - c2 = p->c3 + x2 * p->c4; > - s1 = p->s2 + x2 * p->s3; > + c2 = p->s2c3.c3 + x2 * p->s3c4.c4; > + s1 = p->s2c3.s2 + x2 * p->s3c4.s3; > > /* Swap sin/cos result based on quadrant. */ > float *tmp = (n & 1 ? cosp : sinp); > @@ -71,12 +88,13 @@ sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp, > x5 = x3 * x2; > x6 = x4 * x2; > > - s = x + x3 * p->s1; > - c = c1 + x4 * p->c2; > + s = x + x3 * p->s1c2.s1; > + c = c1 + x4 * p->s1c2.c2; > > *sinp = s + x5 * s1; > *cosp = c + x6 * c2; > } > +#endif I think current trend it to avoid such construction based on define macros and instead use file-based definition which can be override. For this specific change, wouldn't be better to create a generic s_sincosf_poly.h with sincosf_poly and override it on x86_64 folder? Maybe also move s_sincosf_t to its own header, so s_sincosf_poly.h can include it as well or an architecture can use a different layout if it requires so.
On Fri, Nov 30, 2018 at 12:40 PM Adhemerval Zanella <adhemerval.zanella@linaro.org> wrote: > > > > On 21/11/2018 16:11, H.J. Lu wrote: > > Include <s_sincosf.h> in s_sincosf.c, instead of "s_sincosf.h", to allow > > x86-64 s_sincosf.h with vectorized sincosf_poly. Update __sincosf_table > > to allow vectorized load in vectorized sincosf_poly. On Broadwell, > > bench-sincosf shows: > > > > Before After Improvement > > max 160.273 114.198 40% > > min 6.25 5.625 11% > > mean 13.0325 10.6462 22% > > > > Vectorized sincosf_poly shows > > > > Before After Improvement > > max 138.653 114.198 21% > > min 5.004 5.625 -11% > > mean 11.5934 10.6462 9% > > > > * sysdeps/ieee754/flt-32/s_sincosf.c: Include <s_sincosf.h> > > instead of "s_sincosf.h". > > * sysdeps/ieee754/flt-32/s_sincosf.h (sincos_t): Rearranged to > > support vectorized load. > > (sincosf_poly): Don't define if HAVE_SINCOSF_POLY is defined. > > Updated for vectorized load. > > (sinf_poly): Updated for vectorized load. > > * sysdeps/ieee754/flt-32/s_sincosf_data.c (__sincosf_table): > > Rearranged to allow vectorized load. > > * sysdeps/x86_64/fpu/s_sincosf.h: New file. > > * sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: Just include > > <sysdeps/ieee754/flt-32/s_sincosf.c>. > > --- > > sysdeps/ieee754/flt-32/s_sincosf.c | 2 +- > > sysdeps/ieee754/flt-32/s_sincosf.h | 38 ++- > > sysdeps/ieee754/flt-32/s_sincosf_data.c | 18 +- > > sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c | 271 +------------------ > > sysdeps/x86_64/fpu/s_sincosf.h | 57 ++++ > > 5 files changed, 93 insertions(+), 293 deletions(-) > > create mode 100644 sysdeps/x86_64/fpu/s_sincosf.h > > > > diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c > > index f7e3245097..28dd7530c5 100644 > > --- a/sysdeps/ieee754/flt-32/s_sincosf.c > > +++ b/sysdeps/ieee754/flt-32/s_sincosf.c > > @@ -22,7 +22,7 @@ > > #include <math-barriers.h> > > #include <libm-alias-float.h> > > #include "math_config.h" > > -#include "s_sincosf.h" > > +#include <s_sincosf.h> > > > > #ifndef SINCOSF > > # define SINCOSF_FUNC __sincosf > > diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h > > index 1dcb04f235..a1d1639c17 100644 > > --- a/sysdeps/ieee754/flt-32/s_sincosf.h > > +++ b/sysdeps/ieee754/flt-32/s_sincosf.h > > @@ -31,8 +31,24 @@ 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. */ > > + /* Cosine polynomial: c0, c1, c2, c3, c4. > > + Sine polynomial: s1, s2, s3. */ > > + double c0, c1; > > + struct > > + { > > + double s1; > > + double c2; > > + } s1c2; > > + struct > > + { > > + double s2; > > + double c3; > > + } s2c3; > > + struct > > + { > > + double s3; > > + double c4; > > + } s3c4; > > } sincos_t; > > > > I don't think this should be a problem for other architectures, do you see > any possible issue about changing the layout? Not I can think of. > > /* Polynomial data (the cosine polynomial is negated in the 2nd entry). */ > > @@ -48,6 +64,7 @@ abstop12 (float x) > > return (asuint (x) >> 20) & 0x7ff; > > } > > > > +#ifndef HAVE_SINCOSF_POLY > > /* 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. */ > > @@ -59,8 +76,8 @@ sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp, > > > > x4 = x2 * x2; > > x3 = x2 * x; > > - c2 = p->c3 + x2 * p->c4; > > - s1 = p->s2 + x2 * p->s3; > > + c2 = p->s2c3.c3 + x2 * p->s3c4.c4; > > + s1 = p->s2c3.s2 + x2 * p->s3c4.s3; > > > > /* Swap sin/cos result based on quadrant. */ > > float *tmp = (n & 1 ? cosp : sinp); > > @@ -71,12 +88,13 @@ sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp, > > x5 = x3 * x2; > > x6 = x4 * x2; > > > > - s = x + x3 * p->s1; > > - c = c1 + x4 * p->c2; > > + s = x + x3 * p->s1c2.s1; > > + c = c1 + x4 * p->s1c2.c2; > > > > *sinp = s + x5 * s1; > > *cosp = c + x6 * c2; > > } > > +#endif > > I think current trend it to avoid such construction based on define macros > and instead use file-based definition which can be override. For this specific > change, wouldn't be better to create a generic s_sincosf_poly.h with Like this? > sincosf_poly and override it on x86_64 folder? Maybe also move s_sincosf_t > to its own header, so s_sincosf_poly.h can include it as well or an > architecture can use a different layout if it requires so. I don't think that is necessary.
On 30/11/18 21:31, H.J. Lu wrote: > On Fri, Nov 30, 2018 at 12:40 PM Adhemerval Zanella > <adhemerval.zanella@linaro.org> wrote: >> On 21/11/2018 16:11, H.J. Lu wrote: >>> Include <s_sincosf.h> in s_sincosf.c, instead of "s_sincosf.h", to allow >>> x86-64 s_sincosf.h with vectorized sincosf_poly. Update __sincosf_table ... >>> --- a/sysdeps/ieee754/flt-32/s_sincosf.h >>> +++ b/sysdeps/ieee754/flt-32/s_sincosf.h >>> @@ -31,8 +31,24 @@ 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. */ >>> + /* Cosine polynomial: c0, c1, c2, c3, c4. >>> + Sine polynomial: s1, s2, s3. */ >>> + double c0, c1; >>> + struct >>> + { >>> + double s1; >>> + double c2; >>> + } s1c2; >>> + struct >>> + { >>> + double s2; >>> + double c3; >>> + } s2c3; >>> + struct >>> + { >>> + double s3; >>> + double c4; >>> + } s3c4; >>> } sincos_t; >>> >> >> I don't think this should be a problem for other architectures, do you see >> any possible issue about changing the layout? > > Not I can think of. if you change the layout it's likely a measurable regression on aarch64 targets because it was optimized such that ldp (load pair) can be used when loading coeffs (although this should not have a huge effect). i think it should be possible to override the header and data files with target specific definitions (or use ifdefs based on some target macros) so we can have both layouts.
On 03/12/2018 08:47, Szabolcs Nagy wrote: > On 30/11/18 21:31, H.J. Lu wrote: >> On Fri, Nov 30, 2018 at 12:40 PM Adhemerval Zanella >> <adhemerval.zanella@linaro.org> wrote: >>> On 21/11/2018 16:11, H.J. Lu wrote: >>>> Include <s_sincosf.h> in s_sincosf.c, instead of "s_sincosf.h", to allow >>>> x86-64 s_sincosf.h with vectorized sincosf_poly. Update __sincosf_table > ... >>>> --- a/sysdeps/ieee754/flt-32/s_sincosf.h >>>> +++ b/sysdeps/ieee754/flt-32/s_sincosf.h >>>> @@ -31,8 +31,24 @@ 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. */ >>>> + /* Cosine polynomial: c0, c1, c2, c3, c4. >>>> + Sine polynomial: s1, s2, s3. */ >>>> + double c0, c1; >>>> + struct >>>> + { >>>> + double s1; >>>> + double c2; >>>> + } s1c2; >>>> + struct >>>> + { >>>> + double s2; >>>> + double c3; >>>> + } s2c3; >>>> + struct >>>> + { >>>> + double s3; >>>> + double c4; >>>> + } s3c4; >>>> } sincos_t; >>>> >>> >>> I don't think this should be a problem for other architectures, do you see >>> any possible issue about changing the layout? >> >> Not I can think of. > > if you change the layout it's likely a measurable > regression on aarch64 targets because it was optimized > such that ldp (load pair) can be used when loading > coeffs (although this should not have a huge effect). > > i think it should be possible to override the header > and data files with target specific definitions (or > use ifdefs based on some target macros) so we can > have both layouts. I did check on a A53 I saw no regressions with benchtests. Do you see any regressions on other chips or systems? If it is the case one option could be use my suggestion to move s_sincosf_t to its own header.
diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c index f7e3245097..28dd7530c5 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.c +++ b/sysdeps/ieee754/flt-32/s_sincosf.c @@ -22,7 +22,7 @@ #include <math-barriers.h> #include <libm-alias-float.h> #include "math_config.h" -#include "s_sincosf.h" +#include <s_sincosf.h> #ifndef SINCOSF # define SINCOSF_FUNC __sincosf diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h index 1dcb04f235..a1d1639c17 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf.h +++ b/sysdeps/ieee754/flt-32/s_sincosf.h @@ -31,8 +31,24 @@ 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. */ + /* Cosine polynomial: c0, c1, c2, c3, c4. + Sine polynomial: s1, s2, s3. */ + double c0, c1; + struct + { + double s1; + double c2; + } s1c2; + struct + { + double s2; + double c3; + } s2c3; + struct + { + double s3; + double c4; + } s3c4; } sincos_t; /* Polynomial data (the cosine polynomial is negated in the 2nd entry). */ @@ -48,6 +64,7 @@ abstop12 (float x) return (asuint (x) >> 20) & 0x7ff; } +#ifndef HAVE_SINCOSF_POLY /* 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. */ @@ -59,8 +76,8 @@ sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp, x4 = x2 * x2; x3 = x2 * x; - c2 = p->c3 + x2 * p->c4; - s1 = p->s2 + x2 * p->s3; + c2 = p->s2c3.c3 + x2 * p->s3c4.c4; + s1 = p->s2c3.s2 + x2 * p->s3c4.s3; /* Swap sin/cos result based on quadrant. */ float *tmp = (n & 1 ? cosp : sinp); @@ -71,12 +88,13 @@ sincosf_poly (double x, double x2, const sincos_t *p, int n, float *sinp, x5 = x3 * x2; x6 = x4 * x2; - s = x + x3 * p->s1; - c = c1 + x4 * p->c2; + s = x + x3 * p->s1c2.s1; + c = c1 + x4 * p->s1c2.c2; *sinp = s + x5 * s1; *cosp = c + x6 * c2; } +#endif /* 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. */ @@ -88,21 +106,21 @@ sinf_poly (double x, double x2, const sincos_t *p, int n) if ((n & 1) == 0) { x3 = x * x2; - s1 = p->s2 + x2 * p->s3; + s1 = p->s2c3.s2 + x2 * p->s3c4.s3; x7 = x3 * x2; - s = x + x3 * p->s1; + s = x + x3 * p->s1c2.s1; return s + x7 * s1; } else { x4 = x2 * x2; - c2 = p->c3 + x2 * p->c4; + c2 = p->s2c3.c3 + x2 * p->s3c4.c4; c1 = p->c0 + x2 * p->c1; x6 = x4 * x2; - c = c1 + x4 * p->c2; + c = c1 + x4 * p->s1c2.c2; return c + x6 * c2; } diff --git a/sysdeps/ieee754/flt-32/s_sincosf_data.c b/sysdeps/ieee754/flt-32/s_sincosf_data.c index 21fc2b60f9..6d3ed77d4c 100644 --- a/sysdeps/ieee754/flt-32/s_sincosf_data.c +++ b/sysdeps/ieee754/flt-32/s_sincosf_data.c @@ -35,12 +35,9 @@ const sincos_t __sincosf_table[2] = 0x1.921FB54442D18p0, 0x1p0, -0x1.ffffffd0c621cp-2, - 0x1.55553e1068f19p-5, - -0x1.6c087e89a359dp-10, - 0x1.99343027bf8c3p-16, - -0x1.555545995a603p-3, - 0x1.1107605230bc4p-7, - -0x1.994eb3774cf24p-13 + { -0x1.555545995a603p-3, 0x1.55553e1068f19p-5 }, + { 0x1.1107605230bc4p-7, -0x1.6c087e89a359dp-10 }, + { -0x1.994eb3774cf24p-13, 0x1.99343027bf8c3p-16 } }, { { 1.0, -1.0, -1.0, 1.0 }, @@ -52,12 +49,9 @@ const sincos_t __sincosf_table[2] = 0x1.921FB54442D18p0, -0x1p0, 0x1.ffffffd0c621cp-2, - -0x1.55553e1068f19p-5, - 0x1.6c087e89a359dp-10, - -0x1.99343027bf8c3p-16, - -0x1.555545995a603p-3, - 0x1.1107605230bc4p-7, - -0x1.994eb3774cf24p-13 + { -0x1.555545995a603p-3, -0x1.55553e1068f19p-5 }, + { 0x1.1107605230bc4p-7, 0x1.6c087e89a359dp-10 }, + { -0x1.994eb3774cf24p-13, -0x1.99343027bf8c3p-16 } } }; diff --git a/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c b/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c index 0b80c4fe0d..253dab15d8 100644 --- a/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c +++ b/sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c @@ -1,271 +1,2 @@ -/* 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> - #define SINCOSF __sincosf_fma - -#ifndef SINCOSF -# define SINCOSF_FUNC __sincosf -#else -# define SINCOSF_FUNC SINCOSF -#endif - -/* PI/2 with 98 bits of accuracy. */ -static const double PI_2_hi = 0x1.921fb544p+0; -static const double PI_2_lo = 0x1.0b4611a626332p-34; - -static const double SMALL = 0x1p-50; /* 2^-50. */ -static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI. */ - -#define FLOAT_EXPONENT_SHIFT 23 -#define FLOAT_EXPONENT_BIAS 127 - -static const double pio2_table[] = { - 0 * M_PI_2, - 1 * M_PI_2, - 2 * M_PI_2, - 3 * M_PI_2, - 4 * M_PI_2, - 5 * M_PI_2 -}; - -static const double invpio4_table[] = { - 0x0p+0, - 0x1.45f306cp+0, - 0x1.c9c882ap-28, - 0x1.4fe13a8p-58, - 0x1.f47d4dp-85, - 0x1.bb81b6cp-112, - 0x1.4acc9ep-142, - 0x1.0e4107cp-169 -}; - -static const double ones[] = { 1.0, -1.0 }; - -/* 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 +#include <sysdeps/ieee754/flt-32/s_sincosf.c> diff --git a/sysdeps/x86_64/fpu/s_sincosf.h b/sysdeps/x86_64/fpu/s_sincosf.h new file mode 100644 index 0000000000..ecca29db12 --- /dev/null +++ b/sysdeps/x86_64/fpu/s_sincosf.h @@ -0,0 +1,57 @@ +/* x86-64 sincosf_poly for sincosf. + 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/>. */ + +#define HAVE_SINCOSF_POLY +#include_next <s_sincosf.h> +#include <x86intrin.h> + +/* 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) +{ + __v2df vx2x2 = { x2, x2 }; + __v2df vxx2 = { x, x2 }; + __v2df vps1c2 = (__v2df) _mm_loadu_pd (&p->s1c2.s1); + __v2df vps2c3 = (__v2df) _mm_loadu_pd (&p->s2c3.s2); + __v2df vps3c4 = (__v2df) _mm_loadu_pd (&p->s3c4.s3); + __v2df vx3x4, vs1c2; + + vx3x4 = vx2x2 * vxx2; + vs1c2 = vps2c3 + vx2x2 * vps3c4; + + /* Swap sin/cos result based on quadrant. */ + if (n & 1) + { + float *tmp = cosp; + cosp = sinp; + sinp = tmp; + } + + double c1 = p->c0 + x2 * p->c1; + __v2df vxc1 = { x, c1 }; + __v2df vx5x6 = vx3x4 * vx2x2; + + __v2df vsincos = vxc1 + vx3x4 * vps1c2; + vsincos = vsincos + vx5x6 * vs1c2; + __v4sf v4sf = _mm_cvtpd_ps (vsincos); + *sinp = v4sf[0]; + *cosp = v4sf[1]; +}