Optimized generic expf and exp2f with wrappers

Based on new expf and exp2f code from
https://github.com/ARM-software/optimized-routines/

with wrapper on aarch64:
expf reciprocal-throughput: 2.3x faster
expf latency: 1.7x faster
without wrapper on aarch64:
expf reciprocal-throughput: 3.3x faster
expf latency: 1.7x faster
without wrapper on aarch64:
exp2f reciprocal-throughput: 2.8x faster
exp2f latency: 1.3x faster
libm.so size on aarch64:
.text size: -152 bytes
.rodata size: -1740 bytes
expf/exp2f worst case nearest rounding error: 0.502 ulp
worst case non-nearest rounding error: 1 ulp

Error checks are inline and errno setting is in separate tail called
functions, but the wrappers are kept in this patch to handle the
_LIB_VERSION==_SVID_ case.  (So e.g. errno is set twice for expf calls
and once for __expf_finite calls on targets where the new code is used.)

Double precision arithmetics is used which is expected to be faster on
most targets (including soft-float) than using single precision and it
is easier to get good precision result with it.

Const data is kept in a separate translation unit which complicates
maintenance a bit, but is expected to give good code for literal loads
on most targets and allows sharing data across expf, exp2f and powf.
(This data is disabled on i386, m68k and ia64 which have their own
expf, exp2f and powf code.)

Some details may need target specific tweaks:
- best convert and round to int operation in the arg reduction may be
different across targets.
- code was optimized on fma target, optimal polynomial eval may be
different without fma.
- gcc does not always generate good code for fp bit representation
access via unions or it may be inherently slow on some targets.

The libm-test-ulps will need adjustment because..
- The argument reduction ideally uses nearest rounded rint, but that is
not efficient on most targets, so the polynomial can get evaluated on a
wider interval in non-nearest rounding mode making 1 ulp errors common
in that case.
- The polynomial is evaluated such that it may have 1 ulp error on
negative tiny inputs with upward rounding.

	* math/Makefile (type-float-routines): Add math_errf and e_exp2f_data.
	* sysdeps/aarch64/fpu/math_private.h (TOINT_INTRINSICS): Define.
	(roundtoint, converttoint): Likewise.
	* sysdeps/ieee754/flt-32/e_expf.c: New implementation.
	* sysdeps/ieee754/flt-32/e_exp2f.c: New implementation.
	* sysdeps/ieee754/flt-32/e_exp2f_data.c: New file.
	* sysdeps/ieee754/flt-32/math_config.h: New file.
	* sysdeps/ieee754/flt-32/math_errf.c: New file.
	* sysdeps/ieee754/flt-32/t_exp2f.h: Remove.
	* sysdeps/i386/fpu/e_exp2f_data.c: New file.
	* sysdeps/i386/fpu/math_errf.c: New file.
	* sysdeps/ia64/fpu/e_exp2f_data.c: New file.
	* sysdeps/ia64/fpu/math_errf.c: New file.
	* sysdeps/m68k/m680x0/fpu/e_exp2f_data.c: New file.
	* sysdeps/m68k/m680x0/fpu/math_errf.c: New file.

1 files changed, 63 insertions, 107 deletions

diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c
index 567d3ff6d0..72b7d8829f 100644
--- a/sysdeps/ieee754/flt-32/e_exp2f.c
+++ b/sysdeps/ieee754/flt-32/e_exp2f.c
@@ -1,7 +1,6 @@
-/* Single-precision floating point 2^x.
-   Copyright (C) 1997-2017 Free Software Foundation, Inc.
+/* Single-precision 2^x function.
+   Copyright (C) 2017 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
-   Contributed by Geoffrey Keating <geoffk@ozemail.com.au>
 
    The GNU C Library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
@@ -17,116 +16,73 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
-/* The basic design here is from
-   Shmuel Gal and Boris Bachelis, "An Accurate Elementary Mathematical
-   Library for the IEEE Floating Point Standard", ACM Trans. Math. Soft.,
-   17 (1), March 1991, pp. 26-45.
-   It has been slightly modified to compute 2^x instead of e^x, and for
-   single-precision.
-   */
-#ifndef _GNU_SOURCE
-# define _GNU_SOURCE
-#endif
-#include <stdlib.h>
-#include <float.h>
-#include <ieee754.h>
 #include <math.h>
-#include <fenv.h>
-#include <inttypes.h>
-#include <math_private.h>
-
-#include "t_exp2f.h"
-
-static const float TWOM100 = 7.88860905e-31;
-static const float TWO127 = 1.7014118346e+38;
+#include <stdint.h>
+#include "math_config.h"
+
+/*
+EXP2F_TABLE_BITS = 5
+EXP2F_POLY_ORDER = 3
+
+ULP error: 0.502 (nearest rounding.)
+Relative error: 1.69 * 2^-34 in [-1/64, 1/64] (before rounding.)
+Wrong count: 168353 (all nearest rounding wrong results with fma.)
+Non-nearest ULP error: 1 (rounded ULP error)
+*/
+
+#define N (1 << EXP2F_TABLE_BITS)
+#define T __exp2f_data.tab
+#define C __exp2f_data.poly
+#define SHIFT __exp2f_data.shift_scaled
+
+static inline uint32_t
+top12 (float x)
+{
+  return asuint (x) >> 20;
+}
 
 float
 __ieee754_exp2f (float x)
 {
-  static const float himark = (float) FLT_MAX_EXP;
-  static const float lomark = (float) (FLT_MIN_EXP - FLT_MANT_DIG - 1);
-
-  /* Check for usual case.  */
-  if (isless (x, himark) && isgreaterequal (x, lomark))
+  uint32_t abstop;
+  uint64_t ki, t;
+  /* double_t for better performance on targets with FLT_EVAL_METHOD==2.  */
+  double_t kd, xd, z, r, r2, y, s;
+
+  xd = (double_t) x;
+  abstop = top12 (x) & 0x7ff;
+  if (__glibc_unlikely (abstop >= top12 (128.0f)))
     {
-      static const float THREEp14 = 49152.0;
-      int tval, unsafe;
-      float rx, x22, result;
-      union ieee754_float ex2_u, scale_u;
-
-      if (fabsf (x) < FLT_EPSILON / 4.0f)
-	return 1.0f + x;
-
-      {
-	SET_RESTORE_ROUND_NOEXF (FE_TONEAREST);
-
-	/* 1. Argument reduction.
-	   Choose integers ex, -128 <= t < 128, and some real
-	   -1/512 <= x1 <= 1/512 so that
-	   x = ex + t/512 + x1.
-
-	   First, calculate rx = ex + t/256.  */
-	rx = x + THREEp14;
-	rx -= THREEp14;
-	x -= rx;  /* Compute x=x1. */
-	/* Compute tval = (ex*256 + t)+128.
-	   Now, t = (tval mod 256)-128 and ex=tval/256  [that's mod, NOT %;
-	   and /-round-to-nearest not the usual c integer /].  */
-	tval = (int) (rx * 256.0f + 128.0f);
-
-	/* 2. Adjust for accurate table entry.
-	   Find e so that
-	   x = ex + t/256 + e + x2
-	   where -7e-4 < e < 7e-4, and
-	   (float)(2^(t/256+e))
-	   is accurate to one part in 2^-64.  */
-
-	/* 'tval & 255' is the same as 'tval%256' except that it's always
-	   positive.
-	   Compute x = x2.  */
-	x -= __exp2f_deltatable[tval & 255];
-
-	/* 3. Compute ex2 = 2^(t/255+e+ex).  */
-	ex2_u.f = __exp2f_atable[tval & 255];
-	tval >>= 8;
-	/* x2 is an integer multiple of 2^-30; avoid intermediate
-	   underflow from the calculation of x22 * x.  */
-	unsafe = abs(tval) >= -FLT_MIN_EXP - 32;
-	ex2_u.ieee.exponent += tval >> unsafe;
-	scale_u.f = 1.0;
-	scale_u.ieee.exponent += tval - (tval >> unsafe);
-
-	/* 4. Approximate 2^x2 - 1, using a second-degree polynomial,
-	   with maximum error in [-2^-9 - 2^-14, 2^-9 + 2^-14]
-	   less than 1.3e-10.  */
-
-	x22 = (.24022656679f * x + .69314736128f) * ex2_u.f;
-      }
-
-      /* 5. Return (2^x2-1) * 2^(t/512+e+ex) + 2^(t/512+e+ex).  */
-      result = x22 * x + ex2_u.f;
-
-      if (!unsafe)
-	return result;
-      else
-	{
-	  result *= scale_u.f;
-	  math_check_force_underflow_nonneg (result);
-	  return result;
-	}
-    }
-  /* Exceptional cases:  */
-  else if (isless (x, himark))
-    {
-      if (isinf (x))
-	/* e^-inf == 0, with no error.  */
-	return 0;
-      else
-	/* Underflow */
-	return TWOM100 * TWOM100;
+      /* |x| >= 128 or x is nan.  */
+      if (asuint (x) == asuint (-INFINITY))
+	return 0.0f;
+      if (abstop >= top12 (INFINITY))
+	return x + x;
+      if (x > 0.0f)
+	return __math_oflowf (0);
+      if (x <= -150.0f)
+	return __math_uflowf (0);
+#if WANT_ERRNO_UFLOW
+      if (x < -149.0f)
+	return __math_may_uflowf (0);
+#endif
     }
-  else
-    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
-    return TWO127*x;
+
+  /* x = k/N + r with r in [-1/(2N), 1/(2N)] and int k.  */
+  kd = math_narrow_eval ((double) (xd + SHIFT)); /* Needs to be double.  */
+  ki = asuint64 (kd);
+  kd -= SHIFT; /* k/N for int k.  */
+  r = xd - kd;
+
+  /* exp2(x) = 2^(k/N) * 2^r ~= s * (C0*r^3 + C1*r^2 + C2*r + 1) */
+  t = T[ki % N];
+  t += ki << (52 - EXP2F_TABLE_BITS);
+  s = asdouble (t);
+  z = C[0] * r + C[1];
+  r2 = r * r;
+  y = C[2] * r + 1;
+  y = z * r2 + y;
+  y = y * s;
+  return (float) y;
 }
 strong_alias (__ieee754_exp2f, __exp2f_finite)