x86-64: Replace assembly versions of e_expf with generic e_expf.c

This patch replaces x86-64 assembly versions of e_expf with generic
e_expf.c.  For workload-spec2017.wrf, on Nehalem, it improves
performance by:

                           Before            After     Improvement
reciprocal-throughput      36.039           20.7749       73%
latency                    58.8096          40.8715       43%

On Skylake, it improves

                           Before            After     Improvement
reciprocal-throughput      18.4436          11.1693       65%
latency                    47.5162          37.5411       26%

	* sysdeps/x86_64/fpu/e_expf.S: Removed.
	* sysdeps/x86_64/fpu/multiarch/e_expf-fma.S: Likewise.
	* sysdeps/x86_64/fpu/w_expf.c: Likewise.
	* sysdeps/x86_64/fpu/libm-test-ulps: Updated for generic
	e_expf.c.
	* sysdeps/x86_64/fpu/multiarch/Makefile (CFLAGS-e_expf-fma.c):
	New.
	* sysdeps/x86_64/fpu/multiarch/e_expf-fma.c: New file.
	* sysdeps/x86_64/fpu/multiarch/e_expf.c (__redirect_ieee754_expf):
	Renamed to ...
	(__redirect_expf): This.
	(SYMBOL_NAME): Changed to expf.
	(__ieee754_expf): Renamed to ...
	(__expf): This.
	(__GI___expf): This.
	(__ieee754_expf): Add strong_alias.
	(__expf_finite): Likewise.
	(__expf): New.
	Include <sysdeps/ieee754/flt-32/e_expf.c>.

4 files changed, 25 insertions, 189 deletions

diff --git a/sysdeps/x86_64/fpu/multiarch/Makefile b/sysdeps/x86_64/fpu/multiarch/Makefile
index d66055224a..92e91d79af 100644
--- a/sysdeps/x86_64/fpu/multiarch/Makefile
+++ b/sysdeps/x86_64/fpu/multiarch/Makefile
@@ -37,9 +37,10 @@ CFLAGS-slowpow-fma.c = -mfma -mavx2
 CFLAGS-s_sin-fma.c = -mfma -mavx2
 CFLAGS-s_tan-fma.c = -mfma -mavx2
 
-# e_expf-fma.S implements both FMA and SSE2 versions of e_expf.
 libm-sysdep_routines += e_expf-fma
 
+CFLAGS-e_expf-fma.c = -mfma -mavx2
+
 libm-sysdep_routines += e_exp-fma4 e_log-fma4 e_pow-fma4 s_atan-fma4 \
 			e_asin-fma4 e_atan2-fma4 s_sin-fma4 s_tan-fma4 \
 			mplog-fma4 mpa-fma4 slowexp-fma4 slowpow-fma4 \
diff --git a/sysdeps/x86_64/fpu/multiarch/e_expf-fma.S b/sysdeps/x86_64/fpu/multiarch/e_expf-fma.S
deleted file mode 100644
index 43140deced..0000000000
--- a/sysdeps/x86_64/fpu/multiarch/e_expf-fma.S
+++ /dev/null
@@ -1,182 +0,0 @@
-/* FMA/AVX2 version of IEEE 754 expf.
-   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 <sysdep.h>
-
-/* Short algorithm description:
-
-    Let K = 64 (table size).
-         e^x  = 2^(x/log(2)) = 2^n * T[j] * (1 + P(y))
-    where
-         x = m*log(2)/K + y,    y in [0.0..log(2)/K]
-         m = n*K + j,           m,n,j - signed integer, j in [0..K-1]
-         values of 2^(j/K) are tabulated as T[j].
-
-         P(y) is a minimax polynomial approximation of expf(x)-1
-         on small interval [0.0..log(2)/K].
-
-         P(y) = P3*y*y*y*y + P2*y*y*y + P1*y*y + P0*y, calculated as
-         z = y*y;    P(y) = (P3*z + P1)*z + (P2*z + P0)*y
-
-   Special cases:
-    expf(NaN) = NaN
-    expf(+INF) = +INF
-    expf(-INF) = 0
-    expf(x) = 1 for subnormals
-    for finite argument, only expf(0)=1 is exact
-    expf(x) overflows if x>88.7228317260742190
-    expf(x) underflows if x<-103.972076416015620
- */
-
-	.section .text.fma,"ax",@progbits
-ENTRY(__ieee754_expf_fma)
-	/* Input: single precision x in %xmm0 */
-	vcvtss2sd %xmm0, %xmm0, %xmm1	/* Convert x to double precision */
-	vmovd	%xmm0, %ecx		/* Copy x */
-	vmovsd	L(DP_KLN2)(%rip), %xmm2	/* DP K/log(2) */
-	vfmadd213sd L(DP_RD)(%rip), %xmm1, %xmm2 /* DP x*K/log(2)+RD */
-	vmovsd	L(DP_P2)(%rip), %xmm3	/* DP P2 */
-	movl	%ecx, %eax		/* x */
-	andl	$0x7fffffff, %ecx	/* |x| */
-	lea	L(DP_T)(%rip), %rsi	/* address of table T[j] */
-	vmovsd	L(DP_P3)(%rip), %xmm4	/* DP P3 */
-
-	cmpl	$0x42ad496b, %ecx	/* |x|<125*log(2) ? */
-	jae	L(special_paths_fma)
-
-	/* Here if |x|<125*log(2) */
-	cmpl	$0x31800000, %ecx	/* |x|<2^(-28) ? */
-	jb	L(small_arg_fma)
-
-	/* Main path: here if 2^(-28)<=|x|<125*log(2) */
-						/* %xmm2 = SP x*K/log(2)+RS */
-	vmovd	  %xmm2, %eax
-	vsubsd	  L(DP_RD)(%rip), %xmm2, %xmm2 	/* DP t=round(x*K/log(2)) */
-	movl	  %eax, %edx			/* n*K+j with trash */
-	andl	  $0x3f, %eax			/* bits of j */
-	vmovsd	  (%rsi,%rax,8), %xmm5		/* T[j] */
-	andl	  $0xffffffc0, %edx		/* bits of n */
-
-	vfmadd132sd  L(DP_NLN2K)(%rip), %xmm1, %xmm2 /*  DP y=x-t*log(2)/K */
-	vmulsd	    %xmm2, %xmm2, %xmm6		/* DP z=y*y */
-
-
-	vfmadd213sd L(DP_P1)(%rip), %xmm6, %xmm4 /* DP P3*z + P1 */
-	vfmadd213sd L(DP_P0)(%rip), %xmm6, %xmm3 /* DP P2*z+P0 */
-
-	addl	    $0x1fc0, %edx		/* bits of n + SP exponent bias */
-	shll	    $17, %edx			/* SP 2^n */
-	vmovd       %edx, %xmm1			/* SP 2^n */
-
-	vmulsd      %xmm6, %xmm4, %xmm4		/* DP (P3*z+P1)*z */
-
-	vfmadd213sd %xmm4, %xmm3, %xmm2		/* DP P(Y)  (P2*z+P0)*y */
-	vfmadd213sd %xmm5, %xmm5, %xmm2		/* DP T[j]*(P(y)+1) */
-	vcvtsd2ss   %xmm2, %xmm2, %xmm0		/* SP T[j]*(P(y)+1) */
-	vmulss	    %xmm1, %xmm0, %xmm0		/* SP result=2^n*(T[j]*(P(y)+1)) */
-	ret
-
-	.p2align	4
-L(small_arg_fma):
-	/* Here if 0<=|x|<2^(-28) */
-	vaddss	L(SP_ONE)(%rip), %xmm0, %xmm0	/* 1.0 + x */
-	/* Return 1.0 with inexact raised, except for x==0 */
-	ret
-
-	.p2align	4
-L(special_paths_fma):
-	/* Here if 125*log(2)<=|x| */
-	shrl	$31, %eax		/* Get sign bit of x, and depending on it: */
-	lea	L(SP_RANGE)(%rip), %rdx	/* load over/underflow bound */
-	cmpl	(%rdx,%rax,4), %ecx	/* |x|<under/overflow bound ? */
-	jbe	L(near_under_or_overflow_fma)
-
-	/* Here if |x|>under/overflow bound */
-	cmpl	$0x7f800000, %ecx	/* |x| is finite ? */
-	jae	L(arg_inf_or_nan_fma)
-
-	/* Here if |x|>under/overflow bound, and x is finite */
-	testl	%eax, %eax		/* sign of x nonzero ? */
-	je	L(res_overflow_fma)
-
-	/* Here if -inf<x<underflow bound (x<0) */
-	vmovss	L(SP_SMALL)(%rip), %xmm0/* load small value 2^(-100) */
-	vmulss	%xmm0, %xmm0, %xmm0	/* Return underflowed result (zero or subnormal) */
-	ret
-
-	.p2align	4
-L(res_overflow_fma):
-	/* Here if overflow bound<x<inf (x>0) */
-	vmovss	L(SP_LARGE)(%rip), %xmm0/* load large value 2^100 */
-	vmulss	%xmm0, %xmm0, %xmm0	/* Return overflowed result (Inf or max normal) */
-	ret
-
-	.p2align	4
-L(arg_inf_or_nan_fma):
-	/* Here if |x| is Inf or NAN */
-	jne	L(arg_nan_fma)	/* |x| is Inf ? */
-
-	/* Here if |x| is Inf */
-	lea	L(SP_INF_0)(%rip), %rdx	/* depending on sign of x: */
-	vmovss	(%rdx,%rax,4), %xmm0	/* return zero or Inf */
-	ret
-
-	.p2align	4
-L(arg_nan_fma):
-	/* Here if |x| is NaN */
-	vaddss	%xmm0, %xmm0, %xmm0	/* Return x+x (raise invalid) */
-	ret
-
-	.p2align	4
-L(near_under_or_overflow_fma):
-	/* Here if 125*log(2)<=|x|<under/overflow bound */
-	vmovd	%xmm2, %eax		/* bits of n*K+j with trash */
-	vsubsd	L(DP_RD)(%rip), %xmm2, %xmm2 	/* DP t=round(x*K/log(2)) */
-	movl	%eax, %edx		/* n*K+j with trash */
-	andl	$0x3f, %eax		/* bits of j */
-	vmulsd	L(DP_NLN2K)(%rip),%xmm2, %xmm2/* DP -t*log(2)/K */
-	andl	$0xffffffc0, %edx	/* bits of n */
-	vaddsd	%xmm1, %xmm2, %xmm0	/* DP y=x-t*log(2)/K */
-	vmulsd	%xmm0, %xmm0, %xmm2	/* DP z=y*y */
-	addl	$0xffc0, %edx		/* bits of n + DP exponent bias */
-	vfmadd213sd L(DP_P0)(%rip), %xmm2, %xmm3/* DP P2*z+P0 */
-	shlq	$46, %rdx		/* DP 2^n */
-	vfmadd213sd L(DP_P1)(%rip), %xmm2, %xmm4/* DP P3*z+P1 */
-	vmovq	%rdx, %xmm1		/* DP 2^n */
-	vmulsd	%xmm2, %xmm4, %xmm4	/* DP (P3*z+P1)*z */
-	vfmadd213sd %xmm4, %xmm3, %xmm0	/* DP (P2*z+P0)*y */
-	vmovsd	(%rsi,%rax,8), %xmm2
-	vfmadd213sd %xmm2, %xmm2, %xmm0 /* DP T[j]*(P(y)+1) */
-	vmulsd	%xmm1, %xmm0, %xmm0	/* DP result=2^n*(T[j]*(P(y)+1)) */
-	vcvtsd2ss %xmm0, %xmm0, %xmm0	/* convert result to single precision */
-	ret
-END(__ieee754_expf_fma)
-
-	.section .rodata.cst8,"aM",@progbits,8
-	.p2align 3
-L(DP_RD): /* double precision 2^52+2^51 */
-	.long	0x00000000, 0x43380000
-	.type L(DP_RD), @object
-	ASM_SIZE_DIRECTIVE(L(DP_RD))
-
-#define __ieee754_expf __ieee754_expf_sse2
-
-#undef strong_alias
-#define strong_alias(ignored1, ignored2)
-
-#include <sysdeps/x86_64/fpu/e_expf.S>
diff --git a/sysdeps/x86_64/fpu/multiarch/e_expf-fma.c b/sysdeps/x86_64/fpu/multiarch/e_expf-fma.c
new file mode 100644
index 0000000000..4e01cd6a82
--- /dev/null
+++ b/sysdeps/x86_64/fpu/multiarch/e_expf-fma.c
@@ -0,0 +1,3 @@
+#define __expf __expf_fma
+
+#include <sysdeps/ieee754/flt-32/e_expf.c>
diff --git a/sysdeps/x86_64/fpu/multiarch/e_expf.c b/sysdeps/x86_64/fpu/multiarch/e_expf.c
index 096209857e..18451a0f77 100644
--- a/sysdeps/x86_64/fpu/multiarch/e_expf.c
+++ b/sysdeps/x86_64/fpu/multiarch/e_expf.c
@@ -1,4 +1,4 @@
-/* Multiple versions of IEEE 754 expf.
+/* Multiple versions of expf.
    Copyright (C) 2017 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
@@ -16,11 +16,25 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
-extern float __redirect_ieee754_expf (float);
+extern float __redirect_expf (float);
 
-#define SYMBOL_NAME ieee754_expf
+#define SYMBOL_NAME expf
 #include "ifunc-fma.h"
 
-libc_ifunc_redirected (__redirect_ieee754_expf, __ieee754_expf,
-		       IFUNC_SELECTOR ());
-strong_alias (__ieee754_expf, __expf_finite)
+libc_ifunc_redirected (__redirect_expf, __expf, IFUNC_SELECTOR ());
+
+#ifdef SHARED
+__hidden_ver1 (__expf, __GI___expf, __redirect_expf)
+  __attribute__ ((visibility ("hidden")));
+
+# include <shlib-compat.h>
+versioned_symbol (libm, __expf, expf, GLIBC_2_27);
+#else
+weak_alias (__expf, expf)
+#endif
+
+strong_alias (__expf, __ieee754_expf)
+strong_alias (__expf, __expf_finite)
+
+#define __expf __expf_sse2
+#include <sysdeps/ieee754/flt-32/e_expf.c>