Add new exp and exp2 implementations

Optimized exp and exp2 implementations using a lookup table for
fractional powers of 2.  There are several variants, see e_exp_data.c,
they can be selected by modifying math_config.h allowing different
tradeoffs.

The default selection should be acceptable as generic libm code.
Worst case error is 0.509 ULP for exp and 0.507 ULP for exp2, on
aarch64 the rodata size is 2160 bytes, shared between exp and exp2.
On aarch64 .text + .rodata size decreased by 24912 bytes.

The non-nearest rounding error is less than 1 ULP even on targets
without efficient round implementation (although the error rate is
higher in that case).  Targets with single instruction, rounding mode
independent, to nearest integer rounding and conversion can use them
by setting TOINT_INTRINSICS and adding the necessary code to their
math_private.h.

The __exp1 code uses the same algorithm, so the error bound of pow
increased a bit.

New double precision error handling code was added following the
style of the single precision error handling code.

Improvements on Cortex-A72 compared to current glibc master:
exp thruput: 1.61x in [-9.9 9.9]
exp latency: 1.53x in [-9.9 9.9]
exp thruput: 1.13x in [0.5 1]
exp latency: 1.30x in [0.5 1]
exp2 thruput: 2.03x in [-9.9 9.9]
exp2 latency: 1.64x in [-9.9 9.9]

For small (< 1) inputs the current exp code uses a separate algorithm
so the speed up there is less.

Was tested on
aarch64-linux-gnu (TOINT_INTRINSICS, fma contraction) and
arm-linux-gnueabihf (!TOINT_INTRINSICS, no fma contraction) and
x86_64-linux-gnu (!TOINT_INTRINSICS, no fma contraction) and
powerpc64le-linux-gnu (!TOINT_INTRINSICS, fma contraction) targets,
only non-nearest rounding ulp errors increase and they are within
acceptable bounds (ulp updates are in separate patches).

	* NEWS: Mention exp and exp2 improvements.
	* math/Makefile (libm-support): Remove t_exp.
	(type-double-routines): Add math_err and e_exp_data.
	* sysdeps/aarch64/libm-test-ulps: Update.
	* sysdeps/arm/libm-test-ulps: Update.
	* sysdeps/i386/fpu/e_exp_data.c: New file.
	* sysdeps/i386/fpu/math_err.c: New file.
	* sysdeps/i386/fpu/t_exp.c: Remove.
	* sysdeps/ia64/fpu/e_exp_data.c: New file.
	* sysdeps/ia64/fpu/math_err.c: New file.
	* sysdeps/ia64/fpu/t_exp.c: Remove.
	* sysdeps/ieee754/dbl-64/e_exp.c: Rewrite.
	* sysdeps/ieee754/dbl-64/e_exp2.c: Rewrite.
	* sysdeps/ieee754/dbl-64/e_exp_data.c: New file.
	* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Update error bound.
	* sysdeps/ieee754/dbl-64/eexp.tbl: Remove.
	* sysdeps/ieee754/dbl-64/math_config.h: New file.
	* sysdeps/ieee754/dbl-64/math_err.c: New file.
	* sysdeps/ieee754/dbl-64/t_exp.c: Remove.
	* sysdeps/ieee754/dbl-64/t_exp2.h: Remove.
	* sysdeps/ieee754/dbl-64/uexp.h: Remove.
	* sysdeps/ieee754/dbl-64/uexp.tbl: Remove.
	* sysdeps/m68k/m680x0/fpu/e_exp_data.c: New file.
	* sysdeps/m68k/m680x0/fpu/math_err.c: New file.
	* sysdeps/m68k/m680x0/fpu/t_exp.c: Remove.
	* sysdeps/powerpc/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Update.

1 files changed, 154 insertions, 331 deletions

diff --git a/sysdeps/ieee754/dbl-64/e_exp.c b/sysdeps/ieee754/dbl-64/e_exp.c
index 7d8b414034..209f20b972 100644
--- a/sysdeps/ieee754/dbl-64/e_exp.c
+++ b/sysdeps/ieee754/dbl-64/e_exp.c
@@ -1,48 +1,158 @@
-/*
- * IBM Accurate Mathematical Library
- * written by International Business Machines Corp.
- * Copyright (C) 2001-2018 Free Software Foundation, Inc.
- *
- * This program is free software; you can redistribute it and/or modify
- * it under the terms of the GNU Lesser General Public License as published by
- * the Free Software Foundation; either version 2.1 of the License, or
- * (at your option) any later version.
- *
- * This program is distributed in the hope that it will be useful,
- * but WITHOUT ANY WARRANTY; without even the implied warranty of
- * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
- * GNU Lesser General Public License for more details.
- *
- * You should have received a copy of the GNU Lesser General Public License
- * along with this program; if not, see <http://www.gnu.org/licenses/>.
- */
-/***************************************************************************/
-/*  MODULE_NAME:uexp.c                                                     */
-/*                                                                         */
-/*  FUNCTION:uexp                                                          */
-/*           exp1                                                          */
-/*                                                                         */
-/* FILES NEEDED:dla.h endian.h mpa.h mydefs.h uexp.h                       */
-/*                                                                         */
-/* An ultimate exp routine. Given an IEEE double machine number x          */
-/* it computes an almost correctly rounded (to nearest) value of e^x       */
-/* Assumption: Machine arithmetic operations are performed in              */
-/* round to nearest mode of IEEE 754 standard.                             */
-/*                                                                         */
-/***************************************************************************/
+/* Double-precision e^x function.
+   Copyright (C) 2018 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+
+   The GNU C Library is free software; you can redistribute it and/or
+   modify it under the terms of the GNU Lesser General Public
+   License as published by the Free Software Foundation; either
+   version 2.1 of the License, or (at your option) any later version.
+
+   The GNU C Library is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
+   Lesser General Public License for more details.
+
+   You should have received a copy of the GNU Lesser General Public
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
 
 #include <math.h>
-#include "endian.h"
-#include "uexp.h"
-#include "mydefs.h"
-#include "MathLib.h"
-#include "uexp.tbl"
+#include <stdint.h>
 #include <math-barriers.h>
-#include <math_private.h>
-#include <fenv_private.h>
-#include <fenv.h>
-#include <float.h>
-#include "eexp.tbl"
+#include <math-narrow-eval.h>
+#include "math_config.h"
+
+#define N (1 << EXP_TABLE_BITS)
+#define InvLn2N __exp_data.invln2N
+#define NegLn2hiN __exp_data.negln2hiN
+#define NegLn2loN __exp_data.negln2loN
+#define Shift __exp_data.shift
+#define T __exp_data.tab
+#define C2 __exp_data.poly[5 - EXP_POLY_ORDER]
+#define C3 __exp_data.poly[6 - EXP_POLY_ORDER]
+#define C4 __exp_data.poly[7 - EXP_POLY_ORDER]
+#define C5 __exp_data.poly[8 - EXP_POLY_ORDER]
+
+/* Handle cases that may overflow or underflow when computing the result that
+   is scale*(1+TMP) without intermediate rounding.  The bit representation of
+   scale is in SBITS, however it has a computed exponent that may have
+   overflown into the sign bit so that needs to be adjusted before using it as
+   a double.  (int32_t)KI is the k used in the argument reduction and exponent
+   adjustment of scale, positive k here means the result may overflow and
+   negative k means the result may underflow.  */
+static inline double
+specialcase (double_t tmp, uint64_t sbits, uint64_t ki)
+{
+  double_t scale, y;
+
+  if ((ki & 0x80000000) == 0)
+    {
+      /* k > 0, the exponent of scale might have overflowed by <= 460.  */
+      sbits -= 1009ull << 52;
+      scale = asdouble (sbits);
+      y = 0x1p1009 * (scale + scale * tmp);
+      return check_oflow (y);
+    }
+  /* k < 0, need special care in the subnormal range.  */
+  sbits += 1022ull << 52;
+  scale = asdouble (sbits);
+  y = scale + scale * tmp;
+  if (y < 1.0)
+    {
+      /* Round y to the right precision before scaling it into the subnormal
+	 range to avoid double rounding that can cause 0.5+E/2 ulp error where
+	 E is the worst-case ulp error outside the subnormal range.  So this
+	 is only useful if the goal is better than 1 ulp worst-case error.  */
+      double_t hi, lo;
+      lo = scale - y + scale * tmp;
+      hi = 1.0 + y;
+      lo = 1.0 - hi + y + lo;
+      y = math_narrow_eval (hi + lo) - 1.0;
+      /* Avoid -0.0 with downward rounding.  */
+      if (WANT_ROUNDING && y == 0.0)
+	y = 0.0;
+      /* The underflow exception needs to be signaled explicitly.  */
+      math_force_eval (math_opt_barrier (0x1p-1022) * 0x1p-1022);
+    }
+  y = 0x1p-1022 * y;
+  return check_uflow (y);
+}
+
+/* Top 12 bits of a double (sign and exponent bits).  */
+static inline uint32_t
+top12 (double x)
+{
+  return asuint64 (x) >> 52;
+}
+
+/* Computes exp(x+xtail) where |xtail| < 2^-8/N and |xtail| <= |x|.
+   If hastail is 0 then xtail is assumed to be 0 too.  */
+static inline double
+exp_inline (double x, double xtail, int hastail)
+{
+  uint32_t abstop;
+  uint64_t ki, idx, top, sbits;
+  /* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
+  double_t kd, z, r, r2, scale, tail, tmp;
+
+  abstop = top12 (x) & 0x7ff;
+  if (__glibc_unlikely (abstop - top12 (0x1p-54)
+			>= top12 (512.0) - top12 (0x1p-54)))
+    {
+      if (abstop - top12 (0x1p-54) >= 0x80000000)
+	/* Avoid spurious underflow for tiny x.  */
+	/* Note: 0 is common input.  */
+	return WANT_ROUNDING ? 1.0 + x : 1.0;
+      if (abstop >= top12 (1024.0))
+	{
+	  if (asuint64 (x) == asuint64 (-INFINITY))
+	    return 0.0;
+	  if (abstop >= top12 (INFINITY))
+	    return 1.0 + x;
+	  if (asuint64 (x) >> 63)
+	    return __math_uflow (0);
+	  else
+	    return __math_oflow (0);
+	}
+      /* Large x is special cased below.  */
+      abstop = 0;
+    }
+
+  /* exp(x) = 2^(k/N) * exp(r), with exp(r) in [2^(-1/2N),2^(1/2N)].  */
+  /* x = ln2/N*k + r, with int k and r in [-ln2/2N, ln2/2N].  */
+  z = InvLn2N * x;
+#if TOINT_INTRINSICS
+  kd = roundtoint (z);
+  ki = converttoint (z);
+#else
+  /* z - kd is in [-1, 1] in non-nearest rounding modes.  */
+  kd = math_narrow_eval (z + Shift);
+  ki = asuint64 (kd);
+  kd -= Shift;
+#endif
+  r = x + kd * NegLn2hiN + kd * NegLn2loN;
+  /* The code assumes 2^-200 < |xtail| < 2^-8/N.  */
+  if (hastail)
+    r += xtail;
+  /* 2^(k/N) ~= scale * (1 + tail).  */
+  idx = 2 * (ki % N);
+  top = ki << (52 - EXP_TABLE_BITS);
+  tail = asdouble (T[idx]);
+  /* This is only a valid scale when -1023*N < k < 1024*N.  */
+  sbits = T[idx + 1] + top;
+  /* exp(x) = 2^(k/N) * exp(r) ~= scale + scale * (tail + exp(r) - 1).  */
+  /* Evaluation is optimized assuming superscalar pipelined execution.  */
+  r2 = r * r;
+  /* Without fma the worst case error is 0.25/N ulp larger.  */
+  /* Worst case error is less than 0.5+1.11/N+(abs poly error * 2^53) ulp.  */
+  tmp = tail + r + r2 * (C2 + r * C3) + r2 * r2 * (C4 + r * C5);
+  if (__glibc_unlikely (abstop == 0))
+    return specialcase (tmp, sbits, ki);
+  scale = asdouble (sbits);
+  /* Note: tmp == 0 or |tmp| > 2^-200 and scale > 2^-739, so there
+     is no spurious underflow here even without fma.  */
+  return scale + scale * tmp;
+}
 
 #ifndef SECTION
 # define SECTION
@@ -52,187 +162,7 @@ double
 SECTION
 __ieee754_exp (double x)
 {
-  double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
-  double z;
-  mynumber junk1, junk2, binexp = {{0, 0}};
-  int4 i, j, m, n, ex;
-  int4 k;
-  double retval;
-
-  {
-    SET_RESTORE_ROUND (FE_TONEAREST);
-
-    junk1.x = x;
-    m = junk1.i[HIGH_HALF];
-    n = m & hugeint;
-
-    if (n < 0x3ff0a2b2)		/* |x| < 1.03972053527832 */
-      {
-	if (n < 0x3f862e42)	/* |x| < 3/2 ln 2 */
-	  {
-	    if (n < 0x3ed00000)	/* |x| < 1/64 ln 2 */
-	      {
-		if (n < 0x3e300000)	/* |x| < 2^18 */
-		  {
-		    retval = one + junk1.x;
-		    goto ret;
-		  }
-		retval = one + junk1.x * (one + half * junk1.x);
-		goto ret;
-	      }
-	    t = junk1.x * junk1.x;
-	    retval = junk1.x + (t * (half + junk1.x * t2) +
-				(t * t) * (t3 + junk1.x * t4 + t * t5));
-	    retval = one + retval;
-	    goto ret;
-	  }
-
-	/* Find the multiple of 2^-6 nearest x.  */
-	k = n >> 20;
-	j = (0x00100000 | (n & 0x000fffff)) >> (0x40c - k);
-	j = (j - 1) & ~1;
-	if (m < 0)
-	  j += 134;
-	z = junk1.x - TBL2[j];
-	t = z * z;
-	retval = z + (t * (half + (z * t2))
-		      + (t * t) * (t3 + z * t4 + t * t5));
-	retval = TBL2[j + 1] + TBL2[j + 1] * retval;
-	goto ret;
-      }
-
-    if (n < bigint)		/* && |x| >= 1.03972053527832 */
-      {
-	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;
-	/* Maximum relative error is 7.8e-22 (70.1 bits).
-	   Maximum ULP error is 0.500007.  */
-	retval = res * binexp.x;
-	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;
-	    /* Does not underflow: res >= 1.0, binexp >= 0x1p-1022
-	       Maximum relative error is 7.8e-22 (70.1 bits).
-	       Maximum ULP error is 0.500007.  */
-	    retval = res * binexp.x;
-	    goto ret;
-	  }
-	ex = -(1022 + ex);
-	binexp.i[HIGH_HALF] = (1023 - ex) << 20;
-	res *= binexp.x;
-	cor *= binexp.x;
-	t = 1.0 + res;
-	y = ((1.0 - t) + res) + cor;
-	res = t + y;
-	/* Maximum ULP error is 0.5000035.  */
-	binexp.i[HIGH_HALF] = 0x00100000;
-	retval = (res - 1.0) * binexp.x;
-	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;
-	/* Maximum relative error is 7.8e-22 (70.1 bits).
-	   Maximum ULP error is 0.500007.  */
-	retval = res * binexp.x * t256.x;
-	if (isinf (retval))
-	  goto ret_huge;
-	else
-	  goto ret;
-      }
-  }
-ret:
-  return retval;
-
- ret_huge:
-  return hhuge * hhuge;
-
- ret_tiny:
-  return tiny * tiny;
+  return exp_inline (x, 0, 0);
 }
 #ifndef __ieee754_exp
 strong_alias (__ieee754_exp, __exp_finite)
@@ -243,112 +173,5 @@ double
 SECTION
 __exp1 (double x, double xx)
 {
-  double bexp, t, eps, del, base, y, al, bet, res, rem, cor;
-  mynumber junk1, junk2, binexp = {{0, 0}};
-  int4 i, j, m, n, ex;
-
-  junk1.x = x;
-  m = junk1.i[HIGH_HALF];
-  n = m & hugeint;		/* no sign */
-
-  /* fabs (x) > 5.551112e-17 and fabs (x) < 7.080010e+02.  */
-  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) + (xx - 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;
-      /* Maximum relative error before rounding is 8.8e-22 (69.9 bits).
-	 Maximum ULP error is 0.500008.  */
-      return res * binexp.x;
-    }
-
-  if (n <= smallint)
-    return 1.0;			/*  if x->0 e^x=1 */
-
-  if (n >= badint)
-    {
-      if (n > infint)
-	return (zero / zero);	/* x is NaN,  return invalid */
-      if (n < infint)
-	return ((x > 0) ? (hhuge * hhuge) : (tiny * tiny));
-      /* x is finite,  cause either overflow or underflow  */
-      if (junk1.i[LOW_HALF] != 0)
-	return (zero / zero);	/*  x is NaN  */
-      return ((x > 0) ? inf.x : zero);	/* |x| = inf;  return either inf or 0 */
-    }
-
-  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) + (xx - 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)
-    {
-      /* x < 0.  */
-      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;
-	  /* Maximum ULP error is 0.500008.  */
-	  return res * binexp.x;
-	}
-      /* Denormal case - ex < -1022.  */
-      ex = -(1022 + ex);
-      binexp.i[HIGH_HALF] = (1023 - ex) << 20;
-      res *= binexp.x;
-      cor *= binexp.x;
-      t = 1.0 + res;
-      y = ((1.0 - t) + res) + cor;
-      res = t + y;
-      binexp.i[HIGH_HALF] = 0x00100000;
-      /* Maximum ULP error is 0.500004.  */
-      return (res - 1.0) * binexp.x;
-    }
-  else
-    {
-      binexp.i[HIGH_HALF] = (junk1.i[LOW_HALF] + 767) << 20;
-      /* Maximum ULP error is 0.500008.  */
-      return res * binexp.x * t256.x;
-    }
+  return exp_inline (x, xx, 1);
 }