* Re: [PATCH 1/3] Use generic vector computations in s_sincosf.h
@ 2018-12-14 19:44 Wilco Dijkstra
2018-12-14 20:55 ` Carlos O'Donell
0 siblings, 1 reply; 6+ messages in thread
From: Wilco Dijkstra @ 2018-12-14 19:44 UTC (permalink / raw)
To: 'GNU C Library', H.J. Lu; +Cc: nd
Hi,
I finally had a chance to benchmark this on a recent core using the traces posted
here: https://www.sourceware.org/ml/libc-alpha/2018-12/msg00492.html
I used the same trace for sincosf but this requires a small wrapper to accurately
measure latency and throughput - I'll post the patch for this next week.
On Sandy Bridge I get 10% gain in throughput with the vector FMA version of sincosf
(at the cost of ~2% extra latency), but unfortunately it's 5% slower on AArch64
and latency increases by 10% due to GCC using some unnecessary DUPs and
lane inserts. Sinf and cosf are also affected by the data layout change, throughput
and latency are 1.7% and 1.0% worse respectively.
So it seems the best approach would be to make the vector version conditional on
a macro defined by the target so each target can choose the fastest variant.
Wilco
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [PATCH 1/3] Use generic vector computations in s_sincosf.h
2018-12-14 19:44 [PATCH 1/3] Use generic vector computations in s_sincosf.h Wilco Dijkstra
@ 2018-12-14 20:55 ` Carlos O'Donell
2018-12-17 0:16 ` [PATCH] x86: " H.J. Lu
2018-12-17 16:02 ` [PATCH 1/3] " Wilco Dijkstra
0 siblings, 2 replies; 6+ messages in thread
From: Carlos O'Donell @ 2018-12-14 20:55 UTC (permalink / raw)
To: Wilco Dijkstra, 'GNU C Library', H.J. Lu; +Cc: nd
On 12/14/18 2:39 PM, Wilco Dijkstra wrote:
> Hi,
>
> I finally had a chance to benchmark this on a recent core using the traces posted
> here: https://www.sourceware.org/ml/libc-alpha/2018-12/msg00492.html
> I used the same trace for sincosf but this requires a small wrapper to accurately
> measure latency and throughput - I'll post the patch for this next week.
>
> On Sandy Bridge I get 10% gain in throughput with the vector FMA version of sincosf
> (at the cost of ~2% extra latency), but unfortunately it's 5% slower on AArch64
> and latency increases by 10% due to GCC using some unnecessary DUPs and
> lane inserts. Sinf and cosf are also affected by the data layout change, throughput
> and latency are 1.7% and 1.0% worse respectively.
>
> So it seems the best approach would be to make the vector version conditional on
> a macro defined by the target so each target can choose the fastest variant.
In general I agree.
I wonder what is the cost of using a tunable with defaults set differently by
architecture?
This way if a architectural licensee does something different for AArch64 that
changes this back the other way they could tune it differently with a global
env var setting?
--
Cheers,
Carlos.
^ permalink raw reply [flat|nested] 6+ messages in thread
* [PATCH] x86: Use generic vector computations in s_sincosf.h
2018-12-14 20:55 ` Carlos O'Donell
@ 2018-12-17 0:16 ` H.J. Lu
2018-12-26 9:53 ` Siddhesh Poyarekar
2018-12-17 16:02 ` [PATCH 1/3] " Wilco Dijkstra
1 sibling, 1 reply; 6+ messages in thread
From: H.J. Lu @ 2018-12-17 0:16 UTC (permalink / raw)
To: Carlos O'Donell; +Cc: Wilco Dijkstra, 'GNU C Library', nd
On Fri, Dec 14, 2018 at 03:36:00PM -0500, Carlos O'Donell wrote:
> On 12/14/18 2:39 PM, Wilco Dijkstra wrote:
> > Hi,
> >
> > I finally had a chance to benchmark this on a recent core using the traces posted
> > here: https://www.sourceware.org/ml/libc-alpha/2018-12/msg00492.html
> > I used the same trace for sincosf but this requires a small wrapper to accurately
> > measure latency and throughput - I'll post the patch for this next week.
> >
> > On Sandy Bridge I get 10% gain in throughput with the vector FMA version of sincosf
> > (at the cost of ~2% extra latency), but unfortunately it's 5% slower on AArch64
> > and latency increases by 10% due to GCC using some unnecessary DUPs and
> > lane inserts. Sinf and cosf are also affected by the data layout change, throughput
> > and latency are 1.7% and 1.0% worse respectively.
> >
> > So it seems the best approach would be to make the vector version conditional on
> > a macro defined by the target so each target can choose the fastest variant.
>
> In general I agree.
>
> I wonder what is the cost of using a tunable with defaults set differently by
> architecture?
>
> This way if a architectural licensee does something different for AArch64 that
> changes this back the other way they could tune it differently with a global
> env var setting?
>
Here is the updated patch to use generic vector computations only for
x86. Tested on i686 and x86-64. OK for master branch?
Thanks.
H.J.
--
Use generic vector computations in s_sincosf.h to support vectorized
s_sincosf. Update __sincosf_table for vectorized s_sincosf. 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%
Tested on x86-64 and i686 as well as with build-many-glibcs.py.
* sysdeps/ieee754/flt-32/s_sincosf.h: Include <sincosf_poly.h>.
(sincos_t, sincosf_poly, sinf_poly): Moved to ...
* sysdeps/ieee754/flt-32/sincosf_poly.h: Here. New file.
* sysdeps/x86/fpu/s_sincosf_data.c: New file.
* sysdeps/x86/fpu/sincosf_poly.h: Likewise.
* sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: Just include
<sysdeps/ieee754/flt-32/s_sincosf.c>.
---
sysdeps/ieee754/flt-32/s_sincosf.h | 71 +----
sysdeps/ieee754/flt-32/sincosf_poly.h | 87 ++++++
sysdeps/x86/fpu/s_sincosf_data.c | 68 +++++
sysdeps/x86/fpu/sincosf_poly.h | 111 ++++++++
sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c | 271 +------------------
5 files changed, 268 insertions(+), 340 deletions(-)
create mode 100644 sysdeps/ieee754/flt-32/sincosf_poly.h
create mode 100644 sysdeps/x86/fpu/s_sincosf_data.c
create mode 100644 sysdeps/x86/fpu/sincosf_poly.h
diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h
index 1dcb04f235..54a9da543b 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf.h
+++ b/sysdeps/ieee754/flt-32/s_sincosf.h
@@ -19,22 +19,13 @@
#include <stdint.h>
#include <math.h>
#include "math_config.h"
+#include <sincosf_poly.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] attribute_hidden;
@@ -48,66 +39,6 @@ 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
diff --git a/sysdeps/ieee754/flt-32/sincosf_poly.h b/sysdeps/ieee754/flt-32/sincosf_poly.h
new file mode 100644
index 0000000000..f49c75ad78
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/sincosf_poly.h
@@ -0,0 +1,87 @@
+/* Used by sinf, cosf and sincosf functions.
+ 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/>. */
+
+/* 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;
+
+/* 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;
+ }
+}
diff --git a/sysdeps/x86/fpu/s_sincosf_data.c b/sysdeps/x86/fpu/s_sincosf_data.c
new file mode 100644
index 0000000000..b219644884
--- /dev/null
+++ b/sysdeps/x86/fpu/s_sincosf_data.c
@@ -0,0 +1,68 @@
+/* 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 <sysdeps/ieee754/flt-32/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.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 },
+#if TOINT_INTRINSICS
+ 0x1.45F306DC9C883p-1,
+#else
+ 0x1.45F306DC9C883p+23,
+#endif
+ 0x1.921FB54442D18p0,
+ -0x1p0,
+ 0x1.ffffffd0c621cp-2,
+ { -0x1.555545995a603p-3, -0x1.55553e1068f19p-5 },
+ { 0x1.1107605230bc4p-7, 0x1.6c087e89a359dp-10 },
+ { -0x1.994eb3774cf24p-13, -0x1.99343027bf8c3p-16 }
+ }
+};
+
+/* 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/sysdeps/x86/fpu/sincosf_poly.h b/sysdeps/x86/fpu/sincosf_poly.h
new file mode 100644
index 0000000000..396ba8993d
--- /dev/null
+++ b/sysdeps/x86/fpu/sincosf_poly.h
@@ -0,0 +1,111 @@
+/* Used by sinf, cosf and sincosf functions. X86-64 version.
+ 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/>. */
+
+typedef double v2df_t __attribute__ ((vector_size (2 * sizeof (double))));
+
+#ifdef __SSE2_MATH__
+typedef float v4sf_t __attribute__ ((vector_size (4 * sizeof (float))));
+
+static inline void
+v2df_to_sf (v2df_t v2df, float *f0p, float *f1p)
+{
+ v4sf_t v4sf = __builtin_ia32_cvtpd2ps (v2df);
+ *f0p = v4sf[0];
+ *f1p = v4sf[1];
+}
+#else
+static inline void
+v2df_to_sf (v2df_t v2df, float *f0p, float *f1p)
+{
+ *f0p = (float) v2df[0];
+ *f1p = (float) v2df[1];
+}
+#endif
+
+/* 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. */
+ /* Cosine polynomial: c0, c1, c2, c3, c4.
+ Sine polynomial: s1, s2, s3. */
+ double c0, c1;
+ v2df_t s1c2, s2c3, s3c4;
+} sincos_t;
+
+/* 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_t vx2x2 = { x2, x2 };
+ v2df_t vxx2 = { x, x2 };
+ v2df_t vx3x4, vs1c2;
+
+ vx3x4 = vx2x2 * vxx2;
+ vs1c2 = p->s2c3 + vx2x2 * p->s3c4;
+
+ /* 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_t vxc1 = { x, c1 };
+ v2df_t vx5x6 = vx3x4 * vx2x2;
+
+ v2df_t vsincos = vxc1 + vx3x4 * p->s1c2;
+ vsincos = vsincos + vx5x6 * vs1c2;
+ v2df_to_sf (vsincos, sinp, cosp);
+}
+
+/* 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->s2c3[0] + x2 * p->s3c4[0];
+
+ x7 = x3 * x2;
+ s = x + x3 * p->s1c2[0];
+
+ return s + x7 * s1;
+ }
+ else
+ {
+ x4 = x2 * x2;
+ c2 = p->s2c3[1] + x2 * p->s3c4[1];
+ c1 = p->c0 + x2 * p->c1;
+
+ x6 = x4 * x2;
+ c = c1 + x4 * p->s1c2[1];
+
+ return c + x6 * c2;
+ }
+}
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>
--
2.19.2
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [PATCH 1/3] Use generic vector computations in s_sincosf.h
2018-12-14 20:55 ` Carlos O'Donell
2018-12-17 0:16 ` [PATCH] x86: " H.J. Lu
@ 2018-12-17 16:02 ` Wilco Dijkstra
1 sibling, 0 replies; 6+ messages in thread
From: Wilco Dijkstra @ 2018-12-17 16:02 UTC (permalink / raw)
To: Carlos O'Donell, 'GNU C Library', H.J. Lu; +Cc: nd
Hi Carlos,
> I wonder what is the cost of using a tunable with defaults set differently by
> architecture?
>
> This way if a architectural licensee does something different for AArch64 that
> changes this back the other way they could tune it differently with a global
> env var setting?
I'm not sure how it would be useful as a tunable - the overhead is due to
extra scalar/vector conversion instructions that GCC inserts, so it is
not dependent on the micro architecture (and you'd use an ifunc if it was).
We can revisit this when GCC's handling of vectors has been improved.
However I think it's more interesting to go after the big performance gains
first, for example we can get a 50% speedup with an improved interface:
https://sourceware.org/ml/libc-alpha/2018-09/msg00185.html
Wilco
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [PATCH] x86: Use generic vector computations in s_sincosf.h
2018-12-17 0:16 ` [PATCH] x86: " H.J. Lu
@ 2018-12-26 9:53 ` Siddhesh Poyarekar
2018-12-26 18:56 ` H.J. Lu
0 siblings, 1 reply; 6+ messages in thread
From: Siddhesh Poyarekar @ 2018-12-26 9:53 UTC (permalink / raw)
To: H.J. Lu, Carlos O'Donell; +Cc: Wilco Dijkstra, 'GNU C Library', nd
On 17/12/18 4:27 AM, H.J. Lu wrote:
> Here is the updated patch to use generic vector computations only for
> x86. Tested on i686 and x86-64. OK for master branch?
I assume you'll add the original description to this so that there's a
clearer git commit log.
>
> Thanks.
>
> H.J.
> --
> Use generic vector computations in s_sincosf.h to support vectorized
> s_sincosf. Update __sincosf_table for vectorized s_sincosf. 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%
I assume Wilco's performance gain also remains with this patch given
that the crux of the code hasn't changed.
Looks OK to me.
Siddhesh
^ permalink raw reply [flat|nested] 6+ messages in thread
* Re: [PATCH] x86: Use generic vector computations in s_sincosf.h
2018-12-26 9:53 ` Siddhesh Poyarekar
@ 2018-12-26 18:56 ` H.J. Lu
0 siblings, 0 replies; 6+ messages in thread
From: H.J. Lu @ 2018-12-26 18:56 UTC (permalink / raw)
To: Siddhesh Poyarekar
Cc: Carlos O'Donell, Wilco Dijkstra, 'GNU C Library', nd
On Wed, Dec 26, 2018 at 03:21:14PM +0530, Siddhesh Poyarekar wrote:
> On 17/12/18 4:27 AM, H.J. Lu wrote:
> > Here is the updated patch to use generic vector computations only for
> > x86. Tested on i686 and x86-64. OK for master branch?
>
> I assume you'll add the original description to this so that there's a
> clearer git commit log.
>
> >
> > Thanks.
> >
> > H.J.
> > --
> > Use generic vector computations in s_sincosf.h to support vectorized
> > s_sincosf. Update __sincosf_table for vectorized s_sincosf. 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%
>
> I assume Wilco's performance gain also remains with this patch given that
> the crux of the code hasn't changed.
>
> Looks OK to me.
>
This is the patch I am checking in.
Thanks.
H.J.
----
Add <sincosf_poly.h> and include it in s_sincosf.h to allow vectorized
sincosf_poly. Add x86 sincosf_poly.h to vectorize 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%
Tested on x86-64 and i686 as well as with build-many-glibcs.py.
* sysdeps/ieee754/flt-32/s_sincosf.h: Include <sincosf_poly.h>.
(sincos_t, sincosf_poly, sinf_poly): Moved to ...
* sysdeps/ieee754/flt-32/sincosf_poly.h: Here. New file.
* sysdeps/x86/fpu/s_sincosf_data.c: New file.
* sysdeps/x86/fpu/sincosf_poly.h: Likewise.
* sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c: Just include
<sysdeps/ieee754/flt-32/s_sincosf.c>.
---
sysdeps/ieee754/flt-32/s_sincosf.h | 71 +----
sysdeps/ieee754/flt-32/sincosf_poly.h | 87 ++++++
sysdeps/x86/fpu/s_sincosf_data.c | 68 +++++
sysdeps/x86/fpu/sincosf_poly.h | 111 ++++++++
sysdeps/x86_64/fpu/multiarch/s_sincosf-fma.c | 271 +------------------
5 files changed, 268 insertions(+), 340 deletions(-)
create mode 100644 sysdeps/ieee754/flt-32/sincosf_poly.h
create mode 100644 sysdeps/x86/fpu/s_sincosf_data.c
create mode 100644 sysdeps/x86/fpu/sincosf_poly.h
diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h
index 1dcb04f235..54a9da543b 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf.h
+++ b/sysdeps/ieee754/flt-32/s_sincosf.h
@@ -19,22 +19,13 @@
#include <stdint.h>
#include <math.h>
#include "math_config.h"
+#include <sincosf_poly.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] attribute_hidden;
@@ -48,66 +39,6 @@ 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
diff --git a/sysdeps/ieee754/flt-32/sincosf_poly.h b/sysdeps/ieee754/flt-32/sincosf_poly.h
new file mode 100644
index 0000000000..f49c75ad78
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/sincosf_poly.h
@@ -0,0 +1,87 @@
+/* Used by sinf, cosf and sincosf functions.
+ 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/>. */
+
+/* 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;
+
+/* 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;
+ }
+}
diff --git a/sysdeps/x86/fpu/s_sincosf_data.c b/sysdeps/x86/fpu/s_sincosf_data.c
new file mode 100644
index 0000000000..b219644884
--- /dev/null
+++ b/sysdeps/x86/fpu/s_sincosf_data.c
@@ -0,0 +1,68 @@
+/* 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 <sysdeps/ieee754/flt-32/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.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 },
+#if TOINT_INTRINSICS
+ 0x1.45F306DC9C883p-1,
+#else
+ 0x1.45F306DC9C883p+23,
+#endif
+ 0x1.921FB54442D18p0,
+ -0x1p0,
+ 0x1.ffffffd0c621cp-2,
+ { -0x1.555545995a603p-3, -0x1.55553e1068f19p-5 },
+ { 0x1.1107605230bc4p-7, 0x1.6c087e89a359dp-10 },
+ { -0x1.994eb3774cf24p-13, -0x1.99343027bf8c3p-16 }
+ }
+};
+
+/* 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/sysdeps/x86/fpu/sincosf_poly.h b/sysdeps/x86/fpu/sincosf_poly.h
new file mode 100644
index 0000000000..396ba8993d
--- /dev/null
+++ b/sysdeps/x86/fpu/sincosf_poly.h
@@ -0,0 +1,111 @@
+/* Used by sinf, cosf and sincosf functions. X86-64 version.
+ 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/>. */
+
+typedef double v2df_t __attribute__ ((vector_size (2 * sizeof (double))));
+
+#ifdef __SSE2_MATH__
+typedef float v4sf_t __attribute__ ((vector_size (4 * sizeof (float))));
+
+static inline void
+v2df_to_sf (v2df_t v2df, float *f0p, float *f1p)
+{
+ v4sf_t v4sf = __builtin_ia32_cvtpd2ps (v2df);
+ *f0p = v4sf[0];
+ *f1p = v4sf[1];
+}
+#else
+static inline void
+v2df_to_sf (v2df_t v2df, float *f0p, float *f1p)
+{
+ *f0p = (float) v2df[0];
+ *f1p = (float) v2df[1];
+}
+#endif
+
+/* 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. */
+ /* Cosine polynomial: c0, c1, c2, c3, c4.
+ Sine polynomial: s1, s2, s3. */
+ double c0, c1;
+ v2df_t s1c2, s2c3, s3c4;
+} sincos_t;
+
+/* 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_t vx2x2 = { x2, x2 };
+ v2df_t vxx2 = { x, x2 };
+ v2df_t vx3x4, vs1c2;
+
+ vx3x4 = vx2x2 * vxx2;
+ vs1c2 = p->s2c3 + vx2x2 * p->s3c4;
+
+ /* 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_t vxc1 = { x, c1 };
+ v2df_t vx5x6 = vx3x4 * vx2x2;
+
+ v2df_t vsincos = vxc1 + vx3x4 * p->s1c2;
+ vsincos = vsincos + vx5x6 * vs1c2;
+ v2df_to_sf (vsincos, sinp, cosp);
+}
+
+/* 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->s2c3[0] + x2 * p->s3c4[0];
+
+ x7 = x3 * x2;
+ s = x + x3 * p->s1c2[0];
+
+ return s + x7 * s1;
+ }
+ else
+ {
+ x4 = x2 * x2;
+ c2 = p->s2c3[1] + x2 * p->s3c4[1];
+ c1 = p->c0 + x2 * p->c1;
+
+ x6 = x4 * x2;
+ c = c1 + x4 * p->s1c2[1];
+
+ return c + x6 * c2;
+ }
+}
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>
--
2.20.1
^ permalink raw reply [flat|nested] 6+ messages in thread
end of thread, other threads:[~2018-12-26 14:55 UTC | newest]
Thread overview: 6+ messages (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2018-12-14 19:44 [PATCH 1/3] Use generic vector computations in s_sincosf.h Wilco Dijkstra
2018-12-14 20:55 ` Carlos O'Donell
2018-12-17 0:16 ` [PATCH] x86: " H.J. Lu
2018-12-26 9:53 ` Siddhesh Poyarekar
2018-12-26 18:56 ` H.J. Lu
2018-12-17 16:02 ` [PATCH 1/3] " Wilco Dijkstra
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).