1 files changed, 44 insertions, 50 deletions

diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index c86fb9f2aa..b8c366a6f0 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -295,9 +295,13 @@ reduce_and_compute (double x, bool shift_quadrant)
   return retval;
 }
 
+/* Reduce range of x to within PI/2 with abs (x) < 105414350.  The high part
+   is written to *a, the low part to *da.  Range reduction is accurate to 136
+   bits so that when x is large and *a very close to zero, all 53 bits of *a
+   are correct.  */
 static inline int4
 __always_inline
-reduce_sincos_1 (double x, double *a, double *da)
+reduce_sincos (double x, double *a, double *da)
 {
   mynumber v;
 
@@ -306,62 +310,45 @@ reduce_sincos_1 (double x, double *a, double *da)
   v.x = t;
   double y = (x - xn * mp1) - xn * mp2;
   int4 n = v.i[LOW_HALF] & 3;
-  double db = xn * mp3;
-  double b = y - db;
-  db = (y - b) - db;
+
+  double b, db, t1, t2;
+  t1 = xn * pp3;
+  t2 = y - t1;
+  db = (y - t2) - t1;
+
+  t1 = xn * pp4;
+  b = t2 - t1;
+  db += (t2 - b) - t1;
 
   *a = b;
   *da = db;
-
   return n;
 }
 
-/* Compute sin (A + DA).  cos can be computed by passing SHIFT_QUADRANT as
-   true, which results in shifting the quadrant N clockwise.  */
+/* Compute sin or cos (A + DA) for the given quadrant N.  */
 static double
 __always_inline
-do_sincos_1 (double a, double da, double x, int4 n, bool shift_quadrant)
+do_sincos (double a, double da, int4 n)
 {
-  double xx, retval, res, cor;
-  double eps = fabs (x) * 1.2e-30;
+  double retval, cor;
 
-  int k1 = (n + shift_quadrant) & 3;
-  switch (k1)
-    {			/* quarter of unit circle */
-    case 2:
-      a = -a;
-      da = -da;
-      /* Fall through.  */
-    case 0:
-      xx = a * a;
+  if (n & 1)
+    /* Max ULP is 0.513.  */
+    retval = do_cos (a, da, &cor);
+  else
+    {
+      double xx = a * a;
+      /* Max ULP is 0.501 if xx < 0.01588, otherwise ULP is 0.518.  */
       if (xx < 0.01588)
-	{
-	  /* Taylor series.  */
-	  res = TAYLOR_SIN (xx, a, da, cor);
-	  cor = 1.02 * cor + __copysign (eps, cor);
-	  retval = (res == res + cor) ? res : sloww (a, da, x, shift_quadrant);
-	}
+	retval = TAYLOR_SIN (xx, a, da, cor);
       else
-	{
-	  res = do_sin (a, da, &cor);
-	  cor = 1.035 * cor + __copysign (eps, cor);
-	  retval = ((res == res + cor) ? __copysign (res, a)
-		    : sloww1 (a, da, x, shift_quadrant));
-	}
-      break;
-
-    case 1:
-    case 3:
-      res = do_cos (a, da, &cor);
-      cor = 1.025 * cor + __copysign (eps, cor);
-      retval = ((res == res + cor) ? ((n & 2) ? -res : res)
-		: sloww2 (a, da, x, n));
-      break;
+	retval = __copysign (do_sin (a, da, &cor), a);
     }
 
-  return retval;
+  return (n & 2) ? -retval : retval;
 }
 
+
 /*******************************************************************/
 /* An ultimate sin routine. Given an IEEE double machine number x   */
 /* it computes the correctly rounded (to nearest) value of sin(x)  */
@@ -374,13 +361,18 @@ SECTION
 #endif
 __sin (double x)
 {
-  double xx, t, cor;
+#ifndef IN_SINCOS
+  double xx, t, a, da, cor;
   mynumber u;
-  int4 k, m;
+  int4 k, m, n;
   double retval = 0;
 
-#ifndef IN_SINCOS
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
+#else
+  double xx, t, cor;
+  mynumber u;
+  int4 k, m;
+  double retval = 0;
 #endif
 
   u.x = x;
@@ -419,9 +411,8 @@ __sin (double x)
 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
   else if (k < 0x419921FB)
     {
-      double a, da;
-      int4 n = reduce_sincos_1 (x, &a, &da);
-      retval = do_sincos_1 (a, da, x, n, false);
+      n = reduce_sincos (x, &a, &da);
+      retval = do_sincos (a, da, n);
     }				/*   else  if (k <  0x419921FB )    */
 
 /* --------------------105414350 <|x| <2^1024------------------------------*/
@@ -456,7 +447,11 @@ __cos (double x)
 {
   double y, xx, cor, a, da;
   mynumber u;
+#ifndef IN_SINCOS
+  int4 k, m, n;
+#else
   int4 k, m;
+#endif
 
   double retval = 0;
 
@@ -496,9 +491,8 @@ __cos (double x)
 #ifndef IN_SINCOS
   else if (k < 0x419921FB)
     {				/* 2.426265<|x|< 105414350 */
-      double a, da;
-      int4 n = reduce_sincos_1 (x, &a, &da);
-      retval = do_sincos_1 (a, da, x, n, true);
+      n = reduce_sincos (x, &a, &da);
+      retval = do_sincos (a, da, n + 1);
     }				/*   else  if (k <  0x419921FB )    */
 
   /* 105414350 <|x| <2^1024 */