[PATCH 6/6] Remove slow paths from sin/cos

[Wilco Dijkstra <Wilco.Dijkstra@arm.com>]

Restructure the sincos implementation - rather than rely on odd partial inlining
of preprocessed portions from sin and cos, explicitly write out the cases.

ChangeLog:
2018-03-09  Wilco Dijkstra  <wdijkstr@arm.com>

	* sysdeps/ieee754/dbl-64/s_sin.c (__sin): Cleanup ifdefs.
	(__cos): Likewise.
	* sysdeps/ieee754/dbl-64/s_sin.c (__sincos): Reimplement using the same
	logic as sin and cos.

--
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index 91b0abc9c3ac21dae0e673576940ef97bfd20c23..8f804a42e6d94652a62f81b2e0b053135cf9f03a 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -208,12 +208,9 @@ do_sincos (double a, double da, int4 n)
 /* An ultimate sin routine. Given an IEEE double machine number x   */
 /* it computes the correctly rounded (to nearest) value of sin(x)  */
 /*******************************************************************/
-#ifdef IN_SINCOS
-static double
-#else
+#ifndef IN_SINCOS
 double
 SECTION
-#endif
 __sin (double x)
 {
   double xx, t, a, da;
@@ -221,9 +218,7 @@ __sin (double x)
   int4 k, m, n;
   double retval = 0;
 
-#ifndef IN_SINCOS
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
-#endif
 
   u.x = x;
   m = u.i[HIGH_HALF];
@@ -257,7 +252,6 @@ __sin (double x)
       retval = __copysign (do_cos (t, hp1), x);
     }				/*   else  if (k < 0x400368fd)    */
 
-#ifndef IN_SINCOS
 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
   else if (k < 0x419921FB)
     {
@@ -278,12 +272,6 @@ __sin (double x)
 	__set_errno (EDOM);
       retval = x / x;
     }
-#else
-    /* Disable warning... */
-    n = 0, n = n;
-    a = 0, a = a;
-    da = 0, da = da;
-#endif
 
   return retval;
 }
@@ -294,12 +282,8 @@ __sin (double x)
 /* it computes the correctly rounded (to nearest) value of cos(x)  */
 /*******************************************************************/
 
-#ifdef IN_SINCOS
-static double
-#else
 double
 SECTION
-#endif
 __cos (double x)
 {
   double y, xx, a, da;
@@ -308,9 +292,7 @@ __cos (double x)
 
   double retval = 0;
 
-#ifndef IN_SINCOS
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
-#endif
 
   u.x = x;
   m = u.i[HIGH_HALF];
@@ -340,8 +322,6 @@ __cos (double x)
 	retval = __copysign (do_sin (a, da), a);
     }				/*   else  if (k < 0x400368fd)    */
 
-
-#ifndef IN_SINCOS
   else if (k < 0x419921FB)
     {				/* 2.426265<|x|< 105414350 */
       n = reduce_sincos (x, &a, &da);
@@ -361,10 +341,6 @@ __cos (double x)
 	__set_errno (EDOM);
       retval = x / x;		/* |x| > 2^1024 */
     }
-#else
-    /* Disable warning... */
-    n = 0, n = n;
-#endif
 
   return retval;
 }
@@ -375,3 +351,5 @@ libm_alias_double (__cos, cos)
 #ifndef __sin
 libm_alias_double (__sin, sin)
 #endif
+
+#endif
diff --git a/sysdeps/ieee754/dbl-64/s_sincos.c b/sysdeps/ieee754/dbl-64/s_sincos.c
index 4335ecbba3c9894e61c087ac970b392fa73abfab..c04972707b284e37b15e82933a00250cda959985 100644
--- a/sysdeps/ieee754/dbl-64/s_sincos.c
+++ b/sysdeps/ieee754/dbl-64/s_sincos.c
@@ -23,9 +23,7 @@
 #include <math_private.h>
 #include <libm-alias-double.h>
 
-#define __sin __sin_local
-#define __cos __cos_local
-#define IN_SINCOS 1
+#define IN_SINCOS
 #include "s_sin.c"
 
 void
@@ -37,31 +35,79 @@ __sincos (double x, double *sinx, double *cosx)
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
 
   u.x = x;
-  k = 0x7fffffff & u.i[HIGH_HALF];
+  k = u.i[HIGH_HALF] & 0x7fffffff;
 
   if (k < 0x400368fd)
     {
-      *sinx = __sin_local (x);
-      *cosx = __cos_local (x);
-      return;
-    }
-  if (k < 0x419921FB)
-    {
-      double a, da;
-      int4 n = reduce_sincos (x, &a, &da);
-
-      *sinx = do_sincos (a, da, n);
-      *cosx = do_sincos (a, da, n + 1);
+      double t, xx, a, da, y;
+      /* |x| < 2^-27 => cos (x) = 1, sin (x) = x.  */
+      if (k < 0x3e400000)
+	{
+	  if (k < 0x3e500000)
+	    math_check_force_underflow (x);
+	  *sinx = x;
+	  *cosx = 1.0;
+	  return;
+	}
+      /* |x| < 0.855469.  */
+      else if (k < 0x3feb6000)
+	{
+	  /* |x| < 0.25.  */
+	  if (k < 0x3fd00000)
+	    {
+	      xx = x * x;
+	      t = POLYNOMIAL (xx) * (xx * x);
+	      *sinx = x + t;
+	    }
+	  else
+	    *sinx = __copysign (do_sin (x, 0), x);
+	  *cosx = do_cos (x, 0);
+	  return;
+	}
 
+      /* |x| < 2.426265.  */
+      y = hp0 - fabs (x);
+      a = y + hp1;
+      da = (y - a) + hp1;
+      *sinx = __copysign (do_cos (a, da), x);
+      xx = a * a;
+      if (xx < 0.01588)
+	*cosx = TAYLOR_SIN (xx, a, da);
+      else
+	*cosx = __copysign (do_sin (a, da), a);
       return;
     }
+  /* |x| < 2^1024.  */
   if (k < 0x7ff00000)
     {
-      double a, da;
-      int4 n = __branred (x, &a, &da);
+      double a, da, xx;
+      unsigned int n;
 
-      *sinx = do_sincos (a, da, n);
-      *cosx = do_sincos (a, da, n + 1);
+      /* If |x| < 105414350 use simple range reduction.  */
+      n = k < 0x419921FB ? reduce_sincos (x, &a, &da) : __branred (x, &a, &da);
+      n = n & 3;
+
+      if (n == 1 || n == 2)
+	{
+	  a = -a;
+	  da = -da;
+	}
+
+      if (n & 1)
+	{
+	  double *temp = cosx;
+	  cosx = sinx;
+	  sinx = temp;
+	}
+
+      xx = a * a;
+      if (xx < 0.01588)
+	*sinx = TAYLOR_SIN (xx, a, da);
+      else
+	*sinx = __copysign (do_sin (a, da), a);
+      xx = do_cos (a, da);
+      *cosx = (n & 2) ? -xx : xx;
+      return;
     }
 
   if (isinf (x))