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.

26 files changed, 885 insertions, 3672 deletions

diff --git a/NEWS b/NEWS
index 325157c0da..085325ab87 100644
--- a/NEWS
+++ b/NEWS
@@ -16,7 +16,7 @@ Major new features:
   to set the install root if you wish to install into a non-default
   configured location.
 
-* Optimized generic sinf, cosf, sincosf and tanf.
+* Optimized generic exp, exp2, sinf, cosf, sincosf and tanf.
 
 * The reallocarray function is now declared under _DEFAULT_SOURCE, not just
   for _GNU_SOURCE, to match BSD environments.
diff --git a/math/Makefile b/math/Makefile
index bff726fe81..f1ba1a0c36 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -42,7 +42,7 @@ extra-libs-others = $(extra-libs)
 libm-support = s_lib_version s_matherr s_signgam			\
 	       fclrexcpt fgetexcptflg fraiseexcpt fsetexcptflg		\
 	       ftestexcept fegetround fesetround fegetenv feholdexcpt	\
-	       fesetenv feupdateenv t_exp fedisblxcpt feenablxcpt	\
+	       fesetenv feupdateenv fedisblxcpt feenablxcpt	\
 	       fegetexcept fesetexcept fetestexceptflag fegetmode	\
 	       fesetmode
 
@@ -127,7 +127,7 @@ type-ldouble-yes := ldouble
 type-double-suffix :=
 type-double-routines := branred doasin dosincos mpa mpatan2	\
 		       k_rem_pio2 mpatan mpsqrt mptan sincos32	\
-		       sincostab
+		       sincostab math_err e_exp_data
 
 # float support
 type-float-suffix := f
diff --git a/sysdeps/aarch64/libm-test-ulps b/sysdeps/aarch64/libm-test-ulps
index 6a89028459..585e5bbce7 100644
--- a/sysdeps/aarch64/libm-test-ulps
+++ b/sysdeps/aarch64/libm-test-ulps
@@ -716,9 +716,9 @@ ildouble: 2
 ldouble: 2
 
 Function: Imaginary part of "ccos_downward":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
@@ -732,9 +732,9 @@ ildouble: 2
 ldouble: 2
 
 Function: Imaginary part of "ccos_towardzero":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
@@ -772,33 +772,33 @@ ildouble: 1
 ldouble: 1
 
 Function: Real part of "ccosh_downward":
-double: 1
+double: 2
 float: 3
-idouble: 1
+idouble: 2
 ifloat: 3
 ildouble: 2
 ldouble: 2
 
 Function: Imaginary part of "ccosh_downward":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
 
 Function: Real part of "ccosh_towardzero":
-double: 1
+double: 2
 float: 3
-idouble: 1
+idouble: 2
 ifloat: 3
 ildouble: 2
 ldouble: 2
 
 Function: Imaginary part of "ccosh_towardzero":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
@@ -836,33 +836,33 @@ ildouble: 1
 ldouble: 1
 
 Function: Real part of "cexp_downward":
-double: 1
+double: 2
 float: 2
-idouble: 1
+idouble: 2
 ifloat: 2
 ildouble: 2
 ldouble: 2
 
 Function: Imaginary part of "cexp_downward":
-double: 1
+double: 3
 float: 3
-idouble: 1
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
 
 Function: Real part of "cexp_towardzero":
-double: 1
+double: 2
 float: 2
-idouble: 1
+idouble: 2
 ifloat: 2
 ildouble: 2
 ldouble: 2
 
 Function: Imaginary part of "cexp_towardzero":
-double: 1
+double: 3
 float: 3
-idouble: 1
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
@@ -876,9 +876,9 @@ ildouble: 3
 ldouble: 3
 
 Function: Imaginary part of "cexp_upward":
-double: 1
+double: 3
 float: 2
-idouble: 1
+idouble: 3
 ifloat: 2
 ildouble: 3
 ldouble: 3
@@ -1052,25 +1052,25 @@ ildouble: 1
 ldouble: 1
 
 Function: "cosh_downward":
-double: 1
+double: 2
 float: 1
-idouble: 1
+idouble: 2
 ifloat: 1
 ildouble: 1
 ldouble: 2
 
 Function: "cosh_towardzero":
-double: 1
+double: 2
 float: 1
-idouble: 1
+idouble: 2
 ifloat: 1
 ildouble: 1
 ldouble: 2
 
 Function: "cosh_upward":
-double: 1
+double: 2
 float: 2
-idouble: 1
+idouble: 2
 ifloat: 2
 ildouble: 1
 ldouble: 3
@@ -1090,9 +1090,9 @@ ildouble: 1
 ldouble: 1
 
 Function: Real part of "cpow_downward":
-double: 4
+double: 5
 float: 8
-idouble: 4
+idouble: 5
 ifloat: 8
 ildouble: 6
 ldouble: 6
@@ -1106,9 +1106,9 @@ ildouble: 2
 ldouble: 2
 
 Function: Real part of "cpow_towardzero":
-double: 4
+double: 5
 float: 8
-idouble: 4
+idouble: 5
 ifloat: 8
 ildouble: 6
 ldouble: 6
@@ -1150,9 +1150,9 @@ ildouble: 1
 ldouble: 1
 
 Function: Real part of "csin_downward":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
@@ -1166,9 +1166,9 @@ ildouble: 2
 ldouble: 2
 
 Function: Real part of "csin_towardzero":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
@@ -1220,9 +1220,9 @@ ildouble: 2
 ldouble: 2
 
 Function: Imaginary part of "csinh_downward":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
@@ -1236,9 +1236,9 @@ ildouble: 2
 ldouble: 2
 
 Function: Imaginary part of "csinh_towardzero":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 ildouble: 2
 ldouble: 2
@@ -1492,9 +1492,9 @@ ildouble: 2
 ldouble: 2
 
 Function: "erfc_downward":
-double: 3
+double: 4
 float: 4
-idouble: 3
+idouble: 4
 ifloat: 4
 ildouble: 5
 ldouble: 5
@@ -1508,9 +1508,9 @@ ildouble: 4
 ldouble: 4
 
 Function: "erfc_upward":
-double: 3
+double: 4
 float: 4
-idouble: 3
+idouble: 4
 ifloat: 4
 ildouble: 5
 ldouble: 5
diff --git a/sysdeps/arm/libm-test-ulps b/sysdeps/arm/libm-test-ulps
index e62cca3e03..798e3c465c 100644
--- a/sysdeps/arm/libm-test-ulps
+++ b/sysdeps/arm/libm-test-ulps
@@ -530,9 +530,9 @@ idouble: 1
 ifloat: 1
 
 Function: Imaginary part of "ccos_downward":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 
 Function: Real part of "ccos_towardzero":
@@ -542,9 +542,9 @@ idouble: 1
 ifloat: 2
 
 Function: Imaginary part of "ccos_towardzero":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 
 Function: Real part of "ccos_upward":
@@ -572,27 +572,27 @@ idouble: 1
 ifloat: 1
 
 Function: Real part of "ccosh_downward":
-double: 1
+double: 2
 float: 3
-idouble: 1
+idouble: 2
 ifloat: 3
 
 Function: Imaginary part of "ccosh_downward":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 
 Function: Real part of "ccosh_towardzero":
-double: 1
+double: 2
 float: 3
-idouble: 1
+idouble: 2
 ifloat: 3
 
 Function: Imaginary part of "ccosh_towardzero":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 
 Function: Real part of "ccosh_upward":
@@ -620,27 +620,27 @@ idouble: 1
 ifloat: 2
 
 Function: Real part of "cexp_downward":
-double: 1
+double: 2
 float: 2
-idouble: 1
+idouble: 2
 ifloat: 2
 
 Function: Imaginary part of "cexp_downward":
-double: 1
+double: 3
 float: 3
-idouble: 1
+idouble: 3
 ifloat: 3
 
 Function: Real part of "cexp_towardzero":
-double: 1
+double: 2
 float: 2
-idouble: 1
+idouble: 2
 ifloat: 2
 
 Function: Imaginary part of "cexp_towardzero":
-double: 1
+double: 3
 float: 3
-idouble: 1
+idouble: 3
 ifloat: 3
 
 Function: Real part of "cexp_upward":
@@ -650,9 +650,9 @@ idouble: 1
 ifloat: 2
 
 Function: Imaginary part of "cexp_upward":
-double: 1
+double: 3
 float: 2
-idouble: 1
+idouble: 3
 ifloat: 2
 
 Function: Real part of "clog":
@@ -780,21 +780,21 @@ idouble: 1
 ifloat: 1
 
 Function: "cosh_downward":
-double: 1
+double: 2
 float: 1
-idouble: 1
+idouble: 2
 ifloat: 1
 
 Function: "cosh_towardzero":
-double: 1
+double: 2
 float: 1
-idouble: 1
+idouble: 2
 ifloat: 1
 
 Function: "cosh_upward":
-double: 1
+double: 2
 float: 2
-idouble: 1
+idouble: 2
 ifloat: 2
 
 Function: Real part of "cpow":
@@ -808,9 +808,9 @@ float: 2
 ifloat: 2
 
 Function: Real part of "cpow_downward":
-double: 4
+double: 5
 float: 8
-idouble: 4
+idouble: 5
 ifloat: 8
 
 Function: Imaginary part of "cpow_downward":
@@ -820,9 +820,9 @@ idouble: 1
 ifloat: 2
 
 Function: Real part of "cpow_towardzero":
-double: 4
+double: 5
 float: 8
-idouble: 4
+idouble: 5
 ifloat: 8
 
 Function: Imaginary part of "cpow_towardzero":
@@ -850,9 +850,9 @@ idouble: 1
 ifloat: 1
 
 Function: Real part of "csin_downward":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 
 Function: Imaginary part of "csin_downward":
@@ -862,9 +862,9 @@ idouble: 1
 ifloat: 1
 
 Function: Real part of "csin_towardzero":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 
 Function: Imaginary part of "csin_towardzero":
@@ -902,9 +902,9 @@ idouble: 2
 ifloat: 2
 
 Function: Imaginary part of "csinh_downward":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 
 Function: Real part of "csinh_towardzero":
@@ -914,9 +914,9 @@ idouble: 2
 ifloat: 2
 
 Function: Imaginary part of "csinh_towardzero":
-double: 2
+double: 3
 float: 3
-idouble: 2
+idouble: 3
 ifloat: 3
 
 Function: Real part of "csinh_upward":
@@ -1132,15 +1132,15 @@ double: 2
 idouble: 2
 
 Function: "exp10_downward":
-double: 2
+double: 3
 float: 1
-idouble: 2
+idouble: 3
 ifloat: 1
 
 Function: "exp10_towardzero":
-double: 2
+double: 3
 float: 1
-idouble: 2
+idouble: 3
 ifloat: 1
 
 Function: "exp10_upward":
diff --git a/sysdeps/i386/fpu/e_exp_data.c b/sysdeps/i386/fpu/e_exp_data.c
new file mode 100644
index 0000000000..1cc8931700
--- /dev/null
+++ b/sysdeps/i386/fpu/e_exp_data.c
@@ -0,0 +1 @@
+/* Not needed.  */
diff --git a/sysdeps/i386/fpu/math_err.c b/sysdeps/i386/fpu/math_err.c
new file mode 100644
index 0000000000..1cc8931700
--- /dev/null
+++ b/sysdeps/i386/fpu/math_err.c
@@ -0,0 +1 @@
+/* Not needed.  */
diff --git a/sysdeps/i386/fpu/t_exp.c b/sysdeps/i386/fpu/t_exp.c
deleted file mode 100644
index fd37963b05..0000000000
--- a/sysdeps/i386/fpu/t_exp.c
+++ /dev/null
@@ -1 +0,0 @@
-/* Empty.  Not needed. */
diff --git a/sysdeps/ia64/fpu/e_exp_data.c b/sysdeps/ia64/fpu/e_exp_data.c
new file mode 100644
index 0000000000..1cc8931700
--- /dev/null
+++ b/sysdeps/ia64/fpu/e_exp_data.c
@@ -0,0 +1 @@
+/* Not needed.  */
diff --git a/sysdeps/ia64/fpu/math_err.c b/sysdeps/ia64/fpu/math_err.c
new file mode 100644
index 0000000000..1cc8931700
--- /dev/null
+++ b/sysdeps/ia64/fpu/math_err.c
@@ -0,0 +1 @@
+/* Not needed.  */
diff --git a/sysdeps/ia64/fpu/t_exp.c b/sysdeps/ia64/fpu/t_exp.c
deleted file mode 100644
index 41254ae60a..0000000000
--- a/sysdeps/ia64/fpu/t_exp.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 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.                             */