Remove slow paths from pow

Remove the slow paths from pow.  Like several other double precision math
functions, pow is exactly rounded.  This is not required from math functions
and causes major overheads as it requires multiple fallbacks using higher
precision arithmetic if a result is close to 0.5ULP.  Ridiculous slowdowns
of up to 100000x have been reported when the highest precision path triggers.

All GLIBC math tests pass on AArch64 and x64 (with ULP of pow set to 1).
The worst case error is ~0.506ULP.  A simple test over a few hundred million
values shows pow is 10% faster on average.  This fixes BZ #13932.

	[BZ #13932]
	* sysdeps/ieee754/dbl-64/uexp.h (err_1): Remove.
	* benchtests/pow-inputs: Update comment for slow path cases.
	* manual/probes.texi (slowpow_p10): Delete removed probe.
	(slowpow_p10): Likewise.
	* math/Makefile: Remove halfulp.c and slowpow.c.
	* sysdeps/aarch64/libm-test-ulps: Set ULP of pow to 1.
	* sysdeps/generic/math_private.h (__exp1): Remove error argument.
	(__halfulp): Remove.
	(__slowpow): Remove.
	* sysdeps/i386/fpu/halfulp.c: Delete file.
	* sysdeps/i386/fpu/slowpow.c: Likewise.
	* sysdeps/ia64/fpu/halfulp.c: Likewise.
	* sysdeps/ia64/fpu/slowpow.c: Likewise.
	* sysdeps/ieee754/dbl-64/e_exp.c (__exp1): Remove error argument,
	improve comments and add error analysis.
	* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Add error analysis.
	(power1): Remove function:
	(log1): Remove error argument, add error analysis.
	(my_log2): Remove function.
	* sysdeps/ieee754/dbl-64/halfulp.c: Delete file.
	* sysdeps/ieee754/dbl-64/slowpow.c: Likewise.
	* sysdeps/m68k/m680x0/fpu/halfulp.c: Likewise.
	* sysdeps/m68k/m680x0/fpu/slowpow.c: Likewise.
	* sysdeps/powerpc/power4/fpu/Makefile: Remove CPPFLAGS-slowpow.c.
	* sysdeps/x86_64/fpu/libm-test-ulps: Set ULP of pow to 1.
	* sysdeps/x86_64/fpu/multiarch/Makefile: Remove slowpow-fma.c,
	slowpow-fma4.c, halfulp-fma.c, halfulp-fma4.c.
	* sysdeps/x86_64/fpu/multiarch/e_pow-fma.c (__slowpow): Remove define.
	* sysdeps/x86_64/fpu/multiarch/e_pow-fma4.c (__slowpow): Likewise.
	* sysdeps/x86_64/fpu/multiarch/halfulp-fma.c: Delete file.
	* sysdeps/x86_64/fpu/multiarch/halfulp-fma4.c: Likewise.
	* sysdeps/x86_64/fpu/multiarch/slowpow-fma.c: Likewise.
	* sysdeps/x86_64/fpu/multiarch/slowpow-fma4.c: Likewise.

1 files changed, 24 insertions, 148 deletions

diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index f6e5fcdbce..542d03a7e3 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -20,13 +20,9 @@
 /*  MODULE_NAME: upow.c                                                    */
 /*                                                                         */
 /*  FUNCTIONS: upow                                                        */
-/*             power1                                                      */
-/*             my_log2                                                     */
 /*             log1                                                        */
 /*             checkint                                                    */
 /* FILES NEEDED: dla.h endian.h mpa.h mydefs.h                             */
-/*               halfulp.c mpexp.c mplog.c slowexp.c slowpow.c mpa.c       */
-/*                          uexp.c  upow.c				   */
 /*               root.tbl uexp.tbl upow.tbl                                */
 /* An ultimate power routine. Given two IEEE double machine numbers y,x    */
 /* it computes the correctly rounded (to nearest) value of x^y.            */
@@ -50,11 +46,8 @@
 
 static const double huge = 1.0e300, tiny = 1.0e-300;
 
-double __exp1 (double x, double xx, double error);
-static double log1 (double x, double *delta, double *error);
-static double my_log2 (double x, double *delta, double *error);
-double __slowpow (double x, double y, double z);
-static double power1 (double x, double y);
+double __exp1 (double x, double xx);
+static double log1 (double x, double *delta);
 static int checkint (double x);
 
 /* An ultimate power routine. Given two IEEE double machine numbers y, x it
@@ -63,7 +56,7 @@ double
 SECTION
 __ieee754_pow (double x, double y)
 {
-  double z, a, aa, error, t, a1, a2, y1, y2;
+  double z, a, aa, t, a1, a2, y1, y2;
   mynumber u, v;
   int k;
   int4 qx, qy;
@@ -100,7 +93,7 @@ __ieee754_pow (double x, double y)
 	   not matter if |y| <= 2**-64.  */
 	if (fabs (y) < 0x1p-64)
 	  y = y < 0 ? -0x1p-64 : 0x1p-64;
-	z = log1 (x, &aa, &error);	/* x^y  =e^(y log (X)) */
+	z = log1 (x, &aa);	/* x^y  =e^(y log (X)) */
 	t = y * CN;
 	y1 = t - (t - y);
 	y2 = y - y1;
@@ -111,9 +104,16 @@ __ieee754_pow (double x, double y)
 	aa = y2 * a1 + y * a2;
 	a1 = a + aa;
 	a2 = (a - a1) + aa;
-	error = error * fabs (y);
-	t = __exp1 (a1, a2, 1.9e16 * error);	/* return -10 or 0 if wasn't computed exactly */
-	retval = (t > 0) ? t : power1 (x, y);
+
+	/* Maximum relative error RElog of log1 is 1.0e-21 (69.7 bits).
+	   Maximum relative error REexp of __exp1 is 8.8e-22 (69.9 bits).
+	   We actually compute exp ((1 + RElog) * log (x) * y) * (1 + REexp).
+	   Since RElog/REexp are tiny and log (x) * y is at most log (DBL_MAX),
+	   this is equivalent to pow (x, y) * (1 + 710 * RElog + REexp).
+	   So the relative error is 710 * 1.0e-21 + 8.8e-22 = 7.1e-19
+	   (60.2 bits).  The worst-case ULP error is 0.5064.  */
+
+	retval = __exp1 (a1, a2);
       }
 
       if (isinf (retval))
@@ -218,33 +218,11 @@ __ieee754_pow (double x, double y)
 strong_alias (__ieee754_pow, __pow_finite)
 #endif
 
-/* Compute x^y using more accurate but more slow log routine.  */
-static double
-SECTION
-power1 (double x, double y)
-{
-  double z, a, aa, error, t, a1, a2, y1, y2;
-  z = my_log2 (x, &aa, &error);
-  t = y * CN;
-  y1 = t - (t - y);
-  y2 = y - y1;
-  t = z * CN;
-  a1 = t - (t - z);
-  a2 = z - a1;
-  a = y * z;
-  aa = ((y1 * a1 - a) + y1 * a2 + y2 * a1) + y2 * a2 + aa * y;
-  a1 = a + aa;
-  a2 = (a - a1) + aa;
-  error = error * fabs (y);
-  t = __exp1 (a1, a2, 1.9e16 * error);
-  return (t >= 0) ? t : __slowpow (x, y, z);
-}
-
 /* Compute log(x) (x is left argument). The result is the returned double + the
-   parameter DELTA.  The result is bounded by ERROR.  */
+   parameter DELTA.  */
 static double
 SECTION
-log1 (double x, double *delta, double *error)
+log1 (double x, double *delta)
 {
   unsigned int i, j;
   int m;
@@ -260,9 +238,7 @@ log1 (double x, double *delta, double *error)
 
   u.x = x;
   m = u.i[HIGH_HALF];
-  *error = 0;
-  *delta = 0;
-  if (m < 0x00100000)		/*  1<x<2^-1007 */
+  if (m < 0x00100000)		/* Handle denormal x.  */
     {
       x = x * t52.x;
       add = -52.0;
@@ -284,7 +260,7 @@ log1 (double x, double *delta, double *error)
   v.x = u.x + bigu.x;
   uu = v.x - bigu.x;
   i = (v.i[LOW_HALF] & 0x000003ff) << 2;
-  if (two52.i[LOW_HALF] == 1023)	/* nx = 0              */
+  if (two52.i[LOW_HALF] == 1023)	/* Exponent of x is 0.  */
     {
       if (i > 1192 && i < 1208)	/* |x-1| < 1.5*2**-10  */
 	{
@@ -296,8 +272,8 @@ log1 (double x, double *delta, double *error)
 							   * (r7 + t * r8)))))
 		- 0.5 * t2 * (t + t1));
 	  res = e1 + e2;
-	  *error = 1.0e-21 * fabs (t);
 	  *delta = (e1 - res) + e2;
+	  /* Max relative error is 1.464844e-24, so accurate to 79.1 bits.  */
 	  return res;
 	}			/* |x-1| < 1.5*2**-10  */
       else
@@ -316,12 +292,12 @@ log1 (double x, double *delta, double *error)
 	  t2 = ((((t - t1) + e) + (ui.x[i + 3] + vj.x[j + 2])) + e2 + e * e
 		* (p2 + e * (p3 + e * p4)));
 	  res = t1 + t2;
-	  *error = 1.0e-24;
 	  *delta = (t1 - res) + t2;
+	  /* Max relative error is 1.0e-24, so accurate to 79.7 bits.  */
 	  return res;
 	}
-    }				/* nx = 0 */
-  else				/* nx != 0   */
+    }
+  else				/* Exponent of x != 0.  */
     {
       eps = u.x - uu;
       nx = (two52.x - two52e.x) + add;
@@ -334,113 +310,13 @@ log1 (double x, double *delta, double *error)
       t2 = ((((t - t1) + e) + nx * ln2b.x + ui.x[i + 3] + e2) + e * e
 	    * (q2 + e * (q3 + e * (q4 + e * (q5 + e * q6)))));
       res = t1 + t2;
-      *error = 1.0e-21;
-      *delta = (t1 - res) + t2;
-      return res;
-    }				/* nx != 0   */
-}
-
-/* Slower but more accurate routine of log.  The returned result is double +
-   DELTA.  The result is bounded by ERROR.  */
-static double
-SECTION
-my_log2 (double x, double *delta, double *error)
-{
-  unsigned int i, j;
-  int m;
-  double uu, vv, eps, nx, e, e1, e2, t, t1, t2, res, add = 0;
-  double ou1, ou2, lu1, lu2, ov, lv1, lv2, a, a1, a2;
-  double y, yy, z, zz, j1, j2, j7, j8;
-#ifndef DLA_FMS
-  double j3, j4, j5, j6;
-#endif
-  mynumber u, v;
-#ifdef BIG_ENDI
-  mynumber /**/ two52 = {{0x43300000, 0x00000000}};	/* 2**52  */
-#else
-# ifdef LITTLE_ENDI
-  mynumber /**/ two52 = {{0x00000000, 0x43300000}};	/* 2**52  */
-# endif
-#endif
-
-  u.x = x;
-  m = u.i[HIGH_HALF];
-  *error = 0;
-  *delta = 0;
-  add = 0;
-  if (m < 0x00100000)
-    {				/* x < 2^-1022 */
-      x = x * t52.x;
-      add = -52.0;
-      u.x = x;
-      m = u.i[HIGH_HALF];
-    }
-
-  if ((m & 0x000fffff) < 0x0006a09e)
-    {
-      u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3ff00000;
-      two52.i[LOW_HALF] = (m >> 20);
-    }
-  else
-    {
-      u.i[HIGH_HALF] = (m & 0x000fffff) | 0x3fe00000;
-      two52.i[LOW_HALF] = (m >> 20) + 1;
-    }
-
-  v.x = u.x + bigu.x;
-  uu = v.x - bigu.x;
-  i = (v.i[LOW_HALF] & 0x000003ff) << 2;
-  /*------------------------------------- |x-1| < 2**-11-------------------------------  */
-  if ((two52.i[LOW_HALF] == 1023) && (i == 1200))
-    {
-      t = x - 1.0;
-      EMULV (t, s3, y, yy, j1, j2, j3, j4, j5);
-      ADD2 (-0.5, 0, y, yy, z, zz, j1, j2);
-      MUL2 (t, 0, z, zz, y, yy, j1, j2, j3, j4, j5, j6, j7, j8);
-      MUL2 (t, 0, y, yy, z, zz, j1, j2, j3, j4, j5, j6, j7, j8);
-
-      e1 = t + z;
-      e2 = ((((t - e1) + z) + zz) + t * t * t
-	    * (ss3 + t * (s4 + t * (s5 + t * (s6 + t * (s7 + t * s8))))));
-      res = e1 + e2;
-      *error = 1.0e-25 * fabs (t);
-      *delta = (e1 - res) + e2;
-      return res;
-    }
-  /*----------------------------- |x-1| > 2**-11  --------------------------  */
-  else
-    {				/*Computing log(x) according to log table                        */
-      nx = (two52.x - two52e.x) + add;
-      ou1 = ui.x[i];
-      ou2 = ui.x[i + 1];
-      lu1 = ui.x[i + 2];
-      lu2 = ui.x[i + 3];
-      v.x = u.x * (ou1 + ou2) + bigv.x;
-      vv = v.x - bigv.x;
-      j = v.i[LOW_HALF] & 0x0007ffff;
-      j = j + j + j;
-      eps = u.x - uu * vv;
-      ov = vj.x[j];
-      lv1 = vj.x[j + 1];
-      lv2 = vj.x[j + 2];
-      a = (ou1 + ou2) * (1.0 + ov);
-      a1 = (a + 1.0e10) - 1.0e10;
-      a2 = a * (1.0 - a1 * uu * vv);
-      e1 = eps * a1;
-      e2 = eps * a2;
-      e = e1 + e2;
-      e2 = (e1 - e) + e2;
-      t = nx * ln2a.x + lu1 + lv1;
-      t1 = t + e;
-      t2 = ((((t - t1) + e) + (lu2 + lv2 + nx * ln2b.x + e2)) + e * e
-	    * (p2 + e * (p3 + e * p4)));
-      res = t1 + t2;
-      *error = 1.0e-27;
       *delta = (t1 - res) + t2;
+      /* Max relative error is 1.0e-21, so accurate to 69.7 bits.  */
       return res;
     }
 }
 
+
 /* This function receives a double x and checks if it is an integer.  If not,
    it returns 0, else it returns 1 if even or -1 if odd.  */
 static int