Fix x86/x86_64 expm1l inaccuracy and exceptions (bugs 13885, 13923).

9 files changed, 174 insertions, 270 deletions

diff --git a/ChangeLog b/ChangeLog
index 0cd2161ddc..ba2cedf89c 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,26 @@
+2012-05-07  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #13885]
+	[BZ #13923]
+	* sysdeps/i386/fpu/e_expl.S (IEEE754_EXPL): Define conditional on
+	USE_AS_EXPM1L.
+	(EXPL_FINITE): Likewise.
+	(FLDLOG): Likewise.
+	(IEEE754_EXPL) [USE_AS_EXPM1L]: Support use as expm1l.
+	* sysdeps/i386/fpu/s_expm1l.S: Define USE_AS_EXPM1L and include
+	e_expl.S.
+	* sysdeps/x86_64/fpu/e_expl.S (IEEE754_EXPL): Define conditional on
+	USE_AS_EXPM1L.
+	(EXPL_FINITE): Likewise.
+	(FLDLOG): Likewise.
+	(IEEE754_EXPL) [USE_AS_EXPM1L]: Support use as expm1l.
+	* sysdeps/x86_64/fpu/s_expm1l.S: Define USE_AS_EXPM1L and include
+	e_expl.S.
+	* math/libm-test.inc (expm1_test): Add more tests.  Do not disable
+	test of -max_value argument for long double.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2012-05-06  David S. Miller  <davem@davemloft.net>
 
 	* scripts/data/localplt-sparc-linux-gnu.data: Add '?' markers to
diff --git a/NEWS b/NEWS
index 90662f4fa0..b2f8a4e789 100644
--- a/NEWS
+++ b/NEWS
@@ -20,11 +20,11 @@ Version 2.16
   13583, 13592, 13618, 13637, 13656, 13658, 13673, 13691, 13695, 13704,
   13705, 13706, 13726, 13738, 13739, 13758, 13760, 13761, 13775, 13786,
   13787, 13792, 13806, 13824, 13840, 13841, 13844, 13846, 13851, 13852,
-  13854, 13871, 13872, 13873, 13879, 13883, 13884, 13886, 13892, 13895,
-  13908, 13910, 13911, 13912, 13913, 13914, 13915, 13916, 13917, 13918,
-  13919, 13920, 13921, 13922, 13924, 13926, 13927, 13928, 13938, 13941,
-  13942, 13963, 13967, 13970, 13973, 14027, 14033, 14034, 14040, 14049,
-  14055, 14064
+  13854, 13871, 13872, 13873, 13879, 13883, 13884, 13885, 13886, 13892,
+  13895, 13908, 13910, 13911, 13912, 13913, 13914, 13915, 13916, 13917,
+  13918, 13919, 13920, 13921, 13922, 13923, 13924, 13926, 13927, 13928,
+  13938, 13941, 13942, 13963, 13967, 13970, 13973, 14027, 14033, 14034,
+  14040, 14049, 14055, 14064
 
 * ISO C11 support:
 
diff --git a/math/libm-test.inc b/math/libm-test.inc
index 5fe9e5a14e..542131dc7b 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -3551,6 +3551,13 @@ expm1_test (void)
   TEST_f_f (expm1, 1, M_El - 1.0);
   TEST_f_f (expm1, 0.75L, 1.11700001661267466854536981983709561L);
 
+  TEST_f_f (expm1, 50.0L, 5.1847055285870724640864533229334853848275e+21L);
+
+#ifndef TEST_FLOAT
+  TEST_f_f (expm1, 127.0L, 1.4302079958348104463583671072905261080748e+55L);
+  TEST_f_f (expm1, 500.0L, 1.4035922178528374107397703328409120821806e+217L);
+#endif
+
 #if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
   TEST_f_f (expm1, 11356.25L, 9.05128237311923300051376115753226014206e+4931L);
 #endif
@@ -3559,9 +3566,7 @@ expm1_test (void)
   TEST_f_f (expm1, 100000.0, plus_infty, OVERFLOW_EXCEPTION);
   check_int ("errno for expm1(large) == ERANGE", errno, ERANGE, 0, 0, 0);
   TEST_f_f (expm1, max_value, plus_infty, OVERFLOW_EXCEPTION);
-#ifndef TEST_LDOUBLE /* Bug 13923.  */
   TEST_f_f (expm1, -max_value, -1);
-#endif
 
   END (expm1);
 }
diff --git a/sysdeps/i386/fpu/e_expl.S b/sysdeps/i386/fpu/e_expl.S
index 9adf2a489e..bab0a081b8 100644
--- a/sysdeps/i386/fpu/e_expl.S
+++ b/sysdeps/i386/fpu/e_expl.S
@@ -28,6 +28,10 @@
 # define IEEE754_EXPL __ieee754_exp10l
 # define EXPL_FINITE __exp10l_finite
 # define FLDLOG fldl2t
+#elif defined USE_AS_EXPM1L
+# define IEEE754_EXPL __expm1l
+# undef EXPL_FINITE
+# define FLDLOG fldl2e
 #else
 # define IEEE754_EXPL __ieee754_expl
 # define EXPL_FINITE __expl_finite
@@ -69,6 +73,12 @@ csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
 
 	.text
 ENTRY(IEEE754_EXPL)
+#ifdef USE_AS_EXPM1L
+	movzwl	4+8(%esp), %eax
+	xorb	$0x80, %ah	// invert sign bit (now 1 is "positive")
+	cmpl	$0xc006, %eax	// is num positive and exp >= 6 (number is >= 128.0)?
+	jae	HIDDEN_JUMPTARGET (__expl) // (if num is denormal, it is at least >= 64.0)
+#endif
 	fldt	4(%esp)
 /* I added the following ugly construct because expl(+-Inf) resulted
    in NaN.  The ugliness results from the bright minds at Intel.
@@ -78,7 +88,9 @@ ENTRY(IEEE754_EXPL)
 #ifdef PIC
 	LOAD_PIC_REG (cx)
 #endif
+#ifndef USE_AS_EXPM1L
 	movzwl	4+8(%esp), %eax
+#endif
 	andl	$0x7fff, %eax
 	cmpl	$0x400d, %eax
 	jle	3f
@@ -96,7 +108,16 @@ ENTRY(IEEE754_EXPL)
 	andb	$2, %ah
 	jz	3f
 	fchs
-3:	FLDLOG			/* 1  log2(base)      */
+3:
+#ifdef USE_AS_EXPM1L
+	/* Test for +-0 as argument.  */
+	fstsw	%ax
+	movb	$0x45, %dh
+	andb	%ah, %dh
+	cmpb	$0x40, %dh
+	je	2f
+#endif
+	FLDLOG			/* 1  log2(base)      */
 	fmul	%st(1), %st	/* 1  x log2(base)    */
 	frndint			/* 1  i               */
 	fld	%st(1)		/* 2  x               */
@@ -114,17 +135,39 @@ ENTRY(IEEE754_EXPL)
 	fmul	%st(4), %st	/* 4  c1 * x          */
 	faddp	%st, %st(1)	/* 3  f = f + c1 * x  */
 	f2xm1			/* 3 2^(fract(x * log2(base))) - 1 */
+#ifdef USE_AS_EXPM1L
+	fstp	%st(1)		/* 2                  */
+	fscale			/* 2 scale factor is st(1); base^x - 2^i */
+	fxch			/* 2 i                */
+	fld1			/* 3 1.0              */
+	fscale			/* 3 2^i              */
+	fld1			/* 4 1.0              */
+	fsubrp	%st, %st(1)	/* 3 2^i - 1.0        */
+	fstp	%st(1)		/* 2                  */
+	faddp	%st, %st(1)	/* 1 base^x - 1.0     */
+#else
 	fld1			/* 4 1.0              */
 	faddp			/* 3 2^(fract(x * log2(base))) */
 	fstp	%st(1)		/* 2  */
 	fscale			/* 2 scale factor is st(1); base^x */
 	fstp	%st(1)		/* 1  */
+#endif
 	fstp	%st(1)		/* 0  */
 	jmp	2f
 1:	testl	$0x200, %eax	/* Test sign.  */
 	jz	2f		/* If positive, jump.  */
 	fstp	%st
+#ifdef USE_AS_EXPM1L
+	fld1
+	fchs
+#else
 	fldz			/* Set result to 0.  */
+#endif
 2:	ret
 END(IEEE754_EXPL)
+#ifdef USE_AS_EXPM1L
+libm_hidden_def (__expm1l)
+weak_alias (__expm1l, expm1l)
+#else
 strong_alias (IEEE754_EXPL, EXPL_FINITE)
+#endif
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 63235537e0..cdc26a0c63 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -1008,10 +1008,6 @@ Test "cos_upward (9) == -0.9111302618846769883682947111811653112463":
 ildouble: 1
 ldouble: 1
 
-# cosh
-Test "cosh (0.75) == 1.29468328467684468784170818539018176":
-ildouble: 1
-
 # cosh_downward
 Test "cosh_downward (22) == 1792456423.065795780980053377632656584997":
 double: 1
@@ -1899,29 +1895,20 @@ double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 4
-ldouble: 4
+ildouble: 2
+ldouble: 2
 Test "sinh_downward (23) == 4872401723.124451299966006944252978187305":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
+ildouble: 2
+ldouble: 2
 Test "sinh_downward (24) == 13244561064.92173614705070540368454568168":
 float: 1
 ifloat: 1
-ildouble: 5
-ldouble: 5
-
-# sinh_tonearest
-Test "sinh_tonearest (22) == 1792456423.065795780701106568345764104225":
-ildouble: 3
-ldouble: 3
-Test "sinh_tonearest (23) == 4872401723.124451299966006944252978187305":
-ildouble: 1
-ldouble: 1
-Test "sinh_tonearest (24) == 13244561064.92173614705070540368454568168":
-ildouble: 6
-ldouble: 6
+ildouble: 2
+ldouble: 2
 
 # sinh_towardzero
 Test "sinh_towardzero (22) == 1792456423.065795780701106568345764104225":
@@ -1929,31 +1916,31 @@ double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 4
-ldouble: 4
+ildouble: 2
+ldouble: 2
 Test "sinh_towardzero (23) == 4872401723.124451299966006944252978187305":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
+ildouble: 2
+ldouble: 2
 Test "sinh_towardzero (24) == 13244561064.92173614705070540368454568168":
 float: 1
 ifloat: 1
-ildouble: 5
-ldouble: 5
+ildouble: 2
+ldouble: 2
 
 # sinh_upward
 Test "sinh_upward (22) == 1792456423.065795780701106568345764104225":
-ildouble: 16
-ldouble: 16
+ildouble: 1
+ldouble: 1
 Test "sinh_upward (23) == 4872401723.124451299966006944252978187305":
-ildouble: 27
-ldouble: 27
+ildouble: 1
+ldouble: 1
 Test "sinh_upward (24) == 13244561064.92173614705070540368454568168":
 double: 1
 idouble: 1
-ildouble: 7
-ldouble: 7
 
 # tan
 Test "tan (0x1p16383) == 0.422722393732022337800504160054440141575":
@@ -2587,9 +2574,6 @@ ifloat: 1
 ildouble: 1
 ldouble: 1
 
-Function: "cosh":
-ildouble: 1
-
 Function: "cosh_downward":
 double: 1
 float: 1
@@ -2862,26 +2846,22 @@ double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 5
-ldouble: 5
-
-Function: "sinh_tonearest":
-ildouble: 6
-ldouble: 6
+ildouble: 2
+ldouble: 2
 
 Function: "sinh_towardzero":
 double: 1
 float: 1
 idouble: 1
 ifloat: 1
-ildouble: 5
-ldouble: 5
+ildouble: 2
+ldouble: 2
 
 Function: "sinh_upward":
 double: 1
 idouble: 1
-ildouble: 27
-ldouble: 27
+ildouble: 1
+ldouble: 1
 
 Function: "tan":
 double: 1
diff --git a/sysdeps/i386/fpu/s_expm1l.S b/sysdeps/i386/fpu/s_expm1l.S
index e91f18b694..7fbd99b0db 100644
--- a/sysdeps/i386/fpu/s_expm1l.S
+++ b/sysdeps/i386/fpu/s_expm1l.S
@@ -1,89 +1,2 @@
-/* ix87 specific implementation of exp(x)-1.
-   Copyright (C) 1996-1997, 2002, 2005, 2008, 2012 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
-   Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>.
-   Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997.
-
-   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/>.  */
-
-	/* Using: e^x - 1 = 2^(x * log2(e)) - 1 */
-
-#include <sysdep.h>
-#include <machine/asm.h>
-
-	.section .rodata
-
-	.align ALIGNARG(4)
-	ASM_TYPE_DIRECTIVE(minus1,@object)
-minus1:	.double -1.0
-	ASM_SIZE_DIRECTIVE(minus1)
-	ASM_TYPE_DIRECTIVE(one,@object)
-one:	.double 1.0
-	ASM_SIZE_DIRECTIVE(one)
-	ASM_TYPE_DIRECTIVE(l2e,@object)
-l2e:	.tfloat 1.442695040888963407359924681002
-	ASM_SIZE_DIRECTIVE(l2e)
-
-#ifdef PIC
-#define MO(op) op##@GOTOFF(%edx)
-#else
-#define MO(op) op
-#endif
-
-	.text
-ENTRY(__expm1l)
-	movzwl	4+8(%esp), %eax	// load sign bit and 15-bit exponent
-	xorb	$0x80, %ah	// invert sign bit (now 1 is "positive")
-	cmpl	$0xc006, %eax	// is num positive and exp >= 6 (number is >= 128.0)?
-	jae	HIDDEN_JUMPTARGET (__expl) // (if num is denormal, it is at least >= 64.0)
-
-	fldt	4(%esp)		// x
-	fxam			// Is NaN or +-Inf?
-	fstsw	%ax
-	movb	$0x45, %ch
-	andb	%ah, %ch
-	cmpb	$0x40, %ch
-	je	3f		// If +-0, jump.
-#ifdef	PIC
-	LOAD_PIC_REG (dx)
-#endif
-	cmpb	$0x05, %ch
-	je	2f		// If +-Inf, jump.
-
-	fldt	MO(l2e)		// log2(e) : x
-	fmulp			// log2(e)*x
-	fld	%st		// log2(e)*x : log2(e)*x
-	frndint			// int(log2(e)*x) : log2(e)*x
-	fsubr	%st, %st(1)	// int(log2(e)*x) : fract(log2(e)*x)
-	fxch			// fract(log2(e)*x) : int(log2(e)*x)
-	f2xm1			// 2^fract(log2(e)*x)-1 : int(log2(e)*x)
-	fscale			// 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x)
-	fxch			// int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fldl	MO(one)		// 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fscale			// 2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fsubrl	MO(one)		// 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fstp	%st(1)		// 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fsubrp	%st, %st(1)	// 2^(log2(e)*x)
-	ret
-
-2:	testl	$0x200, %eax	// Test sign.
-	jz	3f		// If positive, jump.
-	fstp	%st
-	fldl	MO(minus1)	// Set result to -1.0.
-3:	ret
-END(__expm1l)
-libm_hidden_def (__expm1l)
-weak_alias (__expm1l, expm1l)
+#define USE_AS_EXPM1L
+#include <e_expl.S>
diff --git a/sysdeps/x86_64/fpu/e_expl.S b/sysdeps/x86_64/fpu/e_expl.S
index fd613f91d3..e6b842bf26 100644
--- a/sysdeps/x86_64/fpu/e_expl.S
+++ b/sysdeps/x86_64/fpu/e_expl.S
@@ -28,6 +28,10 @@
 # define IEEE754_EXPL __ieee754_exp10l
 # define EXPL_FINITE __exp10l_finite
 # define FLDLOG fldl2t
+#elif defined USE_AS_EXPM1L
+# define IEEE754_EXPL __expm1l
+# undef EXPL_FINITE
+# define FLDLOG fldl2e
 #else
 # define IEEE754_EXPL __ieee754_expl
 # define EXPL_FINITE __expl_finite
@@ -69,13 +73,21 @@ csat:	.byte 0, 0, 0, 0, 0, 0, 0, 0x80, 0x0e, 0x40
 
 	.text
 ENTRY(IEEE754_EXPL)
+#ifdef USE_AS_EXPM1L
+	movzwl	8+8(%rsp), %eax
+	xorb	$0x80, %ah	// invert sign bit (now 1 is "positive")
+	cmpl	$0xc006, %eax	// is num positive and exp >= 6 (number is >= 128.0)?
+	jae	HIDDEN_JUMPTARGET (__expl) // (if num is denormal, it is at least >= 64.0)
+#endif
 	fldt	8(%rsp)
 /* I added the following ugly construct because expl(+-Inf) resulted
    in NaN.  The ugliness results from the bright minds at Intel.
    For the i686 the code can be written better.
    -- drepper@cygnus.com.  */
 	fxam			/* Is NaN or +-Inf?  */
+#ifndef USE_AS_EXPM1L
 	movzwl	8+8(%rsp), %eax
+#endif
 	andl	$0x7fff, %eax
 	cmpl	$0x400d, %eax
 	jle	3f
@@ -93,7 +105,16 @@ ENTRY(IEEE754_EXPL)
 	andb	$2, %ah
 	jz	3f
 	fchs
-3:	FLDLOG			/* 1  log2(base)      */
+3:
+#ifdef USE_AS_EXPM1L
+	/* Test for +-0 as argument.  */
+	fstsw	%ax
+	movb	$0x45, %dh
+	andb	%ah, %dh
+	cmpb	$0x40, %dh
+	je	2f
+#endif
+	FLDLOG			/* 1  log2(base)      */
 	fmul	%st(1), %st	/* 1  x log2(base)    */
 	frndint			/* 1  i               */
 	fld	%st(1)		/* 2  x               */
@@ -111,17 +132,39 @@ ENTRY(IEEE754_EXPL)
 	fmul	%st(4), %st	/* 4  c1 * x          */
 	faddp	%st, %st(1)	/* 3  f = f + c1 * x  */
 	f2xm1			/* 3 2^(fract(x * log2(base))) - 1 */
+#ifdef USE_AS_EXPM1L
+	fstp	%st(1)		/* 2                  */
+	fscale			/* 2 scale factor is st(1); base^x - 2^i */
+	fxch			/* 2 i                */
+	fld1			/* 3 1.0              */
+	fscale			/* 3 2^i              */
+	fld1			/* 4 1.0              */
+	fsubrp	%st, %st(1)	/* 3 2^i - 1.0        */
+	fstp	%st(1)		/* 2                  */
+	faddp	%st, %st(1)	/* 1 base^x - 1.0     */
+#else
 	fld1			/* 4 1.0              */
 	faddp			/* 3 2^(fract(x * log2(base))) */
 	fstp	%st(1)		/* 2  */
 	fscale			/* 2 scale factor is st(1); base^x */
 	fstp	%st(1)		/* 1  */
+#endif
 	fstp	%st(1)		/* 0  */
 	jmp	2f
 1:	testl	$0x200, %eax	/* Test sign.  */
 	jz	2f		/* If positive, jump.  */
 	fstp	%st
+#ifdef USE_AS_EXPM1L
+	fld1
+	fchs
+#else
 	fldz			/* Set result to 0.  */
+#endif
 2:	ret
 END(IEEE754_EXPL)
+#ifdef USE_AS_EXPM1L
+libm_hidden_def (__expm1l)
+weak_alias (__expm1l, expm1l)
+#else
 strong_alias (IEEE754_EXPL, EXPL_FINITE)
+#endif
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index f33c07f326..9bb26d1fe7 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -1358,6 +1358,9 @@ ifloat: 1
 Test "expm1 (11356.25) == 9.05128237311923300051376115753226014206e+4931":
 ildouble: 1
 ldouble: 1
+Test "expm1 (500.0) == 1.4035922178528374107397703328409120821806e+217":
+double: 1
+idouble: 1
 
 # gamma
 Test "gamma (-0.5) == log(2*sqrt(pi))":
@@ -1831,62 +1834,47 @@ Test "sincos (pi/6, &sin_res, &cos_res) puts 0.866025403784438646763723170752936
 float: 1
 ifloat: 1
 
-# sinh
-Test "sinh (0x8p-32) == 1.86264514923095703232705808926175479e-9":
-ildouble: 1
-ldouble: 1
-
 # sinh_downward
 Test "sinh_downward (22) == 1792456423.065795780701106568345764104225":
 float: 1
 ifloat: 1
-ildouble: 4
-ldouble: 4
+ildouble: 2
+ldouble: 2
 Test "sinh_downward (23) == 4872401723.124451299966006944252978187305":
 float: 1
 ifloat: 1
+ildouble: 2
+ldouble: 2
 Test "sinh_downward (24) == 13244561064.92173614705070540368454568168":
 float: 1
 ifloat: 1
-ildouble: 5
-ldouble: 5
-
-# sinh_tonearest
-Test "sinh_tonearest (22) == 1792456423.065795780701106568345764104225":
-ildouble: 3
-ldouble: 3
-Test "sinh_tonearest (23) == 4872401723.124451299966006944252978187305":
-ildouble: 1
-ldouble: 1
-Test "sinh_tonearest (24) == 13244561064.92173614705070540368454568168":
-ildouble: 6
-ldouble: 6
+ildouble: 2
+ldouble: 2
 
 # sinh_towardzero
 Test "sinh_towardzero (22) == 1792456423.065795780701106568345764104225":
 float: 1
 ifloat: 1
-ildouble: 4
-ldouble: 4
+ildouble: 2
+ldouble: 2
 Test "sinh_towardzero (23) == 4872401723.124451299966006944252978187305":
 float: 1
 ifloat: 1
+ildouble: 2
+ldouble: 2
 Test "sinh_towardzero (24) == 13244561064.92173614705070540368454568168":
 float: 1
 ifloat: 1
-ildouble: 5
-ldouble: 5
+ildouble: 2
+ldouble: 2
 
 # sinh_upward
 Test "sinh_upward (22) == 1792456423.065795780701106568345764104225":
-ildouble: 16
-ldouble: 16
+ildouble: 1
+ldouble: 1
 Test "sinh_upward (23) == 4872401723.124451299966006944252978187305":
-ildouble: 27
-ldouble: 27
-Test "sinh_upward (24) == 13244561064.92173614705070540368454568168":
-ildouble: 7
-ldouble: 7
+ildouble: 1
+ldouble: 1
 
 # tan
 Test "tan (0x1p16383) == 0.422722393732022337800504160054440141575":
@@ -2733,29 +2721,21 @@ ifloat: 1
 ildouble: 1
 ldouble: 1
 
-Function: "sinh":
-ildouble: 1
-ldouble: 1
-
 Function: "sinh_downward":
 float: 1
 ifloat: 1
-ildouble: 5
-ldouble: 5
-
-Function: "sinh_tonearest":
-ildouble: 6
-ldouble: 6
+ildouble: 2
+ldouble: 2
 
 Function: "sinh_towardzero":
 float: 1
 ifloat: 1
-ildouble: 5
-ldouble: 5
+ildouble: 2
+ldouble: 2
 
 Function: "sinh_upward":
-ildouble: 27
-ldouble: 27
+ildouble: 1
+ldouble: 1
 
 Function: "tan":
 double: 1
diff --git a/sysdeps/x86_64/fpu/s_expm1l.S b/sysdeps/x86_64/fpu/s_expm1l.S
index 1380f34f01..7fbd99b0db 100644
--- a/sysdeps/x86_64/fpu/s_expm1l.S
+++ b/sysdeps/x86_64/fpu/s_expm1l.S
@@ -1,85 +1,2 @@
-/* ix87 specific implementation of exp(x)-1.
-   Copyright (C) 1996,1997,2001,2002,2008,2009,2012 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
-   Based on code by John C. Bowman <bowman@ipp-garching.mpg.de>.
-   Corrections by H.J. Lu (hjl@gnu.ai.mit.edu), 1997.
-
-   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/>.  */
-
-	/* Using: e^x - 1 = 2^(x * log2(e)) - 1 */
-
-#include <machine/asm.h>
-
-	.section .rodata
-
-	.align ALIGNARG(4)
-	ASM_TYPE_DIRECTIVE(minus1,@object)
-minus1:	.double -1.0
-	ASM_SIZE_DIRECTIVE(minus1)
-	ASM_TYPE_DIRECTIVE(one,@object)
-one:	.double 1.0
-	ASM_SIZE_DIRECTIVE(one)
-	ASM_TYPE_DIRECTIVE(l2e,@object)
-l2e:	.tfloat 1.442695040888963407359924681002
-	ASM_SIZE_DIRECTIVE(l2e)
-
-#ifdef PIC
-#define MO(op) op##(%rip)
-#else
-#define MO(op) op
-#endif
-
-	.text
-ENTRY(__expm1l)
-	movzwl	8+8(%rsp), %eax	// load sign bit and 15-bit exponent
-	xorb	$0x80, %ah	// invert sign bit (now 1 is "positive")
-	cmpl	$0xc006, %eax	// is num positive and exp >= 6 (number is >= 128.0)?
-	jae     __expl		// (if num is denormal, it is at least >= 64.0)
-
-	fldt	8(%rsp)		// x
-	fxam			// Is NaN or +-Inf?
-	fstsw	%ax
-	movb	$0x45, %ch
-	andb	%ah, %ch
-	cmpb	$0x40, %ch
-	je	3f		// If +-0, jump.
-	cmpb	$0x05, %ch
-	je	2f		// If +-Inf, jump.
-
-	fldt	MO(l2e)		// log2(e) : x
-	fmulp			// log2(e)*x
-	fld	%st		// log2(e)*x : log2(e)*x
-	frndint			// int(log2(e)*x) : log2(e)*x
-	fsubr	%st, %st(1)	// int(log2(e)*x) : fract(log2(e)*x)
-	fxch			// fract(log2(e)*x) : int(log2(e)*x)
-	f2xm1			// 2^fract(log2(e)*x)-1 : int(log2(e)*x)
-	fscale			// 2^(log2(e)*x)-2^int(log2(e)*x) : int(log2(e)*x)
-	fxch			// int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fldl	MO(one)		// 1 : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fscale			// 2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fsubrl	MO(one)		// 1-2^int(log2(e)*x) : int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fstp	%st(1)		// 1-2^int(log2(e)*x) : 2^(log2(e)*x)-2^int(log2(e)*x)
-	fsubrp	%st, %st(1)	// 2^(log2(e)*x)-1
-	ret
-
-2:	testl	$0x200, %eax	// Test sign.
-	jz	3f		// If positive, jump.
-	fstp	%st
-	fldl	MO(minus1)	// Set result to -1.0.
-3:	ret
-END(__expm1l)
-libm_hidden_def (__expm1l)
-weak_alias (__expm1l, expm1l)
+#define USE_AS_EXPM1L
+#include <e_expl.S>