Improve __ieee754_exp() performance by greater than 5x on sparc/x86.

These changes will be active for all platforms that don't provide
their own exp() routines. They will also be active for ieee754
versions of ccos, ccosh, cosh, csin, csinh, sinh, exp10, gamma, and
erf.

Typical performance gains is typically around 5x when measured on
Sparc s7 for common values between exp(1) and exp(40).

Using the glibc perf tests on sparc,
      sparc (nsec)    x86 (nsec)
      old     new     old     new
max   17629   395    5173     144
min     399    54      15      13
mean   5317   200    1349      23

The extreme max times for the old (ieee754) exp are due to the
multiprecision computation in the old algorithm when the true value is
very near 0.5 ulp away from an value representable in double
precision. The new algorithm does not take special measures for those
cases. The current glibc exp perf tests overrepresent those values.
Informal testing suggests approximately one in 200 cases might
invoke the high cost computation. The performance advantage of the new
algorithm for other values is still large but not as large as indicated
by the chart above.

Glibc correctness tests for exp() and expf() were run. Within the
test suite 3 input values were found to cause 1 bit differences (ulp)
when "FE_TONEAREST" rounding mode is set. No differences in exp() were
seen for the tested values for the other rounding modes.
Typical example:
exp(-0x1.760cd2p+0)  (-1.46113312244415283203125)
 new code:    2.31973271630014299393707e-01   0x1.db14cd799387ap-3
 old code:    2.31973271630014271638132e-01   0x1.db14cd7993879p-3
    exp    =  2.31973271630014285508337 (high precision)
Old delta: off by 0.49 ulp
New delta: off by 0.51 ulp

In addition, because ieee754_exp() is used by other routines, cexp()
showed test results with very small imaginary input values where the
imaginary portion of the result was off by 3 ulp when in upward
rounding mode, but not in the other rounding modes.  For x86, tgamma
showed a few values where the ulp increased to 6 (max ulp for tgamma
is 5). Sparc tgamma did not show these failures.  I presume the tgamma
differences are due to compiler optimization differences within the
gamma function.The gamma function is known to be difficult to compute
accurately.

	* sysdeps/ieee754/dbl-64/e_exp.c: Include <math-svid-compat.h> and
	<errno.h>.  Include "eexp.tbl".
	(half): New constant.
	(one): Likewise.
	(__ieee754_exp): Rewrite.
	(__slowexp): Remove prototype.
	* sysdeps/ieee754/dbl-64/eexp.tbl: New file.
	* sysdeps/ieee754/dbl-64/slowexp.c: Remove file.
	* sysdeps/i386/fpu/slowexp.c: Likewise.
	* sysdeps/ia64/fpu/slowexp.c: Likewise.
	* sysdeps/m68k/m680x0/fpu/slowexp.c: Likewise.
	* sysdeps/x86_64/fpu/multiarch/slowexp-avx.c: Likewise.
	* sysdeps/x86_64/fpu/multiarch/slowexp-fma.c: Likewise.
	* sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c: Likewise.
	* sysdeps/generic/math_private.h (__slowexp): Remove prototype.
	* sysdeps/ieee754/dbl-64/e_pow.c: Remove mention of slowexp.c in
	comment.
	* sysdeps/powerpc/power4/fpu/Makefile [$(subdir) = math]
	(CPPFLAGS-slowexp.c): Remove variable.
	* sysdeps/x86_64/fpu/multiarch/Makefile (libm-sysdep_routines):
	Remove slowexp-fma, slowexp-fma4 and slowexp-avx.
	(CFLAGS-slowexp-fma.c): Remove variable.
	(CFLAGS-slowexp-fma4.c): Likewise.
	(CFLAGS-slowexp-avx.c): Likewise.
	* sysdeps/x86_64/fpu/multiarch/e_exp-avx.c (__slowexp): Do not
	define as macro.
	* sysdeps/x86_64/fpu/multiarch/e_exp-fma.c (__slowexp): Likewise.
	* sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c (__slowexp): Likewise.
	* math/Makefile (type-double-routines): Remove slowexp.
	* manual/probes.texi (slowexp_p6): Remove.
	(slowexp_p32): Likewise.

19 files changed, 509 insertions, 326 deletions

diff --git a/ChangeLog b/ChangeLog
index 4b88a8601d..36d93bb30b 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,37 @@
+2017-12-19  Patrick McGehearty  <patrick.mcgehearty@oracle.com>
+
+	* sysdeps/ieee754/dbl-64/e_exp.c: Include <math-svid-compat.h> and
+	<errno.h>.  Include "eexp.tbl".
+	(half): New constant.
+	(one): Likewise.
+	(__ieee754_exp): Rewrite.
+	(__slowexp): Remove prototype.
+	* sysdeps/ieee754/dbl-64/eexp.tbl: New file.
+	* sysdeps/ieee754/dbl-64/slowexp.c: Remove file.
+	* sysdeps/i386/fpu/slowexp.c: Likewise.
+	* sysdeps/ia64/fpu/slowexp.c: Likewise.
+	* sysdeps/m68k/m680x0/fpu/slowexp.c: Likewise.
+	* sysdeps/x86_64/fpu/multiarch/slowexp-avx.c: Likewise.
+	* sysdeps/x86_64/fpu/multiarch/slowexp-fma.c: Likewise.
+	* sysdeps/x86_64/fpu/multiarch/slowexp-fma4.c: Likewise.
+	* sysdeps/generic/math_private.h (__slowexp): Remove prototype.
+	* sysdeps/ieee754/dbl-64/e_pow.c: Remove mention of slowexp.c in
+	comment.
+	* sysdeps/powerpc/power4/fpu/Makefile [$(subdir) = math]
+	(CPPFLAGS-slowexp.c): Remove variable.
+	* sysdeps/x86_64/fpu/multiarch/Makefile (libm-sysdep_routines):
+	Remove slowexp-fma, slowexp-fma4 and slowexp-avx.
+	(CFLAGS-slowexp-fma.c): Remove variable.
+	(CFLAGS-slowexp-fma4.c): Likewise.
+	(CFLAGS-slowexp-avx.c): Likewise.
+	* sysdeps/x86_64/fpu/multiarch/e_exp-avx.c (__slowexp): Do not
+	define as macro.
+	* sysdeps/x86_64/fpu/multiarch/e_exp-fma.c (__slowexp): Likewise.
+	* sysdeps/x86_64/fpu/multiarch/e_exp-fma4.c (__slowexp): Likewise.
+	* math/Makefile (type-double-routines): Remove slowexp.
+	* manual/probes.texi (slowexp_p6): Remove.
+	(slowexp_p32): Likewise.
+
 2017-12-19  Adhemerval Zanella  <adhemerval.zanella@linaro.org>
 	    James Clarke <jrtc27@jrtc27.com>
 
diff --git a/manual/probes.texi b/manual/probes.texi
index 8ab67562d7..f8ae64be33 100644
--- a/manual/probes.texi
+++ b/manual/probes.texi
@@ -258,20 +258,6 @@ Unless explicitly mentioned otherwise, a precision of 1 implies 24 bits of
 precision in the mantissa of the multiple precision number.  Hence, a precision
 level of 32 implies 768 bits of precision in the mantissa.
 
-@deftp Probe slowexp_p6 (double @var{$arg1}, double @var{$arg2})
-This probe is triggered when the @code{exp} function is called with an
-input that results in multiple precision computation with precision
-6.  Argument @var{$arg1} is the input value and @var{$arg2} is the
-computed output.
-@end deftp
-
-@deftp Probe slowexp_p32 (double @var{$arg1}, double @var{$arg2})
-This probe is triggered when the @code{exp} function is called with an
-input that results in multiple precision computation with precision
-32.  Argument @var{$arg1} is the input value and @var{$arg2} is the
-computed output.
-@end deftp
-
 @deftp Probe slowpow_p10 (double @var{$arg1}, double @var{$arg2}, double @var{$arg3}, double @var{$arg4})
 This probe is triggered when the @code{pow} function is called with
 inputs that result in multiple precision computation with precision
diff --git a/math/Makefile b/math/Makefile
index cba9ce9405..c8cd7cff85 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -114,7 +114,7 @@ type-ldouble-yes := ldouble
 # double support
 type-double-suffix :=
 type-double-routines := branred doasin dosincos halfulp mpa mpatan2	\
-		       mpatan mpexp mplog mpsqrt mptan sincos32 slowexp	\
+		       mpatan mpexp mplog mpsqrt mptan sincos32	\
 		       slowpow sincostab k_rem_pio2
 
 # float support
diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
index f29898c19c..689dc54767 100644
--- a/sysdeps/generic/math_private.h
+++ b/sysdeps/generic/math_private.h
@@ -262,7 +262,6 @@ 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, 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/i386/fpu/slowexp.c b/sysdeps/i386/fpu/slowexp.c
deleted file mode 100644
index 1cc8931700..0000000000
--- a/sysdeps/i386/fpu/slowexp.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
diff --git a/sysdeps/ia64/fpu/slowexp.c b/sysdeps/ia64/fpu/slowexp.c
deleted file mode 100644
index 1cc8931700..0000000000
--- a/sysdeps/ia64/fpu/slowexp.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Not needed.  */
diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c
index 6757a14ce1..6a7122f585 100644
--- a/sysdeps/ieee754/dbl-64/e_exp.c
+++ b/sysdeps/ieee754/dbl-64/e_exp.c
@@ -1,3 +1,4 @@
+/* EXP function - Compute double precision exponential */
 /*
  * IBM Accurate Mathematical Library
  * written by International Business Machines Corp.
@@ -23,7 +24,7 @@
 /*           exp1                                                          */
 /*                                                                         */
 /* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h                       */
-/*              mpa.c mpexp.x slowexp.c                                    */
+/*              mpa.c mpexp.x                                              */
 /*                                                                         */
 /* An ultimate exp routine. Given an IEEE double machine number x          */
 /* it computes the correctly rounded (to nearest) value of e^x             */
@@ -32,207 +33,238 @@
 /*                                                                         */
 /***************************************************************************/
 
+/*  IBM exp(x) replaced by following exp(x) in 2017. IBM exp1(x,xx) remains.  */
+/* exp(x)
+   Hybrid algorithm of Peter Tang's Table driven method (for large
+   arguments) and an accurate table (for small arguments).
+   Written by K.C. Ng, November 1988.
+   Revised by Patrick McGehearty, Nov 2017 to use j/64 instead of j/32
+   Method (large arguments):
+	1. Argument Reduction: given the input x, find r and integer k
+	   and j such that
+	             x = (k+j/64)*(ln2) + r,  |r| <= (1/128)*ln2
+
+	2. exp(x) = 2^k * (2^(j/64) + 2^(j/64)*expm1(r))
+	   a. expm1(r) is approximated by a polynomial:
+	      expm1(r) ~ r + t1*r^2 + t2*r^3 + ... + t5*r^6
+	      Here t1 = 1/2 exactly.
+	   b. 2^(j/64) is represented to twice double precision
+	      as TBL[2j]+TBL[2j+1].
+
+   Note: If divide were fast enough, we could use another approximation
+	 in 2.a:
+	      expm1(r) ~ (2r)/(2-R), R = r - r^2*(t1 + t2*r^2)
+	      (for the same t1 and t2 as above)
+
+   Special cases:
+	exp(INF) is INF, exp(NaN) is NaN;
+	exp(-INF)=  0;
+	for finite argument, only exp(0)=1 is exact.
+
+   Accuracy:
+	According to an error analysis, the error is always less than
+	an ulp (unit in the last place).  The largest errors observed
+	are less than 0.55 ulp for normal results and less than 0.75 ulp
+	for subnormal results.
+
+   Misc. info.
+	For IEEE double
+		if x >  7.09782712893383973096e+02 then exp(x) overflow
+		if x < -7.45133219101941108420e+02 then exp(x) underflow.  */
+
 #include <math.h>
+#include <math-svid-compat.h>
+#include <math_private.h>
+#include <errno.h>
 #include "endian.h"
 #include "uexp.h"
+#include "uexp.tbl"
 #include "mydefs.h"
 #include "MathLib.h"
-#include "uexp.tbl"
-#include <math_private.h>
 #include <fenv.h>
 #include <float.h>
 
-#ifndef SECTION
-# define SECTION
-#endif
+extern double __ieee754_exp (double);
+
+#include "eexp.tbl"
+
+static const double
+  half = 0.5,
+  one = 1.0;
 
-double __slowexp (double);
 
-/* An ultimate exp routine. Given an IEEE double machine number x it computes
-   the correctly rounded (to nearest) value of e^x.  */
 double
-SECTION
-__ieee754_exp (double x)
+__ieee754_exp (double x_arg)
 {
-  double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
-  mynumber junk1, junk2, binexp = {{0, 0}};
-  int4 i, j, m, n, ex;
+  double z, t;
   double retval;
-
+  int hx, ix, k, j, m;
+  int fe_val;
+  union
   {
-    SET_RESTORE_ROUND (FE_TONEAREST);
-
-    junk1.x = x;
-    m = junk1.i[HIGH_HALF];
-    n = m & hugeint;
-
-    if (n > smallint && n < bigint)
-      {
-	y = x * log2e.x + three51.x;
-	bexp = y - three51.x;	/*  multiply the result by 2**bexp        */
-
-	junk1.x = y;
-
-	eps = bexp * ln_two2.x;	/* x = bexp*ln(2) + t - eps               */
-	t = x - bexp * ln_two1.x;
-
-	y = t + three33.x;
-	base = y - three33.x;	/* t rounded to a multiple of 2**-18      */
-	junk2.x = y;
-	del = (t - base) - eps;	/*  x = bexp*ln(2) + base + del           */
-	eps = del + del * del * (p3.x * del + p2.x);
-
-	binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 1023) << 20;
-
-	i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
-	j = (junk2.i[LOW_HALF] & 511) << 1;
-
-	al = coar.x[i] * fine.x[j];
-	bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
-	       + coar.x[i + 1] * fine.x[j + 1]);
-
-	rem = (bet + bet * eps) + al * eps;
-	res = al + rem;
-	cor = (al - res) + rem;
-	if (res == (res + cor * err_0))
-	  {
-	    retval = res * binexp.x;
-	    goto ret;
-	  }
-	else
-	  {
-	    retval = __slowexp (x);
-	    goto ret;
-	  }			/*if error is over bound */
-      }
-
-    if (n <= smallint)
-      {
-	retval = 1.0;
-	goto ret;
-      }
-
-    if (n >= badint)
-      {
-	if (n > infint)
-	  {
-	    retval = x + x;
-	    goto ret;
-	  }			/* x is NaN */
-	if (n < infint)
-	  {
-	    if (x > 0)
-	      goto ret_huge;
-	    else
-	      goto ret_tiny;
-	  }
-	/* x is finite,  cause either overflow or underflow  */
-	if (junk1.i[LOW_HALF] != 0)
-	  {
-	    retval = x + x;
-	    goto ret;
-	  }			/*  x is NaN  */
-	retval = (x > 0) ? inf.x : zero;	/* |x| = inf;  return either inf or 0 */
-	goto ret;
-      }
-
-    y = x * log2e.x + three51.x;
-    bexp = y - three51.x;
-    junk1.x = y;
-    eps = bexp * ln_two2.x;
-    t = x - bexp * ln_two1.x;
-    y = t + three33.x;
-    base = y - three33.x;
-    junk2.x = y;
-    del = (t - base) - eps;
-    eps = del + del * del * (p3.x * del + p2.x);
-    i = ((junk2.i[LOW_HALF] >> 8) & 0xfffffffe) + 356;
-    j = (junk2.i[LOW_HALF] & 511) << 1;
-    al = coar.x[i] * fine.x[j];
-    bet = ((coar.x[i] * fine.x[j + 1] + coar.x[i + 1] * fine.x[j])
-	   + coar.x[i + 1] * fine.x[j + 1]);
-    rem = (bet + bet * eps) + al * eps;
-    res = al + rem;
-    cor = (al - res) + rem;
-    if (m >> 31)
-      {
-	ex = junk1.i[LOW_HALF];
-	if (res < 1.0)
-	  {
-	    res += res;
-	    cor += cor;
-	    ex -= 1;
-	  }
-	if (ex >= -1022)
-	  {
-	    binexp.i[HIGH_HALF] = (1023 + ex) << 20;
-	    if (res == (res + cor * err_0))
-	      {
-		retval = res * binexp.x;
-		goto ret;
-	      }
-	    else
-	      {
-		retval = __slowexp (x);
-		goto check_uflow_ret;
-	      }			/*if error is over bound */
-	  }
-	ex = -(1022 + ex);
-	binexp.i[HIGH_HALF] = (1023 - ex) << 20;
-	res *= binexp.x;
-	cor *= binexp.x;
-	eps = 1.0000000001 + err_0 * binexp.x;
-	t = 1.0 + res;
-	y = ((1.0 - t) + res) + cor;
-	res = t + y;
-	cor = (t - res) + y;
-	if (res == (res + eps * cor))
-	  {
-	    binexp.i[HIGH_HALF] = 0x00100000;
-	    retval = (res - 1.0) * binexp.x;
-	    goto check_uflow_ret;
-	  }
-	else
-	  {
-	    retval = __slowexp (x);
-	    goto check_uflow_ret;
-	  }			/*   if error is over bound    */
-      check_uflow_ret:
-	if (retval < DBL_MIN)
-	  {
-	    double force_underflow = tiny * tiny;
-	    math_force_eval (force_underflow);
-	  }
-	if (retval == 0)
-	  goto ret_tiny;
-	goto ret;
-      }
-    else
-      {
-	binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
-	if (res == (res + cor * err_0))
-	  retval = res * binexp.x * t256.x;
-	else
-	  retval = __slowexp (x);
-	if (isinf (retval))
-	  goto ret_huge;
-	else
-	  goto ret;
-      }
-  }
-ret:
-  return retval;
-
- ret_huge:
-  return hhuge * hhuge;
-
- ret_tiny:
-  return tiny * tiny;
+    int i_part[2];
+    double x;
+  } xx;
+  union
+  {
+    int y_part[2];
+    double y;
+  } yy;
+  xx.x = x_arg;
+
+  ix = xx.i_part[HIGH_HALF];
+  hx = ix & ~0x80000000;
+
+  if (hx < 0x3ff0a2b2)
+    {				/* |x| < 3/2 ln 2 */
+      if (hx < 0x3f862e42)
+	{			/* |x| < 1/64 ln 2 */
+	  if (hx < 0x3ed00000)
+	    {			/* |x| < 2^-18 */
+	      if (hx < 0x3e300000)
+		{
+		  retval = one + xx.x;
+		  return retval;
+		}
+	      retval = one + xx.x * (one + half * xx.x);
+	      return retval;
+	    }
+	  /* Use FE_TONEAREST rounding mode for computing yy.y.
+	     Avoid set/reset of rounding mode if in FE_TONEAREST mode.  */
+	  fe_val = get_rounding_mode ();
+	  if (fe_val == FE_TONEAREST)
+	    {
+	      t = xx.x * xx.x;
+	      yy.y = xx.x + (t * (half + xx.x * t2)
+			     + (t * t) * (t3 + xx.x * t4 + t * t5));
+	      retval = one + yy.y;
+	    }
+	  else
+	    {
+	      libc_fesetround (FE_TONEAREST);
+	      t = xx.x * xx.x;
+	      yy.y = xx.x + (t * (half + xx.x * t2)
+			     + (t * t) * (t3 + xx.x * t4 + t * t5));
+	      retval = one + yy.y;
+	      libc_fesetround (fe_val);
+	    }
+	  return retval;
+	}
+
+      /* Find the multiple of 2^-6 nearest x.  */
+      k = hx >> 20;
+      j = (0x00100000 | (hx & 0x000fffff)) >> (0x40c - k);
+      j = (j - 1) & ~1;
+      if (ix < 0)
+	j += 134;
+      /* Use FE_TONEAREST rounding mode for computing yy.y.
+	 Avoid set/reset of rounding mode if in FE_TONEAREST mode.  */
+      fe_val = get_rounding_mode ();
+      if (fe_val == FE_TONEAREST)
+	{
+	  z = xx.x - TBL2[j];
+	  t = z * z;
+	  yy.y = z + (t * (half + (z * t2))
+		      + (t * t) * (t3 + z * t4 + t * t5));
+	  retval = TBL2[j + 1] + TBL2[j + 1] * yy.y;
+	}
+      else
+	{
+	  libc_fesetround (FE_TONEAREST);
+	  z = xx.x - TBL2[j];
+	  t = z * z;
+	  yy.y = z + (t * (half + (z * t2))
+		      + (t * t) * (t3 + z * t4 + t * t5));
+	  retval = TBL2[j + 1] + TBL2[j + 1] * yy.y;
+	  libc_fesetround (fe_val);
+	}
+      return retval;
+    }
+
+  if (hx >= 0x40862e42)
+    {				/* x is large, infinite, or nan.  */
+      if (hx >= 0x7ff00000)
+	{
+	  if (ix == 0xfff00000 && xx.i_part[LOW_HALF] == 0)
+	    return zero;	/* exp(-inf) = 0.  */
+	  return (xx.x * xx.x);	/* exp(nan/inf) is nan or inf.  */
+	}
+      if (xx.x > threshold1)
+	{			/* Set overflow error condition.  */
+	  retval = hhuge * hhuge;
+	  return retval;
+	}
+      if (-xx.x > threshold2)
+	{			/* Set underflow error condition.  */
+	  double force_underflow = tiny * tiny;
+	  math_force_eval (force_underflow);
+	  retval = force_underflow;
+	  return retval;
+	}
+    }
+
+  /* Use FE_TONEAREST rounding mode for computing yy.y.
+     Avoid set/reset of rounding mode if already in FE_TONEAREST mode.  */
+  fe_val = get_rounding_mode ();
+  if (fe_val == FE_TONEAREST)
+    {
+      t = invln2_64 * xx.x;
+      if (ix < 0)
+	t -= half;
+      else
+	t += half;
+      k = (int) t;
+      j = (k & 0x3f) << 1;
+      m = k >> 6;
+      z = (xx.x - k * ln2_64hi) - k * ln2_64lo;
+
+      /* z is now in primary range.  */
+      t = z * z;
+      yy.y = z + (t * (half + z * t2) + (t * t) * (t3 + z * t4 + t * t5));
+      yy.y = TBL[j] + (TBL[j + 1] + TBL[j] * yy.y);
+    }
+  else
+    {
+      libc_fesetround (FE_TONEAREST);
+      t = invln2_64 * xx.x;
+      if (ix <