diff --git a/math/Makefile b/math/Makefile
index 0601f3ac43..04586156f8 100644
--- a/math/Makefile
+++ b/math/Makefile
@@ -115,7 +115,7 @@ type-double-routines := branred doasin dosincos halfulp mpa mpatan2	\
 
 # float support
 type-float-suffix := f
-type-float-routines := k_rem_pio2f
+type-float-routines := k_rem_pio2f math_errf e_exp2f_data
 
 # _Float128 support
 type-float128-suffix := f128
diff --git a/math/w_expf_compat.c b/math/w_expf_compat.c
new file mode 100644
index 0000000000..8a1fa51e46
--- /dev/null
+++ b/math/w_expf_compat.c
@@ -0,0 +1,35 @@
+/* Copyright (C) 2011-2017 Free Software Foundation, Inc.
+   This file is part of the GNU C Library.
+   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
+
+   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 <math_private.h>
+#include <math-svid-compat.h>
+
+/* wrapper expf */
+float
+__expf (float x)
+{
+  float z = __ieee754_expf (x);
+  if (__builtin_expect (!isfinite (z) || z == 0, 0)
+      && isfinite (x) && _LIB_VERSION != _IEEE_)
+    return __kernel_standard_f (x, x, 106 + !!signbit (x));
+
+  return z;
+}
+hidden_def (__expf)
+weak_alias (__expf, expf)
diff --git a/sysdeps/aarch64/fpu/math_private.h b/sysdeps/aarch64/fpu/math_private.h
index 807111ea5a..775237ffb2 100644
--- a/sysdeps/aarch64/fpu/math_private.h
+++ b/sysdeps/aarch64/fpu/math_private.h
@@ -319,6 +319,26 @@ libc_feresetround_noex_aarch64_ctx (struct rm_ctx *ctx)
 #define libc_feresetround_noexf_ctx	libc_feresetround_noex_aarch64_ctx
 #define libc_feresetround_noexl_ctx	libc_feresetround_noex_aarch64_ctx
 
+/* Hack: only include the large arm_neon.h when needed.  */
+#ifdef _MATH_CONFIG_H
+#include <arm_neon.h>
+
+/* ACLE intrinsics for frintn and fcvtns instructions.  */
+#define TOINT_INTRINSICS 1
+
+static inline double_t
+roundtoint (double_t x)
+{
+  return vget_lane_f64 (vrndn_f64 (vld1_f64 (&x)), 0);
+}
+
+static inline uint64_t
+converttoint (double_t x)
+{
+  return vcvtnd_s64_f64 (x);
+}
+#endif
+
 #include_next <math_private.h>
 
 #endif
diff --git a/sysdeps/ieee754/flt-32/e_exp2f.c b/sysdeps/ieee754/flt-32/e_exp2f.c
index 567d3ff6d0..dfadd85f93 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,75 @@
    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)
+__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 (__builtin_expect (abstop >= top12 (128.0f), 0))
     {
-      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 = (double) (xd + SHIFT); /* Rounding to double precision is required.  */
+  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)
+weak_alias (__exp2f, exp2f)
+strong_alias (__exp2f, __ieee754_exp2f)
+strong_alias (__exp2f, __exp2f_finite)
diff --git a/sysdeps/ieee754/flt-32/e_exp2f_data.c b/sysdeps/ieee754/flt-32/e_exp2f_data.c
new file mode 100644
index 0000000000..390dcae333
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/e_exp2f_data.c
@@ -0,0 +1,44 @@
+/* Shared data between expf, exp2f and powf.
+   Copyright (C) 2017 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_config.h"
+
+#define N (1 << EXP2F_TABLE_BITS)
+
+const struct exp2f_data __exp2f_data = {
+  /* tab[i] = uint(2^(i/N)) - (i << 52-BITS)
+     used for computing 2^(k/N) for an int |k| < 150 N as
+     double(tab[k%N] + (k << 52-BITS)) */
+  .tab = {
+0x3ff0000000000000, 0x3fefd9b0d3158574, 0x3fefb5586cf9890f, 0x3fef9301d0125b51,
+0x3fef72b83c7d517b, 0x3fef54873168b9aa, 0x3fef387a6e756238, 0x3fef1e9df51fdee1,
+0x3fef06fe0a31b715, 0x3feef1a7373aa9cb, 0x3feedea64c123422, 0x3feece086061892d,
+0x3feebfdad5362a27, 0x3feeb42b569d4f82, 0x3feeab07dd485429, 0x3feea47eb03a5585,
+0x3feea09e667f3bcd, 0x3fee9f75e8ec5f74, 0x3feea11473eb0187, 0x3feea589994cce13,
+0x3feeace5422aa0db, 0x3feeb737b0cdc5e5, 0x3feec49182a3f090, 0x3feed503b23e255d,
+0x3feee89f995ad3ad, 0x3feeff76f2fb5e47, 0x3fef199bdd85529c, 0x3fef3720dcef9069,
+0x3fef5818dcfba487, 0x3fef7c97337b9b5f, 0x3fefa4afa2a490da, 0x3fefd0765b6e4540,
+  },
+  .shift_scaled = 0x1.8p+52 / N,
+  .poly = { 0x1.c6af84b912394p-5, 0x1.ebfce50fac4f3p-3, 0x1.62e42ff0c52d6p-1 },
+  .shift = 0x1.8p+52,
+  .invln2_scaled = 0x1.71547652b82fep+0 * N,
+  .poly_scaled = {
+0x1.c6af84b912394p-5/N/N/N, 0x1.ebfce50fac4f3p-3/N/N, 0x1.62e42ff0c52d6p-1/N
+  },
+};
diff --git a/sysdeps/ieee754/flt-32/e_expf.c b/sysdeps/ieee754/flt-32/e_expf.c
index 782072f213..d1eabdae08 100644
--- a/sysdeps/ieee754/flt-32/e_expf.c
+++ b/sysdeps/ieee754/flt-32/e_expf.c
@@ -1,7 +1,6 @@
-/* Single-precision floating point e^x.
-   Copyright (C) 1997-2017 Free Software Foundation, Inc.
+/* Single-precision e^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,117 +16,90 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
-/* How this works:
-
-   The input value, x, is written as
-
-   x = n * ln(2) + t/512 + delta[t] + x;
-
-   where:
-   - n is an integer, 127 >= n >= -150;
-   - t is an integer, 177 >= t >= -177
-   - delta is based on a table entry, delta[t] < 2^-28
-   - x is whatever is left, |x| < 2^-10
-
-   Then e^x is approximated as
-
-   e^x = 2^n ( e^(t/512 + delta[t])
-	       + ( e^(t/512 + delta[t])
-		   * ( p(x + delta[t] + n * ln(2)) - delta ) ) )
-
-   where
-   - p(x) is a polynomial approximating e(x)-1;
-   - e^(t/512 + delta[t]) is obtained from a table.
-
-   The table used is the same one as for the double precision version;
-   since we have the table, we might as well use it.
-
-   It turns out to be faster to do calculations in double precision than
-   to perform an 'accurate table method' expf, because of the range reduction
-   overhead (compare exp2f).
-   */
-#include <float.h>
-#include <ieee754.h>
 #include <math.h>
-#include <fenv.h>
-#include <inttypes.h>
-#include <math_private.h>
-
-extern const float __exp_deltatable[178];
-extern const double __exp_atable[355] /* __attribute__((mode(DF))) */;
-
-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 [-ln2/64, ln2/64] (before rounding.)
+Wrong count: 170635 (all nearest rounding wrong results with fma.)
+Non-nearest ULP error: 1 (rounded ULP error)
+*/
+
+#define N (1 << EXP2F_TABLE_BITS)
+#define InvLn2N __exp2f_data.invln2_scaled
+#define T __exp2f_data.tab
+#define C __exp2f_data.poly_scaled
+
+static inline uint32_t
+top12 (float x)
+{
+  return asuint (x) >> 20;
+}
 
 float
-__ieee754_expf (float x)
+__expf (float x)
 {
-  static const float himark = 88.72283935546875;
-  static const float lomark = -103.972084045410;
-  /* Check for usual case.  */
-  if (isless (x, himark) && isgreater (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 (__builtin_expect (abstop >= top12 (88.0f), 0))
     {
-      static const float THREEp42 = 13194139533312.0;
-      static const float THREEp22 = 12582912.0;
-      /* 1/ln(2).  */
-#undef M_1_LN2
-      static const float M_1_LN2 = 1.44269502163f;
-      /* ln(2) */
-#undef M_LN2
-      static const double M_LN2 = .6931471805599452862;
-
-      int tval;
-      double x22, t, result, dx;
-      float n, delta;
-      union ieee754_double ex2_u;
-
-      {
-	SET_RESTORE_ROUND_NOEXF (FE_TONEAREST);
-
-	/* Calculate n.  */
-	n = x * M_1_LN2 + THREEp22;
-	n -= THREEp22;
-	dx = x - n*M_LN2;
-
-	/* Calculate t/512.  */
-	t = dx + THREEp42;
-	t -= THREEp42;
-	dx -= t;
-
-	/* Compute tval = t.  */
-	tval = (int) (t * 512.0);
-
-	if (t >= 0)
-	  delta = - __exp_deltatable[tval];
-	else
-	  delta = __exp_deltatable[-tval];
-
-	/* Compute ex2 = 2^n e^(t/512+delta[t]).  */
-	ex2_u.d = __exp_atable[tval+177];
-	ex2_u.ieee.exponent += (int) n;
-
-	/* Approximate e^(dx+delta) - 1, using a second-degree polynomial,
-	   with maximum error in [-2^-10-2^-28,2^-10+2^-28]
-	   less than 5e-11.  */
-	x22 = (0.5000000496709180453 * dx + 1.0000001192102037084) * dx + delta;
-      }
-
-      /* Return result.  */
-      result = x22 * ex2_u.d + ex2_u.d;
-      return (float) result;
+      /* |x| >= 88 or x is nan.  */
+      if (asuint (x) == asuint (-INFINITY))
+	return 0.0f;
+      if (abstop >= top12 (INFINITY))
+	return x + x;
+      if (x > 0x1.62e42ep6f) /* x > log(0x1p128) ~= 88.72 */
+	return __math_oflowf (0);
+      if (x < -0x1.9fe368p6f) /* x < log(0x1p-150) ~= -103.97 */
+	return __math_uflowf (0);
+#if WANT_ERRNO_UFLOW
+      if (x < -0x1.9d1d9ep6f) /* x < log(0x1p-149) ~= -103.28 */
+	return __math_may_uflowf (0);
+#endif
     }
-  /* Exceptional cases:  */
-  else if (isless (x, himark))
-    {
-      if (isinf (x))
-	/* e^-inf == 0, with no error.  */
-	return 0;
-      else
-	/* Underflow */
-	return TWOM100 * TWOM100;
-    }
-  else
-    /* Return x, if x is a NaN or Inf; or overflow, otherwise.  */
-    return TWO127*x;
+
+  /* x*N/Ln2 = k + r with r in [-1/2, 1/2] and int k.  */
+  z = InvLn2N * xd;
+
+  /* Round and convert z to int, the result is in [-150*N, 128*N] and
+     ideally ties-to-even rule is used, otherwise the magnitude of r
+     can be bigger which gives larger approximation error.  */
+#if TOINT_INTRINSICS
+  kd = roundtoint (z);
+  ki = converttoint (z);
+#elif TOINT_RINT
+  kd = rint (z);
+  ki = (long) kd;
+#elif TOINT_SHIFT
+# define SHIFT __exp2f_data.shift
+  kd = (double) (z + SHIFT); /* Rounding to double precision is required.  */
+  ki = asuint64 (kd);
+  kd -= SHIFT;
+#endif
+  r = z - kd;
+
+  /* exp(x) = 2^(k/N) * 2^(r/N) ~= 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_expf, __expf_finite)
+hidden_def (__expf)
+weak_alias (__expf, expf)
+strong_alias (__expf, __ieee754_expf)
+strong_alias (__expf, __expf_finite)
diff --git a/sysdeps/ieee754/flt-32/math_config.h b/sysdeps/ieee754/flt-32/math_config.h
new file mode 100644
index 0000000000..aaaa90067d
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/math_config.h
@@ -0,0 +1,128 @@
+/* Configuration for math routines.
+   Copyright (C) 2017 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/>.  */
+
+#ifndef _MATH_CONFIG_H
+#define _MATH_CONFIG_H
+
+#include <math.h>
+#include <math_private.h>
+#include <stdint.h>
+
+#ifndef WANT_ROUNDING
+/* Correct special case results in non-nearest rounding modes.  */
+#define WANT_ROUNDING 1
+#endif
+#ifndef WANT_ERRNO
+/* Set errno according to ISO C with (math_errhandling & MATH_ERRNO) != 0.  */
+#define WANT_ERRNO 1
+#endif
+#ifndef WANT_ERRNO_UFLOW
+/* Set errno to ERANGE if result underflows to 0 (in all rounding modes).  */
+#define WANT_ERRNO_UFLOW (WANT_ROUNDING && WANT_ERRNO)
+#endif
+
+#ifndef TOINT_INTRINSICS
+#define TOINT_INTRINSICS 0
+#endif
+#ifndef TOINT_RINT
+#define TOINT_RINT 0
+#endif
+#ifndef TOINT_SHIFT
+#define TOINT_SHIFT 1
+#endif
+
+static inline uint32_t
+asuint (float f)
+{
+  union
+  {
+    float f;
+    uint32_t i;
+  } u = {f};
+  return u.i;
+}
+
+static inline float
+asfloat (uint32_t i)
+{
+  union
+  {
+    uint32_t i;
+    float f;
+  } u = {i};
+  return u.f;
+}
+
+static inline uint64_t
+asuint64 (double f)
+{
+  union
+  {
+    double f;
+    uint64_t i;
+  } u = {f};
+  return u.i;
+}
+
+static inline double
+asdouble (uint64_t i)
+{
+  union
+  {
+    uint64_t i;
+    double f;
+  } u = {i};
+  return u.f;
+}
+
+static inline int
+ieee_2008_issignaling (float x)
+{
+  uint32_t ix = asuint (x);
+  ix ^= 0x00400000; /* IEEE 754-2008 snan bit.  */
+  return 2 * ix > 2u * 0x7fc00000;
+}
+
+#ifdef __GNUC__
+#define HIDDEN __attribute__ ((__visibility__ ("hidden")))
+#define NOINLINE __attribute__ ((noinline))
+#else
+#define HIDDEN
+#define NOINLINE
+#endif
+
+HIDDEN float __math_oflowf (unsigned long);
+HIDDEN float __math_uflowf (unsigned long);
+HIDDEN float __math_may_uflowf (unsigned long);
+HIDDEN float __math_divzerof (unsigned long);
+HIDDEN float __math_invalidf (float);
+
+/* Shared between expf, exp2f and powf.  */
+#define EXP2F_TABLE_BITS 5
+#define EXP2F_POLY_ORDER 3
+extern const struct exp2f_data
+{
+  uint64_t tab[1 << EXP2F_TABLE_BITS];
+  double shift_scaled;
+  double poly[EXP2F_POLY_ORDER];
+  double shift;
+  double invln2_scaled;
+  double poly_scaled[EXP2F_POLY_ORDER];
+} __exp2f_data;
+
+#endif
diff --git a/sysdeps/ieee754/flt-32/math_errf.c b/sysdeps/ieee754/flt-32/math_errf.c
new file mode 100644
index 0000000000..11fbdc0b92
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/math_errf.c
@@ -0,0 +1,77 @@
+/* Single-precision math error handling.
+   Copyright (C) 2017 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_config.h"
+
+#if WANT_ERRNO
+#include <errno.h>
+/* NOINLINE reduces code size and avoids making math functions non-leaf
+   when the error handling is inlined.  */
+NOINLINE static float
+with_errnof (float y, int e)
+{
+  errno = e;
+  return y;
+}
+#else
+#define with_errnof(x, e) (x)
+#endif
+
+/* NOINLINE prevents fenv semantics breaking optimizations.  */
+NOINLINE static float
+xflowf (unsigned long sign, float y)
+{
+  y = (sign ? -y : y) * y;
+  return with_errnof (y, ERANGE);
+}
+
+HIDDEN float
+__math_uflowf (unsigned long sign)
+{
+  return xflowf (sign, 0x1p-95f);
+}
+
+#if WANT_ERRNO_UFLOW
+/* Underflows to zero in some non-nearest rounding mode, setting errno
+   is valid even if the result is non-zero, but in the subnormal range.  */
+HIDDEN float
+__math_may_uflowf (unsigned long sign)
+{
+  return xflowf (sign, 0x1.4p-75f);
+}
+#endif
+
+HIDDEN float
+__math_oflowf (unsigned long sign)
+{
+  return xflowf (sign, 0x1p97f);
+}
+
+HIDDEN float
+__math_divzerof (unsigned long sign)
+{
+  float y = 0;
+  return with_errnof ((sign ? -1 : 1) / y, ERANGE);
+}
+
+HIDDEN float
+__math_invalidf (float x)
+{
+  float y = (x - x) / (x - x);
+  return isnan (x) ? y : with_errnof (y, EDOM);
+}
diff --git a/sysdeps/ieee754/flt-32/w_exp2f_compat.c b/sysdeps/ieee754/flt-32/w_exp2f_compat.c
new file mode 100644
index 0000000000..c37de4df43
--- /dev/null
+++ b/sysdeps/ieee754/flt-32/w_exp2f_compat.c
@@ -0,0 +1 @@
+/* Empty file to disable the exp2f wrapper.  */
diff --git a/sysdeps/ieee754/flt-32/w_expf_compat.c b/sysdeps/ieee754/flt-32/w_expf_compat.c
index 8a1fa51e46..578cf4d4a7 100644
--- a/sysdeps/ieee754/flt-32/w_expf_compat.c
+++ b/sysdeps/ieee754/flt-32/w_expf_compat.c
@@ -1,35 +1 @@
-/* Copyright (C) 2011-2017 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@gmail.com>, 2011.
-
-   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 <math_private.h>
-#include <math-svid-compat.h>
-
-/* wrapper expf */
-float
-__expf (float x)
-{
-  float z = __ieee754_expf (x);
-  if (__builtin_expect (!isfinite (z) || z == 0, 0)
-      && isfinite (x) && _LIB_VERSION != _IEEE_)
-    return __kernel_standard_f (x, x, 106 + !!signbit (x));
-
-  return z;
-}
-hidden_def (__expf)
-weak_alias (__expf, expf)
+/* Empty file to disable the expf wrapper.  */
diff --git a/sysdeps/x86_64/fpu/w_expf_compat.c b/sysdeps/x86_64/fpu/w_expf_compat.c
new file mode 100644
index 0000000000..dc454638c3
--- /dev/null
+++ b/sysdeps/x86_64/fpu/w_expf_compat.c
@@ -0,0 +1 @@
+#include <math/w_expf_compat.c>
