Vector exp for x86_64 and tests.

Here is implementation of vectorized exp containing SSE, AVX,
AVX2 and AVX512 versions according to Vector ABI
<https://groups.google.com/forum/#!topic/x86-64-abi/LmppCfN1rZ4>.

    * bits/libm-simd-decl-stubs.h: Added stubs for exp.
    * math/bits/mathcalls.h: Added exp declaration with __MATHCALL_VEC.
    * sysdeps/unix/sysv/linux/x86_64/libmvec.abilist: New versions added.
    * sysdeps/x86/fpu/bits/math-vector.h: Added SIMD declaration and asm
    redirections for exp.
    * sysdeps/x86_64/fpu/Makefile (libmvec-support): Added new files.
    * sysdeps/x86_64/fpu/Versions: New versions added.
    * sysdeps/x86_64/fpu/libm-test-ulps: Regenerated.
    * sysdeps/x86_64/fpu/multiarch/Makefile (libmvec-sysdep_routines): Added
    build of SSE, AVX2 and AVX512 IFUNC versions.
    * sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core.S: New file.
    * sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core_sse4.S: New file.
    * sysdeps/x86_64/fpu/multiarch/svml_d_exp4_core.S: New file.
    * sysdeps/x86_64/fpu/multiarch/svml_d_exp4_core_avx2.S: New file.
    * sysdeps/x86_64/fpu/multiarch/svml_d_exp8_core.S: New file.
    * sysdeps/x86_64/fpu/multiarch/svml_d_exp8_core_avx512.S: New file.
    * sysdeps/x86_64/fpu/svml_d_exp2_core.S: New file.
    * sysdeps/x86_64/fpu/svml_d_exp4_core.S: New file.
    * sysdeps/x86_64/fpu/svml_d_exp4_core_avx.S: New file.
    * sysdeps/x86_64/fpu/svml_d_exp8_core.S: New file.
    * sysdeps/x86_64/fpu/svml_d_exp_data.S: New file.
    * sysdeps/x86_64/fpu/svml_d_exp_data.h: New file.
    * sysdeps/x86_64/fpu/test-double-vlen2-wrappers.c: Added vector exp test.
    * sysdeps/x86_64/fpu/test-double-vlen2.c: Likewise.
    * sysdeps/x86_64/fpu/test-double-vlen4-avx2-wrappers.c: Likewise.
    * sysdeps/x86_64/fpu/test-double-vlen4-avx2.c: Likewise.
    * sysdeps/x86_64/fpu/test-double-vlen4-wrappers.c: Likewise.
    * sysdeps/x86_64/fpu/test-double-vlen4.c: Likewise.
    * sysdeps/x86_64/fpu/test-double-vlen8-wrappers.c: Likewise.
    * sysdeps/x86_64/fpu/test-double-vlen8.c: Likewise.
    * NEWS: Mention addition of x86_64 vector exp.

30 files changed, 2329 insertions, 4 deletions

diff --git a/ChangeLog b/ChangeLog
index a0c2a9a0fb..62ddb2f43a 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,37 @@
 2015-06-17  Andrew Senkevich  <andrew.senkevich@intel.com>
 
+	* bits/libm-simd-decl-stubs.h: Added stubs for exp.
+	* math/bits/mathcalls.h: Added exp declaration with __MATHCALL_VEC.
+	* sysdeps/unix/sysv/linux/x86_64/libmvec.abilist: New versions added.
+	* sysdeps/x86/fpu/bits/math-vector.h: Added SIMD declaration and asm
+	redirections for exp.
+	* sysdeps/x86_64/fpu/Makefile (libmvec-support): Added new files.
+	* sysdeps/x86_64/fpu/Versions: New versions added.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Regenerated.
+	* sysdeps/x86_64/fpu/multiarch/Makefile (libmvec-sysdep_routines): Added
+	build of SSE, AVX2 and AVX512 IFUNC versions.
+	* sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core.S: New file.
+	* sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core_sse4.S: New file.
+	* sysdeps/x86_64/fpu/multiarch/svml_d_exp4_core.S: New file.
+	* sysdeps/x86_64/fpu/multiarch/svml_d_exp4_core_avx2.S: New file.
+	* sysdeps/x86_64/fpu/multiarch/svml_d_exp8_core.S: New file.
+	* sysdeps/x86_64/fpu/multiarch/svml_d_exp8_core_avx512.S: New file.
+	* sysdeps/x86_64/fpu/svml_d_exp2_core.S: New file.
+	* sysdeps/x86_64/fpu/svml_d_exp4_core.S: New file.
+	* sysdeps/x86_64/fpu/svml_d_exp4_core_avx.S: New file.
+	* sysdeps/x86_64/fpu/svml_d_exp8_core.S: New file.
+	* sysdeps/x86_64/fpu/svml_d_exp_data.S: New file.
+	* sysdeps/x86_64/fpu/svml_d_exp_data.h: New file.
+	* sysdeps/x86_64/fpu/test-double-vlen2-wrappers.c: Added vector exp test.
+	* sysdeps/x86_64/fpu/test-double-vlen2.c: Likewise.
+	* sysdeps/x86_64/fpu/test-double-vlen4-avx2-wrappers.c: Likewise.
+	* sysdeps/x86_64/fpu/test-double-vlen4-avx2.c: Likewise.
+	* sysdeps/x86_64/fpu/test-double-vlen4-wrappers.c: Likewise.
+	* sysdeps/x86_64/fpu/test-double-vlen4.c: Likewise.
+	* sysdeps/x86_64/fpu/test-double-vlen8-wrappers.c: Likewise.
+	* sysdeps/x86_64/fpu/test-double-vlen8.c: Likewise.
+	* NEWS: Mention addition of x86_64 vector exp.
+
 	* sysdeps/unix/sysv/linux/x86_64/libmvec.abilist: New symbols added.
 	* sysdeps/x86/fpu/bits/math-vector.h: Added SIMD declaration and asm
 	redirections for logf.
diff --git a/NEWS b/NEWS
index c361309b4c..bfb4487b46 100644
--- a/NEWS
+++ b/NEWS
@@ -53,7 +53,7 @@ Version 2.22
   condition in some applications.
 
 * Added vector math library named libmvec with the following vectorized x86_64
-  implementations: cos, cosf, sin, sinf, log, logf.
+  implementations: cos, cosf, sin, sinf, log, logf, exp.
   The library can be disabled with --disable-mathvec. Use of the functions is
   enabled with -fopenmp -ffast-math starting from -O1 for GCC version >= 4.9.0.
   The library is linked in as needed when using -lm (no need to specify -lmvec
diff --git a/bits/libm-simd-decl-stubs.h b/bits/libm-simd-decl-stubs.h
index 6367b775a3..1a8bf6f262 100644
--- a/bits/libm-simd-decl-stubs.h
+++ b/bits/libm-simd-decl-stubs.h
@@ -45,4 +45,8 @@
 #define __DECL_SIMD_logf
 #define __DECL_SIMD_logl
 
+#define __DECL_SIMD_exp
+#define __DECL_SIMD_expf
+#define __DECL_SIMD_expl
+
 #endif
diff --git a/math/bits/mathcalls.h b/math/bits/mathcalls.h
index c41946fa0e..1d0c6bd347 100644
--- a/math/bits/mathcalls.h
+++ b/math/bits/mathcalls.h
@@ -97,7 +97,7 @@ __END_NAMESPACE_C99
 
 _Mdouble_BEGIN_NAMESPACE
 /* Exponential function of X.  */
-__MATHCALL (exp,, (_Mdouble_ __x));
+__MATHCALL_VEC (exp,, (_Mdouble_ __x));
 
 /* Break VALUE into a normalized fraction and an integral power of 2.  */
 __MATHCALL (frexp,, (_Mdouble_ __x, int *__exponent));
diff --git a/sysdeps/unix/sysv/linux/x86_64/libmvec.abilist b/sysdeps/unix/sysv/linux/x86_64/libmvec.abilist
index 3593edcbfa..ff9431fa10 100644
--- a/sysdeps/unix/sysv/linux/x86_64/libmvec.abilist
+++ b/sysdeps/unix/sysv/linux/x86_64/libmvec.abilist
@@ -1,18 +1,21 @@
 GLIBC_2.22
  GLIBC_2.22 A
  _ZGVbN2v_cos F
+ _ZGVbN2v_exp F
  _ZGVbN2v_log F
  _ZGVbN2v_sin F
  _ZGVbN4v_cosf F
  _ZGVbN4v_logf F
  _ZGVbN4v_sinf F
  _ZGVcN4v_cos F
+ _ZGVcN4v_exp F
  _ZGVcN4v_log F
  _ZGVcN4v_sin F
  _ZGVcN8v_cosf F
  _ZGVcN8v_logf F
  _ZGVcN8v_sinf F
  _ZGVdN4v_cos F
+ _ZGVdN4v_exp F
  _ZGVdN4v_log F
  _ZGVdN4v_sin F
  _ZGVdN8v_cosf F
@@ -22,5 +25,6 @@ GLIBC_2.22
  _ZGVeN16v_logf F
  _ZGVeN16v_sinf F
  _ZGVeN8v_cos F
+ _ZGVeN8v_exp F
  _ZGVeN8v_log F
  _ZGVeN8v_sin F
diff --git a/sysdeps/x86/fpu/bits/math-vector.h b/sysdeps/x86/fpu/bits/math-vector.h
index 5c3e492ef9..9a353bc62a 100644
--- a/sysdeps/x86/fpu/bits/math-vector.h
+++ b/sysdeps/x86/fpu/bits/math-vector.h
@@ -40,6 +40,8 @@
 #  define __DECL_SIMD_log __DECL_SIMD_x86_64
 #  undef __DECL_SIMD_logf
 #  define __DECL_SIMD_logf __DECL_SIMD_x86_64
+#  undef __DECL_SIMD_exp
+#  define __DECL_SIMD_exp __DECL_SIMD_x86_64
 
 /* Workaround to exclude unnecessary symbol aliases in libmvec
    while GCC creates the vector names based on scalar asm name.
@@ -53,6 +55,10 @@ __asm__ ("_ZGVbN4v___logf_finite = _ZGVbN4v_logf");
 __asm__ ("_ZGVcN8v___logf_finite = _ZGVcN8v_logf");
 __asm__ ("_ZGVdN8v___logf_finite = _ZGVdN8v_logf");
 __asm__ ("_ZGVeN16v___logf_finite = _ZGVeN16v_logf");
+__asm__ ("_ZGVbN2v___exp_finite = _ZGVbN2v_exp");
+__asm__ ("_ZGVcN4v___exp_finite = _ZGVcN4v_exp");
+__asm__ ("_ZGVdN4v___exp_finite = _ZGVdN4v_exp");
+__asm__ ("_ZGVeN8v___exp_finite = _ZGVeN8v_exp");
 
 # endif
 #endif
diff --git a/sysdeps/x86_64/fpu/Makefile b/sysdeps/x86_64/fpu/Makefile
index b610e3faf5..bd6d693d39 100644
--- a/sysdeps/x86_64/fpu/Makefile
+++ b/sysdeps/x86_64/fpu/Makefile
@@ -10,7 +10,8 @@ libmvec-support += svml_d_cos2_core svml_d_cos4_core_avx \
 		   svml_d_log2_core svml_d_log4_core_avx svml_d_log4_core \
 		   svml_d_log8_core svml_d_log_data svml_s_logf4_core \
 		   svml_s_logf8_core_avx svml_s_logf8_core svml_s_logf16_core \
-		   svml_s_logf_data \
+		   svml_s_logf_data svml_d_exp2_core svml_d_exp4_core_avx \
+		   svml_d_exp4_core svml_d_exp8_core svml_d_exp_data \
 		   init-arch
 endif
 
diff --git a/sysdeps/x86_64/fpu/Versions b/sysdeps/x86_64/fpu/Versions
index ecd1b7023e..00e34e771b 100644
--- a/sysdeps/x86_64/fpu/Versions
+++ b/sysdeps/x86_64/fpu/Versions
@@ -3,6 +3,7 @@ libmvec {
     _ZGVbN2v_cos; _ZGVcN4v_cos; _ZGVdN4v_cos; _ZGVeN8v_cos;
     _ZGVbN2v_sin; _ZGVcN4v_sin; _ZGVdN4v_sin; _ZGVeN8v_sin;
     _ZGVbN2v_log; _ZGVcN4v_log; _ZGVdN4v_log; _ZGVeN8v_log;
+    _ZGVbN2v_exp; _ZGVcN4v_exp; _ZGVdN4v_exp; _ZGVeN8v_exp;
     _ZGVbN4v_cosf; _ZGVcN8v_cosf; _ZGVdN8v_cosf; _ZGVeN16v_cosf;
     _ZGVbN4v_sinf; _ZGVcN8v_sinf; _ZGVdN8v_sinf; _ZGVeN16v_sinf;
     _ZGVbN4v_logf; _ZGVcN8v_logf; _ZGVdN8v_logf; _ZGVeN16v_logf;
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index 1812370ecc..45ebc04e58 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1535,6 +1535,18 @@ idouble: 1
 ildouble: 1
 ldouble: 1
 
+Function: "exp_vlen2":
+double: 1
+
+Function: "exp_vlen4":
+double: 1
+
+Function: "exp_vlen4_avx2":
+double: 1
+
+Function: "exp_vlen8":
+double: 1
+
 Function: "expm1":
 double: 1
 float: 1
diff --git a/sysdeps/x86_64/fpu/multiarch/Makefile b/sysdeps/x86_64/fpu/multiarch/Makefile
index 5fc6ea3d23..d6355ae98d 100644
--- a/sysdeps/x86_64/fpu/multiarch/Makefile
+++ b/sysdeps/x86_64/fpu/multiarch/Makefile
@@ -62,5 +62,6 @@ libmvec-sysdep_routines += svml_d_cos2_core_sse4 svml_d_cos4_core_avx2 \
 			   svml_s_cosf16_core_avx512 svml_s_sinf4_core_sse4 \
 			   svml_s_sinf8_core_avx2 svml_s_sinf16_core_avx512 \
 			   svml_s_logf4_core_sse4 svml_s_logf8_core_avx2 \
-			   svml_s_logf16_core_avx512
+			   svml_s_logf16_core_avx512 svml_d_exp2_core_sse4 \
+			   svml_d_exp4_core_avx2 svml_d_exp8_core_avx512
 endif
diff --git a/sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core.S b/sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core.S
new file mode 100644
index 0000000000..ef3dc49a1c
--- /dev/null
+++ b/sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core.S
@@ -0,0 +1,38 @@
+/* Multiple versions of vectorized exp.
+   Copyright (C) 2014-2015 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>
+#include <init-arch.h>
+
+	.text
+ENTRY (_ZGVbN2v_exp)
+        .type   _ZGVbN2v_exp, @gnu_indirect_function
+        cmpl    $0, KIND_OFFSET+__cpu_features(%rip)
+        jne     1f
+        call    __init_cpu_features
+1:      leaq    _ZGVbN2v_exp_sse4(%rip), %rax
+        testl   $bit_SSE4_1, __cpu_features+CPUID_OFFSET+index_SSE4_1(%rip)
+        jz      2f
+        ret
+2:      leaq    _ZGVbN2v_exp_sse2(%rip), %rax
+        ret
+END (_ZGVbN2v_exp)
+libmvec_hidden_def (_ZGVbN2v_exp)
+
+#define _ZGVbN2v_exp _ZGVbN2v_exp_sse2
+#include "../svml_d_exp2_core.S"
diff --git a/sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core_sse4.S b/sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core_sse4.S
new file mode 100644
index 0000000000..1f5445924a
--- /dev/null
+++ b/sysdeps/x86_64/fpu/multiarch/svml_d_exp2_core_sse4.S
@@ -0,0 +1,225 @@
+/* Function exp vectorized with SSE4.
+   Copyright (C) 2014-2015 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>
+#include "svml_d_exp_data.h"
+
+	.text
+ENTRY (_ZGVbN2v_exp_sse4)
+/*
+   ALGORITHM DESCRIPTION:
+
+     Argument representation:
+     N = rint(X*2^k/ln2) = 2^k*M+j
+     X = N*ln2/2^k + r = M*ln2 + ln2*(j/2^k) + r
+     then -ln2/2^(k+1) < r < ln2/2^(k+1)
+     Alternatively:
+     N = trunc(X*2^k/ln2)
+     then 0 < r < ln2/2^k
+
+     Result calculation:
+     exp(X) = exp(M*ln2 + ln2*(j/2^k) + r)
+     = 2^M * 2^(j/2^k) * exp(r)
+     2^M is calculated by bit manipulation
+     2^(j/2^k) is stored in table
+     exp(r) is approximated by polynomial.
+
+     The table lookup is skipped if k = 0.  */
+
+        pushq     %rbp
+        cfi_adjust_cfa_offset (8)
+        cfi_rel_offset (%rbp, 0)
+        movq      %rsp, %rbp
+        cfi_def_cfa_register (%rbp)
+        andq      $-64, %rsp
+        subq      $320, %rsp
+        movaps    %xmm0, %xmm3
+        movq      __svml_dexp_data@GOTPCREL(%rip), %r8
+
+/* iAbsX = (int)(lX>>32), lX = *(longlong*)&X */
+        pshufd    $221, %xmm3, %xmm7
+        movups __dbInvLn2(%r8), %xmm0
+
+/* dK = X*dbInvLn2 */
+        mulpd     %xmm3, %xmm0
+        movq __iAbsMask(%r8), %xmm5
+        movq __iDomainRange(%r8), %xmm6
+
+/* iAbsX = iAbsX&iAbsMask */
+        pand      %xmm5, %xmm7
+
+/* iRangeMask = (iAbsX>iDomainRange) */
+        pcmpgtd   %xmm6, %xmm7
+
+/* Mask = iRangeMask?1:0, set mask for overflow/underflow */
+        movmskps  %xmm7, %eax
+
+/* dN = rint(X*2^k/Ln2) */
+        xorps     %xmm7, %xmm7
+        movups __dbLn2hi(%r8), %xmm5
+        movups __dbLn2lo(%r8), %xmm6
+        roundpd   $0, %xmm0, %xmm7
+
+/* dR = X - dN*dbLn2hi, dbLn2hi is 52-8-k hi bits of ln2/2^k */
+        mulpd     %xmm7, %xmm5
+
+/* dR = dR - dN*dbLn2lo, dbLn2lo is 40..94 bits of lo part of ln2/2^k */
+        mulpd     %xmm6, %xmm7
+        movups __dbShifter(%r8), %xmm4
+
+/* dM = X*dbInvLn2+dbShifter */
+        addpd     %xmm0, %xmm4
+        movaps    %xmm3, %xmm0
+        subpd     %xmm5, %xmm0
+        subpd     %xmm7, %xmm0
+        movups __dPC2(%r8), %xmm5
+
+/* exp(r) = b0+r*(b0+r*(b1+r*b2)) */
+        mulpd     %xmm0, %xmm5
+        addpd __dPC1(%r8), %xmm5
+        mulpd     %xmm0, %xmm5
+        movups __dPC0(%r8), %xmm6
+        addpd     %xmm6, %xmm5
+        mulpd     %xmm5, %xmm0
+        movdqu __lIndexMask(%r8), %xmm2
+
+/* lIndex = (*(longlong*)&dM)&lIndexMask, lIndex is the lower K bits of lM */
+        movdqa    %xmm2, %xmm1
+
+/* lM = (*(longlong*)&dM)&(~lIndexMask) */
+        pandn     %xmm4, %xmm2
+        pand      %xmm4, %xmm1
+
+/* lM = lM<<(52-K), 2^M */
+        psllq     $42, %xmm2
+
+/* table lookup for dT[j] = 2^(j/2^k) */
+        movd      %xmm1, %edx
+        pextrw    $4, %xmm1, %ecx
+        addpd     %xmm0, %xmm6
+        shll      $3, %edx
+        shll      $3, %ecx
+        movq      (%r8,%rdx), %xmm0
+        andl      $3, %eax
+        movhpd    (%r8,%rcx), %xmm0
+
+/* 2^(j/2^k) * exp(r) */
+        mulpd     %xmm6, %xmm0
+
+/* multiply by 2^M through integer add */
+        paddq     %xmm2, %xmm0
+        jne       .LBL_1_3
+
+.LBL_1_2:
+        cfi_remember_state