Consolidate range reduction in sincos for x > 281474976710656

Range reduction needs to be done only once for sin and cos, so copy
over all of the relevant functions (__sin, __cos, reduce_and_compute)
and consolidate common code.

2 files changed, 57 insertions, 2 deletions

diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index a635a86933..b619905dcb 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -497,6 +497,7 @@ __sin (double x)
 	}
     }				/*   else  if (k <  0x42F00000 )   */
 
+#ifndef IN_SINCOS
 /* -----------------281474976710656 <|x| <2^1024----------------------------*/
   else if (k < 0x7ff00000)
     retval = reduce_and_compute (x, 0);
@@ -508,6 +509,7 @@ __sin (double x)
 	__set_errno (EDOM);
       retval = x / x;
     }
+#endif
 
   return retval;
 }
@@ -729,6 +731,7 @@ __cos (double x)
 	}
     }				/*   else  if (k <  0x42F00000 )    */
 
+#ifndef IN_SINCOS
   /* 281474976710656 <|x| <2^1024 */
   else if (k < 0x7ff00000)
     retval = reduce_and_compute (x, 1);
@@ -739,6 +742,7 @@ __cos (double x)
 	__set_errno (EDOM);
       retval = x / x;		/* |x| > 2^1024 */
     }
+#endif
 
   return retval;
 }
diff --git a/sysdeps/ieee754/dbl-64/s_sincos.c b/sysdeps/ieee754/dbl-64/s_sincos.c
index 2a3fc06386..f50ffa625f 100644
--- a/sysdeps/ieee754/dbl-64/s_sincos.c
+++ b/sysdeps/ieee754/dbl-64/s_sincos.c
@@ -27,13 +27,64 @@
 #define IN_SINCOS 1
 #include "s_sin.c"
 
+/* Consolidated version of reduce_and_compute in s_sin.c that does range
+   reduction only once and computes sin and cos together.  */
+static inline void
+__always_inline
+reduce_and_compute_sincos (double x, double *sinx, double *cosx)
+{
+  double a, da;
+  unsigned int n = __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;
+    }
+
+  if (a * a < 0.01588)
+    *sinx = bsloww (a, da, x, n);
+  else
+    *sinx = bsloww1 (a, da, x, n);
+  *cosx = bsloww2 (a, da, x, n);
+}
+
 void
 __sincos (double x, double *sinx, double *cosx)
 {
+  mynumber u;
+  int k;
+
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
 
-  *sinx = __sin (x);
-  *cosx = __cos (x);
+  u.x = x;
+  k = 0x7fffffff & u.i[HIGH_HALF];
+
+  if (k < 0x42F00000)
+    {
+      *sinx = __sin_local (x);
+      *cosx = __cos_local (x);
+      return;
+    }
+  if (k < 0x7ff00000)
+    {
+      reduce_and_compute_sincos (x, sinx, cosx);
+      return;
+    }
+
+  if (isinf (x))
+    __set_errno (EDOM);
+
+  *sinx = *cosx = x / x;
 }
 weak_alias (__sincos, sincos)
 #ifdef NO_LONG_DOUBLE