Consolidate multiple precision sin/cos functions

4 files changed, 140 insertions, 129 deletions

diff --git a/ChangeLog b/ChangeLog
index 5c8fef0..54f712d 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,30 @@
+2013-10-08  Siddhesh Poyarekar  <siddhesh@redhat.com>
+
+	* sysdeps/generic/math_private.h (__mpsin1): Remove
+	declaration.
+	(__mpcos1): Likewise.
+	(__mpsin): New argument __range_reduce.
+	(__mpcos): Likewise.
+	* sysdeps/ieee754/dbl-64/s_sin.c: Likewise.
+	(slow): Use __mpsin and __mpcos.
+	(slow1): Likewise.
+	(slow2): Likewise.
+	(sloww): Likewise.
+	(sloww1): Likewise.
+	(sloww2): Likewise.
+	(bsloww): Likewise.
+	(bsloww1): Likewise.
+	(bsloww2): Likewise.
+	(cslow2): Likewise.
+	(csloww): Likewise.
+	(csloww1): Likewise.
+	(csloww2): Likewise.
+	* sysdeps/ieee754/dbl-64/sincos32.c (__mpsin): Add argument
+	range_reduce.  Merge in __mpsin1.
+	(__mpcos): Likewise.
+	(__mpsin1): Remove.
+	(__mpcos1): Likewise.
+
 2013-10-07  Joseph Myers  <joseph@codesourcery.com>
 
 	* locale/loadlocale.c (_nl_intern_locale_data): Use
diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
index c0fc03d..9b881a3 100644
--- a/sysdeps/generic/math_private.h
+++ b/sysdeps/generic/math_private.h
@@ -356,10 +356,8 @@ extern void __dubcos (double __x, double __dx, double __v[]);
 extern double __halfulp (double __x, double __y);
 extern double __sin32 (double __x, double __res, double __res1);
 extern double __cos32 (double __x, double __res, double __res1);
-extern double __mpsin (double __x, double __dx);
-extern double __mpcos (double __x, double __dx);
-extern double __mpsin1 (double __x);
-extern double __mpcos1 (double __x);
+extern double __mpsin (double __x, double __dx, bool __range_reduce);
+extern double __mpcos (double __x, double __dx, bool __range_reduce);
 extern double __slowexp (double __x);
 extern double __slowpow (double __x, double __y, double __z);
 extern void __docos (double __x, double __dx, double __v[]);
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index 2be8fe3..b3dfef2 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -127,10 +127,8 @@ static const double
 
 void __dubsin (double x, double dx, double w[]);
 void __docos (double x, double dx, double w[]);
-double __mpsin (double x, double dx);
-double __mpcos (double x, double dx);
-double __mpsin1 (double x);
-double __mpcos1 (double x);
+double __mpsin (double x, double dx, bool reduce_range);
+double __mpcos (double x, double dx, bool reduce_range);
 static double slow (double x);
 static double slow1 (double x);
 static double slow2 (double x);
@@ -722,7 +720,7 @@ slow (double x)
       if (w[0] == w[0] + 1.000000001 * w[1])
 	return (x > 0) ? w[0] : -w[0];
       else
-	return (x > 0) ? __mpsin (x, 0) : -__mpsin (-x, 0);
+	return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
     }
 }
 
@@ -762,7 +760,7 @@ slow1 (double x)
       if (w[0] == w[0] + 1.000000005 * w[1])
 	return (x > 0) ? w[0] : -w[0];
       else
-	return (x > 0) ? __mpsin (x, 0) : -__mpsin (-x, 0);
+	return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
     }
 }
 
@@ -815,7 +813,7 @@ slow2 (double x)
       if (w[0] == w[0] + 1.000000005 * w[1])
 	return (x > 0) ? w[0] : -w[0];
       else
-	return (x > 0) ? __mpsin (x, 0) : -__mpsin (-x, 0);
+	return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false);
     }
 }
 
@@ -882,7 +880,7 @@ sloww (double x, double dx, double orig)
 	  if (w[0] == w[0] + cor)
 	    return (a > 0) ? w[0] : -w[0];
 	  else
-	    return __mpsin1 (orig);
+	    return __mpsin (orig, 0, true);
 	}
     }
 }
@@ -939,7 +937,7 @@ sloww1 (double x, double dx, double orig)
       if (w[0] == w[0] + cor)
 	return (x > 0) ? w[0] : -w[0];
       else
-	return __mpsin1 (orig);
+	return __mpsin (orig, 0, true);
     }
 }
 
@@ -996,7 +994,7 @@ sloww2 (double x, double dx, double orig, int n)
       if (w[0] == w[0] + cor)
 	return (n & 2) ? -w[0] : w[0];
       else
-	return __mpsin1 (orig);
+	return __mpsin (orig, 0, true);
     }
 }
 
@@ -1028,7 +1026,7 @@ bsloww (double x, double dx, double orig, int n)
       if (w[0] == w[0] + cor)
 	return (x > 0) ? w[0] : -w[0];
       else
-	return (n & 1) ? __mpcos1 (orig) : __mpsin1 (orig);
+	return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
     }
 }
 
@@ -1079,7 +1077,7 @@ bsloww1 (double x, double dx, double orig, int n)
       if (w[0] == w[0] + cor)
 	return (x > 0) ? w[0] : -w[0];
       else
-	return (n & 1) ? __mpcos1 (orig) : __mpsin1 (orig);
+	return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true);
     }
 }
 
@@ -1131,7 +1129,7 @@ bsloww2 (double x, double dx, double orig, int n)
       if (w[0] == w[0] + cor)
 	return (n & 2) ? -w[0] : w[0];
       else
-	return (n & 1) ? __mpsin1 (orig) : __mpcos1 (orig);
+	return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true);
     }
 }
 
@@ -1173,7 +1171,7 @@ cslow2 (double x)
       if (w[0] == w[0] + 1.000000005 * w[1])
 	return w[0];
       else
-	return __mpcos (x, 0);
+	return __mpcos (x, 0, false);
     }
 }
 
@@ -1246,7 +1244,7 @@ csloww (double x, double dx, double orig)
 	  if (w[0] == w[0] + cor)
 	    return (a > 0) ? w[0] : -w[0];
 	  else
-	    return __mpcos1 (orig);
+	    return __mpcos (orig, 0, true);
 	}
     }
 }
@@ -1301,7 +1299,7 @@ csloww1 (double x, double dx, double orig)
       if (w[0] == w[0] + cor)
 	return (x > 0) ? w[0] : -w[0];
       else
-	return __mpcos1 (orig);
+	return __mpcos (orig, 0, true);
     }
 }
 
@@ -1357,7 +1355,7 @@ csloww2 (double x, double dx, double orig, int n)
       if (w[0] == w[0] + cor)
 	return (n) ? -w[0] : w[0];
       else
-	return __mpcos1 (orig);
+	return __mpcos (orig, 0, true);
     }
 }
 
diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c
index ac27fcd..f253b8c 100644
--- a/sysdeps/ieee754/dbl-64/sincos32.c
+++ b/sysdeps/ieee754/dbl-64/sincos32.c
@@ -187,50 +187,119 @@ __cos32 (double x, double res, double res1)
     return (res < res1) ? res : res1;
 }
 
-/* Compute sin(x+dx) as Multi Precision number and return result as double.  */
+/* Compute sin() of double-length number (X + DX) as Multi Precision number and
+   return result as double.  If REDUCE_RANGE is true, X is assumed to be the
+   original input and DX is ignored.  */
 double
 SECTION
-__mpsin (double x, double dx)
+__mpsin (double x, double dx, bool reduce_range)
 {
-  int p;
   double y;
-  mp_no a, b, c;
-  p = 32;
-  __dbl_mp (x, &a, p);
-  __dbl_mp (dx, &b, p);
-  __add (&a, &b, &c, p);
-  if (x > 0.8)
+  mp_no a, b, c, s;
+  int n;
+  int p = 32;
+
+  if (reduce_range)
     {
-      __sub (&hp, &c, &a, p);
-      __c32 (&a, &b, &c, p);
+      n = __mpranred (x, &a, p);	/* n is 0, 1, 2 or 3.  */
+      __c32 (&a, &c, &s, p);
     }
   else
-    __c32 (&c, &a, &b, p);	/* b = sin(x+dx)  */
-  __mp_dbl (&b, &y, p);
+    {
+      n = -1;
+      __dbl_mp (x, &b, p);
+      __dbl_mp (dx, &c, p);
+      __add (&b, &c, &a, p);
+      if (x > 0.8)
+        {
+          __sub (&hp, &a, &b, p);
+          __c32 (&b, &s, &c, p);
+        }
+      else
+        __c32 (&a, &c, &s, p);	/* b = sin(x+dx)  */
+    }
+
+  /* Convert result based on which quarter of unit circle y is in.  */
+  switch (n)
+    {
+    case 1:
+      __mp_dbl (&c, &y, p);
+      break;
+
+    case 3:
+      __mp_dbl (&c, &y, p);
+      y = -y;
+      break;
+
+    case 2:
+      __mp_dbl (&s, &y, p);
+      y = -y;
+      break;
+
+    /* Quadrant not set, so the result must be sin (X + DX), which is also in
+       S.  */
+    case 0:
+    default:
+      __mp_dbl (&s, &y, p);
+    }
   return y;
 }
 
-/* Compute cos() of double-length number (x+dx) as Multi Precision number and
-   return result as double.  */
+/* Compute cos() of double-length number (X + DX) as Multi Precision number and
+   return result as double.  If REDUCE_RANGE is true, X is assumed to be the
+   original input and DX is ignored.  */
 double
 SECTION
-__mpcos (double x, double dx)
+__mpcos (double x, double dx, bool reduce_range)
 {
-  int p;
   double y;
-  mp_no a, b, c;
-  p = 32;
-  __dbl_mp (x, &a, p);
-  __dbl_mp (dx, &b, p);
-  __add (&a, &b, &c, p);
-  if (x > 0.8)
+  mp_no a, b, c, s;
+  int n;
+  int p = 32;
+
+  if (reduce_range)
     {
-      __sub (&hp, &c, &b, p);
-      __c32 (&b, &c, &a, p);
+      n = __mpranred (x, &a, p);	/* n is 0, 1, 2 or 3.  */
+      __c32 (&a, &c, &s, p);
     }
   else
-    __c32 (&c, &a, &b, p);	/* a = cos(x+dx)     */
-  __mp_dbl (&a, &y, p);
+    {
+      n = -1;
+      __dbl_mp (x, &b, p);
+      __dbl_mp (dx, &c, p);
+      __add (&b, &c, &a, p);
+      if (x > 0.8)
+        {
+          __sub (&hp, &a, &b, p);
+          __c32 (&b, &s, &c, p);
+        }
+      else
+        __c32 (&a, &c, &s, p);	/* a = cos(x+dx)     */
+    }
+
+  /* Convert result based on which quarter of unit circle y is in.  */
+  switch (n)
+    {
+    case 1:
+      __mp_dbl (&s, &y, p);
+      y = -y;
+      break;
+
+    case 3:
+      __mp_dbl (&s, &y, p);
+      break;
+
+    case 2:
+      __mp_dbl (&c, &y, p);
+      y = -y;
+      break;
+
+    /* Quadrant not set, so the result must be cos (X + DX), which is also
+       stored in C.  */
+    case 0:
+    default:
+      __mp_dbl (&c, &y, p);
+    }
   return y;
 }
 
@@ -294,84 +363,3 @@ __mpranred (double x, mp_no *y, int p)
       return (n & 3);
     }
 }
-
-/* Multi-Precision sin() function subroutine, for p = 32.  It is based on the
-   routines mpranred() and c32().  */
-double
-SECTION
-__mpsin1 (double x)
-{
-  int p;
-  int n;
-  mp_no u, s, c;
-  double y;
-  p = 32;
-  n = __mpranred (x, &u, p);	/* n is 0, 1, 2 or 3.  */
-  __c32 (&u, &c, &s, p);
-  /* Convert result based on which quarter of unit circle y is in.  */
-  switch (n)
-    {
-    case 0:
-      __mp_dbl (&s, &y, p);
-      return y;
-      break;
-
-    case 2:
-      __mp_dbl (&s, &y, p);
-      return -y;
-      break;
-
-    case 1:
-      __mp_dbl (&c, &y, p);
-      return y;
-      break;
-
-    case 3:
-      __mp_dbl (&c, &y, p);
-      return -y;
-      break;
-    }
-  /* Unreachable, to make the compiler happy.  */
-  return 0;
-}
-
-/* Multi-Precision cos() function subroutine, for p = 32.  It is based on the
-   routines mpranred() and c32().  */
-double
-SECTION
-__mpcos1 (double x)
-{
-  int p;
-  int n;
-  mp_no u, s, c;
-  double y;
-
-  p = 32;
-  n = __mpranred (x, &u, p);	/* n is 0, 1, 2 or 3.  */
-  __c32 (&u, &c, &s, p);
-  /* Convert result based on which quarter of unit circle y is in.  */
-  switch (n)
-    {
-    case 0:
-      __mp_dbl (&c, &y, p);
-      return y;
-      break;
-
-    case 2:
-      __mp_dbl (&c, &y, p);
-      return -y;
-      break;
-
-    case 1:
-      __mp_dbl (&s, &y, p);
-      return -y;
-      break;
-
-    case 3:
-      __mp_dbl (&s, &y, p);
-      return y;
-      break;
-    }
-  /* Unreachable, to make the compiler happy.  */
-  return 0;
-}