x86-64: Vectorize sincosf_poly and update s_sincosf-fma.c

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>.

2 files changed, 88 insertions, 70 deletions

diff --git a/sysdeps/ieee754/flt-32/s_sincosf.h b/sysdeps/ieee754/flt-32/s_sincosf.h
index 1dcb04f..54a9da5 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 0000000..f49c75a
--- /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;
+    }
+}