Patchwork [3/6] Remove slow paths from sin/cos

login
register
mail settings
Submitter Wilco Dijkstra
Date March 9, 2018, 3:48 p.m.
Message ID <DB6PR0801MB205347ED2AF8612B96B257AE83DE0@DB6PR0801MB2053.eurprd08.prod.outlook.com>
Download mbox | patch
Permalink /patch/26252/
State New
Headers show

Comments

Wilco Dijkstra - March 9, 2018, 3:48 p.m.
This patch improves the accuracy of the range reduction.  When the input is
large (2^27) and very close to a multiple of PI/2, using 110 bits of PI is not
enough.  Improve range reduction accuracy to 136 bits.  As a result the special
checks for results close to zero can be removed.  The ULP of the polynomials is
at worst 0.55ULP, so there is no reason for the slow functions, and they can be
removed.

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

	* sysdeps/ieee754/dbl-64/s_sin.c (reduce_sincos_1): Rename to sincos_1,
	improve accuracy to 136 bits.
	(do_sincos_1): Rename to do_sincos, remove fallbacks to slow functions.
	(__sin): Use improved reduction and simplified do_sincos calculation.
	(__cos): Likewise.
	* sysdeps/ieee754/dbl-64/s_sincos.c (__sincos): Likewise.

--
Joseph Myers - March 15, 2018, 4:49 p.m.
On Fri, 9 Mar 2018, Wilco Dijkstra wrote:

> +#else
> +    /* Disable warning... */
> +    n = 0, n = n;
> +    a = 0, a = a;
> +    da = 0, da = da;
>  #endif

I don't think we should use this approach for warning avoidance - rather, 
the variable declarations should be conditioned appropriately, or if 
what's conditional is whether the variables are used or only set, there 
should be conditional __attribute__ ((unused)) on the declarations.

(I realise patch 6 would remove this code, but think each patch should 
leave the tree in a sensible state - especially given patch 6 is evidently 
one of the more controversial pieces of this series.)

> +#else
> +    /* Disable warning... */
> +    n = 0, n = n;
>  #endif

Likewise.

Patch

diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index 1f98e29278183d1fccd7c2b3fd467d6b16c245ed..b48b8627a7a801dfafecc920062aaaac51969a8a 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,9 +361,9 @@  SECTION
 #endif
 __sin (double x)
 {
-  double xx, t, cor;
+  double xx, t, a, da, cor;
   mynumber u;
-  int4 k, m;
+  int4 k, m, n;
   double retval = 0;
 
 #ifndef IN_SINCOS
@@ -419,9 +406,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------------------------------*/
@@ -435,6 +421,11 @@  __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;
@@ -456,7 +447,7 @@  __cos (double x)
 {
   double y, xx, cor, a, da;
   mynumber u;
-  int4 k, m;
+  int4 k, m, n;
 
   double retval = 0;
 
@@ -496,9 +487,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 */
@@ -511,6 +501,9 @@  __cos (double x)
 	__set_errno (EDOM);
       retval = x / x;		/* |x| > 2^1024 */
     }
+#else
+    /* Disable warning... */
+    n = 0, n = n;
 #endif
 
   return retval;
diff --git a/sysdeps/ieee754/dbl-64/s_sincos.c b/sysdeps/ieee754/dbl-64/s_sincos.c
index a9af8ce526bfe78c06cfafa65de0815ec69585c5..4f032d2e42593ccde22169b374728386dd8fca8e 100644
--- a/sysdeps/ieee754/dbl-64/s_sincos.c
+++ b/sysdeps/ieee754/dbl-64/s_sincos.c
@@ -79,10 +79,10 @@  __sincos (double x, double *sinx, double *cosx)
   if (k < 0x419921FB)
     {
       double a, da;
-      int4 n = reduce_sincos_1 (x, &a, &da);
+      int4 n = reduce_sincos (x, &a, &da);
 
-      *sinx = do_sincos_1 (a, da, x, n, false);
-      *cosx = do_sincos_1 (a, da, x, n, true);
+      *sinx = do_sincos (a, da, n);
+      *cosx = do_sincos (a, da, n + 1);
 
       return;
     }