Remove inaccurate x86 cexp implementations (bug 13883).

5 files changed, 88 insertions, 786 deletions

diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 1d87514e9b..4d61635f25 100644
--- a/sysdeps/i386/fpu/libm-test-ulps
+++ b/sysdeps/i386/fpu/libm-test-ulps
@@ -431,15 +431,44 @@ idouble: 1
 ifloat: 1
 
 # cexp
+Test "Real part of: cexp (-10000 + 0x1p16383 i) == 1.045876464564882298442774542991176546722e-4343 + 4.421154026488516836023811173959413420548e-4344 i":
+ildouble: 1
+ldouble: 1
 Test "Real part of: cexp (-2.0 - 3.0 i) == -0.13398091492954261346140525546115575 - 0.019098516261135196432576240858800925 i":
 ildouble: 1
 ldouble: 1
 Test "Imaginary part of: cexp (-2.0 - 3.0 i) == -0.13398091492954261346140525546115575 - 0.019098516261135196432576240858800925 i":
+float: 1
+ifloat: 1
 ildouble: 1
 ldouble: 1
+Test "Imaginary part of: cexp (0 + 0x1p65 i) == 0.99888622066058013610642172179340364209972 - 0.047183876212354673805106149805700013943218 i":
+float: 1
+ifloat: 1
+Test "Imaginary part of: cexp (0 - 0x1p65 i) == 0.99888622066058013610642172179340364209972 + 0.047183876212354673805106149805700013943218 i":
+float: 1
+ifloat: 1
 Test "Real part of: cexp (0.75 + 1.25 i) == 0.667537446429131586942201977015932112 + 2.00900045494094876258347228145863909 i":
+float: 1
+ifloat: 1
 ildouble: 1
 ldouble: 1
+Test "Imaginary part of: cexp (0.75 + 1.25 i) == 0.667537446429131586942201977015932112 + 2.00900045494094876258347228145863909 i":
+ildouble: 1
+ldouble: 1
+Test "Real part of: cexp (50 + 0x1p127 i) == 4.053997150228616856622417636046265337193e21 + 3.232070315463388524466674772633810238819e21 i":
+double: 2
+idouble: 2
+Test "Imaginary part of: cexp (50 + 0x1p127 i) == 4.053997150228616856622417636046265337193e21 + 3.232070315463388524466674772633810238819e21 i":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
+ildouble: 1
+ldouble: 1
+Test "Real part of: cexp (500 + 0x1p1023 i) == -1.159886268932754433233243794561351783426e217 + 7.904017694554466595359379965081774849708e216 i":
+double: 1
+idouble: 1
 
 # clog
 Test "Real part of: clog (0.75 + 1.25 i) == 0.376885901188190075998919126749298416 + 1.03037682652431246378774332703115153 i":
@@ -815,18 +844,22 @@ ifloat: 1
 ildouble: 1
 ldouble: 1
 Test "Imaginary part of: cpow (0.75 + 1.25 i, 0.0 + 1.0 i) == 0.331825439177608832276067945276730566 + 0.131338600281188544930936345230903032 i":
+float: 1
+ifloat: 1
 ildouble: 1
 ldouble: 1
 Test "Real part of: cpow (0.75 + 1.25 i, 0.75 + 1.25 i) == 0.117506293914473555420279832210420483 + 0.346552747708338676483025352060418001 i":
-float: 3
-ifloat: 3
+double: 1
+float: 4
+idouble: 1
+ifloat: 4
 ildouble: 6
 ldouble: 6
 Test "Imaginary part of: cpow (0.75 + 1.25 i, 0.75 + 1.25 i) == 0.117506293914473555420279832210420483 + 0.346552747708338676483025352060418001 i":
 float: 1
 ifloat: 1
-ildouble: 1
-ldouble: 1
+ildouble: 2
+ldouble: 2
 Test "Real part of: cpow (0.75 + 1.25 i, 1.0 + 0.0 i) == 0.75 + 1.25 i":
 ildouble: 1
 ldouble: 1
@@ -834,22 +867,27 @@ Test "Imaginary part of: cpow (0.75 + 1.25 i, 1.0 + 0.0 i) == 0.75 + 1.25 i":
 float: 1
 ifloat: 1
 Test "Real part of: cpow (0.75 + 1.25 i, 1.0 + 1.0 i) == 0.0846958290317209430433805274189191353 + 0.513285749182902449043287190519090481 i":
-double: 1
+double: 2
 float: 3
-idouble: 1
+idouble: 2
 ifloat: 3
 ildouble: 3
 ldouble: 3
+Test "Real part of: cpow (2 + 0 i, 10 + 0 i) == 1024.0 + 0.0 i":
+ildouble: 1
+ldouble: 1
 Test "Real part of: cpow (2 + 3 i, 4 + 0 i) == -119.0 - 120.0 i":
 double: 1
-float: 4
+float: 5
 idouble: 1
-ifloat: 4
+ifloat: 5
+ildouble: 1
+ldouble: 1
 Test "Imaginary part of: cpow (2 + 3 i, 4 + 0 i) == -119.0 - 120.0 i":
-float: 1
-ifloat: 1
-ildouble: 2
-ldouble: 2
+float: 2
+ifloat: 2
+ildouble: 4
+ldouble: 4
 Test "Imaginary part of: cpow (e + 0 i, 0 + 2 * M_PIl i) == 1.0 + 0.0 i":
 double: 2
 float: 3
@@ -2124,10 +2162,18 @@ ildouble: 1
 ldouble: 1
 
 Function: Real part of "cexp":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
 ildouble: 1
 ldouble: 1
 
 Function: Imaginary part of "cexp":
+double: 1
+float: 1
+idouble: 1
+ifloat: 1
 ildouble: 1
 ldouble: 1
 
@@ -2212,20 +2258,20 @@ double: 1
 ldouble: 1
 
 Function: Real part of "cpow":
-double: 1
-float: 4
-idouble: 1
-ifloat: 4
-ildouble: 6
-ldouble: 6
+double: 2
+float: 5
+idouble: 2
+ifloat: 5
+ildouble: 5
+ldouble: 5
 
 Function: Imaginary part of "cpow":
 double: 2
 float: 3
 idouble: 2
 ifloat: 3
-ildouble: 2
-ldouble: 2
+ildouble: 4
+ldouble: 4
 
 Function: Real part of "csin":
 float: 1
diff --git a/sysdeps/i386/fpu/s_cexp.S b/sysdeps/i386/fpu/s_cexp.S
deleted file mode 100644
index e5fdb7d735..0000000000
--- a/sysdeps/i386/fpu/s_cexp.S
+++ /dev/null
@@ -1,253 +0,0 @@
-/* ix87 specific implementation of complex exponential function for double.
-   Copyright (C) 1997, 2005, 2012 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 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/>.  */
-
-#include <sysdep.h>
-
-	.section .rodata
-
-	.align ALIGNARG(4)
-	ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
-huge_nan_null_null:
-	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
-	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
-	.double	0.0
-zero:	.double	0.0
-infinity:
-	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
-	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
-	.double 0.0
-	.byte 0, 0, 0, 0, 0, 0, 0, 0x80
-	ASM_SIZE_DIRECTIVE(huge_nan_null_null)
-
-	ASM_TYPE_DIRECTIVE(twopi,@object)
-twopi:
-	.byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
-	.byte 0, 0, 0, 0, 0, 0
-	ASM_SIZE_DIRECTIVE(twopi)
-
-	ASM_TYPE_DIRECTIVE(l2e,@object)
-l2e:
-	.byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
-	.byte 0, 0, 0, 0, 0, 0
-	ASM_SIZE_DIRECTIVE(l2e)
-
-	ASM_TYPE_DIRECTIVE(one,@object)
-one:	.double 1.0
-	ASM_SIZE_DIRECTIVE(one)
-
-
-#ifdef PIC
-#define MO(op) op##@GOTOFF(%ecx)
-#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
-#else
-#define MO(op) op
-#define MOX(op,x,f) op(,x,f)
-#endif
-
-	.text
-ENTRY(__cexp)
-	fldl	8(%esp)			/* x */
-	fxam
-	fnstsw
-	fldl	16(%esp)		/* y : x */
-#ifdef  PIC
-	LOAD_PIC_REG (cx)
-#endif
-	movb	%ah, %dh
-	andb	$0x45, %ah
-	cmpb	$0x05, %ah
-	je	1f			/* Jump if real part is +-Inf */
-	cmpb	$0x01, %ah
-	je	2f			/* Jump if real part is NaN */
-
-	fxam				/* y : x */
-	fnstsw
-	/* If the imaginary part is not finite we return NaN+i NaN, as
-	   for the case when the real part is NaN.  A test for +-Inf and
-	   NaN would be necessary.  But since we know the stack register
-	   we applied `fxam' to is not empty we can simply use one test.
-	   Check your FPU manual for more information.  */
-	andb	$0x01, %ah
-	cmpb	$0x01, %ah
-	je	20f
-
-	/* We have finite numbers in the real and imaginary part.  Do
-	   the real work now.  */
-	fxch			/* x : y */
-	fldt	MO(l2e)		/* log2(e) : x : y */
-	fmulp			/* x * log2(e) : y */
-	fld	%st		/* x * log2(e) : x * log2(e) : y */
-	frndint			/* int(x * log2(e)) : x * log2(e) : y */
-	fsubr	%st, %st(1)	/* int(x * log2(e)) : frac(x * log2(e)) : y */
-	fxch			/* frac(x * log2(e)) : int(x * log2(e)) : y */
-	f2xm1			/* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
-	faddl	MO(one)		/* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
-	fscale			/* e^x : int(x * log2(e)) : y */
-	fst	%st(1)		/* e^x : e^x : y */
-	fxch	%st(2)		/* y : e^x : e^x */
-	fsincos			/* cos(y) : sin(y) : e^x : e^x */
-	fnstsw
-	testl	$0x400, %eax
-	jnz	7f
-	fmulp	%st, %st(3)	/* sin(y) : e^x : e^x * cos(y) */
-	fmulp	%st, %st(1)	/* e^x * sin(y) : e^x * cos(y) */
-	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-	fstpl	8(%eax)
-	fstpl	(%eax)
-	ret	$4
-
-	/* We have to reduce the argument to fsincos.  */
-	.align ALIGNARG(4)
-7:	fldt	MO(twopi)	/* 2*pi : y : e^x : e^x */
-	fxch			/* y : 2*pi : e^x : e^x */
-8:	fprem1			/* y%(2*pi) : 2*pi : e^x : e^x */
-	fnstsw
-	testl	$0x400, %eax
-	jnz	8b
-	fstp	%st(1)		/* y%(2*pi) : e^x : e^x */
-	fsincos			/* cos(y) : sin(y) : e^x : e^x */
-	fmulp	%st, %st(3)
-	fmulp	%st, %st(1)
-	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-	fstpl	8(%eax)
-	fstpl	(%eax)
-	ret	$4
-
-	/* The real part is +-inf.  We must make further differences.  */
-	.align ALIGNARG(4)
-1:	fxam			/* y : x */
-	fnstsw
-	movb	%ah, %dl
-	testb	$0x01, %ah	/* See above why 0x01 is usable here.  */
-	jne	3f
-
-
-	/* The real part is +-Inf and the imaginary part is finite.  */
-	andl	$0x245, %edx
-	cmpb	$0x40, %dl	/* Imaginary part == 0?  */
-	je	4f		/* Yes ->  */
-
-	fxch			/* x : y */
-	shrl	$5, %edx
-	fstp	%st(0)		/* y */ /* Drop the real part.  */
-	andl	$16, %edx	/* This puts the sign bit of the real part
-				   in bit 4.  So we can use it to index a
-				   small array to select 0 or Inf.  */
-	fsincos			/* cos(y) : sin(y) */
-	fnstsw
-	testl	$0x0400, %eax
-	jnz	5f
-	fldl	MOX(huge_nan_null_null,%edx,1)
-	movl	4(%esp), %edx		/* Pointer to memory for result.  */
-	fstl	8(%edx)
-	fstpl	(%edx)
-	ftst
-	fnstsw
-	shll	$23, %eax
-	andl	$0x80000000, %eax
-	orl	%eax, 4(%edx)
-	fstp	%st(0)
-	ftst
-	fnstsw
-	shll	$23, %eax
-	andl	$0x80000000, %eax
-	orl	%eax, 12(%edx)
-	fstp	%st(0)
-	ret	$4
-	/* We must reduce the argument to fsincos.  */
-	.align ALIGNARG(4)
-5:	fldt	MO(twopi)
-	fxch
-6:	fprem1
-	fnstsw
-	testl	$0x400, %eax
-	jnz	6b
-	fstp	%st(1)
-	fsincos
-	fldl	MOX(huge_nan_null_null,%edx,1)
-	movl	4(%esp), %edx		/* Pointer to memory for result.  */
-	fstl	8(%edx)
-	fstpl	(%edx)
-	ftst
-	fnstsw
-	shll	$23, %eax
-	andl	$0x80000000, %eax
-	orl	%eax, 4(%edx)
-	fstp	%st(0)
-	ftst
-	fnstsw
-	shll	$23, %eax
-	andl	$0x80000000, %eax
-	orl	%eax, 12(%edx)
-	fstp	%st(0)
-	ret	$4
-
-	/* The real part is +-Inf and the imaginary part is +-0.  So return
-	   +-Inf+-0i.  */
-	.align ALIGNARG(4)
-4:	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-	fstpl	8(%eax)
-	shrl	$5, %edx
-	fstp	%st(0)
-	andl	$16, %edx
-	fldl	MOX(huge_nan_null_null,%edx,1)
-	fstpl	(%eax)
-	ret	$4
-
-	/* The real part is +-Inf, the imaginary is also is not finite.  */
-	.align ALIGNARG(4)
-3:	fstp	%st(0)
-	fstp	%st(0)		/* <empty> */
-	andb	$0x45, %ah
-	andb	$0x47, %dh
-	xorb	%dh, %ah
-	jnz	30f
-	fldl	MO(infinity)	/* Raise invalid exception.  */
-	fmull	MO(zero)
-	fstp	%st(0)
-30:	movl	%edx, %eax
-	shrl	$5, %edx
-	shll	$4, %eax
-	andl	$16, %edx
-	andl	$32, %eax
-	orl	%eax, %edx
-	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-
-	fldl	MOX(huge_nan_null_null,%edx,1)
-	fldl	MOX(huge_nan_null_null+8,%edx,1)
-	fxch
-	fstpl	(%eax)
-	fstpl	8(%eax)
-	ret	$4
-
-	/* The real part is NaN.  */
-	.align ALIGNARG(4)
-20:	fldl	MO(infinity)		/* Raise invalid exception.  */
-	fmull	MO(zero)
-	fstp	%st(0)
-2:	fstp	%st(0)
-	fstp	%st(0)
-	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-	fldl	MO(huge_nan_null_null+8)
-	fstl	(%eax)
-	fstpl	8(%eax)
-	ret	$4
-
-END(__cexp)
-weak_alias (__cexp, cexp)
diff --git a/sysdeps/i386/fpu/s_cexpf.S b/sysdeps/i386/fpu/s_cexpf.S
deleted file mode 100644
index 6ed66e6d04..0000000000
--- a/sysdeps/i386/fpu/s_cexpf.S
+++ /dev/null
@@ -1,257 +0,0 @@
-/* ix87 specific implementation of complex exponential function for double.
-   Copyright (C) 1997, 2005, 2012 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 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/>.  */
-
-#include <sysdep.h>
-
-	.section .rodata
-
-	.align ALIGNARG(4)
-	ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
-huge_nan_null_null:
-	.byte 0, 0, 0x80, 0x7f
-	.byte 0, 0, 0xc0, 0x7f
-	.float	0.0
-zero:	.float	0.0
-infinity:
-	.byte 0, 0, 0x80, 0x7f
-	.byte 0, 0, 0xc0, 0x7f
-	.float 0.0
-	.byte 0, 0, 0, 0x80
-	ASM_SIZE_DIRECTIVE(huge_nan_null_null)
-
-	ASM_TYPE_DIRECTIVE(twopi,@object)
-twopi:
-	.byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
-	.byte 0, 0, 0, 0, 0, 0
-	ASM_SIZE_DIRECTIVE(twopi)
-
-	ASM_TYPE_DIRECTIVE(l2e,@object)
-l2e:
-	.byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
-	.byte 0, 0, 0, 0, 0, 0
-	ASM_SIZE_DIRECTIVE(l2e)
-
-	ASM_TYPE_DIRECTIVE(one,@object)
-one:	.double 1.0
-	ASM_SIZE_DIRECTIVE(one)
-
-
-#ifdef PIC
-#define MO(op) op##@GOTOFF(%ecx)
-#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
-#else
-#define MO(op) op
-#define MOX(op,x,f) op(,x,f)
-#endif
-
-	.text
-ENTRY(__cexpf)
-	flds	4(%esp)			/* x */
-	fxam
-	fnstsw
-	flds	8(%esp)			/* y : x */
-#ifdef  PIC
-	LOAD_PIC_REG (cx)
-#endif
-	movb	%ah, %dh
-	andb	$0x45, %ah
-	cmpb	$0x05, %ah
-	je	1f			/* Jump if real part is +-Inf */
-	cmpb	$0x01, %ah
-	je	2f			/* Jump if real part is NaN */
-
-	fxam				/* y : x */
-	fnstsw
-	/* If the imaginary part is not finite we return NaN+i NaN, as
-	   for the case when the real part is NaN.  A test for +-Inf and
-	   NaN would be necessary.  But since we know the stack register
-	   we applied `fxam' to is not empty we can simply use one test.
-	   Check your FPU manual for more information.  */
-	andb	$0x01, %ah
-	cmpb	$0x01, %ah
-	je	20f
-
-	/* We have finite numbers in the real and imaginary part.  Do
-	   the real work now.  */
-	fxch			/* x : y */
-	fldt	MO(l2e)		/* log2(e) : x : y */
-	fmulp			/* x * log2(e) : y */
-	fld	%st		/* x * log2(e) : x * log2(e) : y */
-	frndint			/* int(x * log2(e)) : x * log2(e) : y */
-	fsubr	%st, %st(1)	/* int(x * log2(e)) : frac(x * log2(e)) : y */
-	fxch			/* frac(x * log2(e)) : int(x * log2(e)) : y */
-	f2xm1			/* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
-	faddl	MO(one)		/* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
-	fscale			/* e^x : int(x * log2(e)) : y */
-	fst	%st(1)		/* e^x : e^x : y */
-	fxch	%st(2)		/* y : e^x : e^x */
-	fsincos			/* cos(y) : sin(y) : e^x : e^x */
-	fnstsw
-	testl	$0x400, %eax
-	jnz	7f
-	fmulp	%st, %st(3)	/* sin(y) : e^x : e^x * cos(y) */
-	fmulp	%st, %st(1)	/* e^x * sin(y) : e^x * cos(y) */
-	subl	$8, %esp
-	cfi_adjust_cfa_offset (8)
-	fstps	4(%esp)
-	fstps	(%esp)
-	popl	%eax
-	cfi_adjust_cfa_offset (-4)
-	popl	%edx
-	cfi_adjust_cfa_offset (-4)
-	ret
-
-	/* We have to reduce the argument to fsincos.  */
-	.align ALIGNARG(4)
-7:	fldt	MO(twopi)	/* 2*pi : y : e^x : e^x */
-	fxch			/* y : 2*pi : e^x : e^x */
-8:	fprem1			/* y%(2*pi) : 2*pi : e^x : e^x */
-	fnstsw
-	testl	$0x400, %eax
-	jnz	8b
-	fstp	%st(1)		/* y%(2*pi) : e^x : e^x */
-	fsincos			/* cos(y) : sin(y) : e^x : e^x */
-	fmulp	%st, %st(3)
-	fmulp	%st, %st(1)
-	subl	$8, %esp
-	cfi_adjust_cfa_offset (8)
-	fstps	4(%esp)
-	fstps	(%esp)
-	popl	%eax
-	cfi_adjust_cfa_offset (-4)
-	popl	%edx
-	cfi_adjust_cfa_offset (-4)
-	ret
-
-	/* The real part is +-inf.  We must make further differences.  */
-	.align ALIGNARG(4)
-1:	fxam			/* y : x */
-	fnstsw
-	movb	%ah, %dl
-	testb	$0x01, %ah	/* See above why 0x01 is usable here.  */
-	jne	3f
-
-
-	/* The real part is +-Inf and the imaginary part is finite.  */
-	andl	$0x245, %edx
-	cmpb	$0x40, %dl	/* Imaginary part == 0?  */
-	je	4f		/* Yes ->  */
-
-	fxch			/* x : y */
-	shrl	$6, %edx
-	fstp	%st(0)		/* y */ /* Drop the real part.  */
-	andl	$8, %edx	/* This puts the sign bit of the real part
-				   in bit 3.  So we can use it to index a
-				   small array to select 0 or Inf.  */
-	fsincos			/* cos(y) : sin(y) */
-	fnstsw
-	testl	$0x0400, %eax
-	jnz	5f
-	fxch
-	ftst
-	fnstsw
-	fstp	%st(0)
-	shll	$23, %eax
-	andl	$0x80000000, %eax
-	orl	MOX(huge_nan_null_null,%edx,1), %eax
-	movl	MOX(huge_nan_null_null,%edx,1), %ecx
-	movl	%eax, %edx
-	ftst
-	fnstsw
-	fstp	%st(0)
-	shll	$23, %eax
-	andl	$0x80000000, %eax
-	orl	%ecx, %eax
-	ret
-	/* We must reduce the argument to fsincos.  */
-	.align ALIGNARG(4)
-5:	fldt	MO(twopi)
-	fxch
-6:	fprem1
-	fnstsw
-	testl	$0x400, %eax
-	jnz	6b
-	fstp	%st(1)
-	fsincos
-	fxch
-	ftst
-	fnstsw
-	fstp	%st(0)
-	shll	$23, %eax
-	andl	$0x80000000, %eax
-	orl	MOX(huge_nan_null_null,%edx,1), %eax
-	movl	MOX(huge_nan_null_null,%edx,1), %ecx
-	movl	%eax, %edx
-	ftst
-	fnstsw
-	fstp	%st(0)
-	shll	$23, %eax
-	andl	$0x80000000, %eax
-	orl	%ecx, %eax
-	ret
-
-	/* The real part is +-Inf and the imaginary part is +-0.  So return
-	   +-Inf+-0i.  */
-	.align ALIGNARG(4)
-4:	subl	$4, %esp
-	cfi_adjust_cfa_offset (4)
-	fstps	(%esp)
-	shrl	$6, %edx
-	fstp	%st(0)
-	andl	$8, %edx
-	movl	MOX(huge_nan_null_null,%edx,1), %eax
-	popl	%edx
-	cfi_adjust_cfa_offset (-4)
-	ret
-
-	/* The real part is +-Inf, the imaginary is also is not finite.  */
-	.align ALIGNARG(4)
-3:	fstp	%st(0)
-	fstp	%st(0)		/* <empty> */
-	andb	$0x45, %ah
-	andb	$0x47, %dh
-	xorb	%dh, %ah
-	jnz	30f
-	flds	MO(infinity)	/* Raise invalid exception.  */
-	fmuls	MO(zero)
-	fstp	%st(0)
-30:	movl	%edx, %eax
-	shrl	$6, %edx
-	shll	$3, %eax
-	andl	$8, %edx
-	andl	$16, %eax
-	orl	%eax, %edx
-
-	movl	MOX(huge_nan_null_null,%edx,1), %eax
-	movl	MOX(huge_nan_null_null+4,%edx,1), %edx
-	ret
-
-	/* The real part is NaN.  */
-	.align ALIGNARG(4)
-20:	flds	MO(infinity)		/* Raise invalid exception.  */
-	fmuls	MO(zero)
-	fstp	%st(0)
-2:	fstp	%st(0)
-	fstp	%st(0)
-	movl	MO(huge_nan_null_null+4), %eax
-	movl	%eax, %edx
-	ret
-
-END(__cexpf)
-weak_alias (__cexpf, cexpf)
diff --git a/sysdeps/i386/fpu/s_cexpl.S b/sysdeps/i386/fpu/s_cexpl.S
deleted file mode 100644
index ab02a172ad..0000000000
--- a/sysdeps/i386/fpu/s_cexpl.S
+++ /dev/null
@@ -1,256 +0,0 @@
-/* ix87 specific implementation of complex exponential function for double.
-   Copyright (C) 1997, 2005, 2012 Free Software Foundation, Inc.
-   This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 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/>.  */
-
-#include <sysdep.h>
-
-	.section .rodata
-
-	.align ALIGNARG(4)
-	ASM_TYPE_DIRECTIVE(huge_nan_null_null,@object)
-huge_nan_null_null:
-	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
-	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
-	.double	0.0
-zero:	.double	0.0
-infinity:
-	.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
-	.byte 0, 0, 0, 0, 0, 0, 0xff, 0x7f
-	.double 0.0
-	.byte 0, 0, 0, 0, 0, 0, 0, 0x80
-	ASM_SIZE_DIRECTIVE(huge_nan_null_null)
-
-	ASM_TYPE_DIRECTIVE(twopi,@object)
-twopi:
-	.byte 0x35, 0xc2, 0x68, 0x21, 0xa2, 0xda, 0xf, 0xc9, 0x1, 0x40
-	.byte 0, 0, 0, 0, 0, 0
-	ASM_SIZE_DIRECTIVE(twopi)
-
-	ASM_TYPE_DIRECTIVE(l2e,@object)
-l2e:
-	.byte 0xbc, 0xf0, 0x17, 0x5c, 0x29, 0x3b, 0xaa, 0xb8, 0xff, 0x3f
-	.byte 0, 0, 0, 0, 0, 0
-	ASM_SIZE_DIRECTIVE(l2e)
-
-	ASM_TYPE_DIRECTIVE(one,@object)
-one:	.double 1.0
-	ASM_SIZE_DIRECTIVE(one)
-
-
-#ifdef PIC
-#define MO(op) op##@GOTOFF(%ecx)
-#define MOX(op,x,f) op##@GOTOFF(%ecx,x,f)
-#else
-#define MO(op) op
-#define MOX(op,x,f) op(,x,f)
-#endif
-
-	.text
-ENTRY(__cexpl)
-	fldt	8(%esp)			/* x */
-	fxam
-	fnstsw
-	fldt	20(%esp)		/* y : x */
-#ifdef  PIC
-	LOAD_PIC_REG (cx)
-#endif
-	movb	%ah, %dh
-	andb	$0x45, %ah
-	cmpb	$0x05, %ah
-	je	1f			/* Jump if real part is +-Inf */
-	cmpb	$0x01, %ah
-	je	2f			/* Jump if real part is NaN */
-
-	fxam				/* y : x */
-	fnstsw
-	/* If the imaginary part is not finite we return NaN+i NaN, as
-	   for the case when the real part is NaN.  A test for +-Inf and
-	   NaN would be necessary.  But since we know the stack register
-	   we applied `fxam' to is not empty we can simply use one test.
-	   Check your FPU manual for more information.  */
-	andb	$0x01, %ah
-	cmpb	$0x01, %ah
-	je	20f
-
-	/* We have finite numbers in the real and imaginary part.  Do
-	   the real work now.  */
-	fxch			/* x : y */
-	fldt	MO(l2e)		/* log2(e) : x : y */
-	fmulp			/* x * log2(e) : y */
-	fld	%st		/* x * log2(e) : x * log2(e) : y */
-	frndint			/* int(x * log2(e)) : x * log2(e) : y */
-	fsubr	%st, %st(1)	/* int(x * log2(e)) : frac(x * log2(e)) : y */
-	fxch			/* frac(x * log2(e)) : int(x * log2(e)) : y */
-	f2xm1			/* 2^frac(x * log2(e))-1 : int(x * log2(e)) : y */
-	faddl	MO(one)		/* 2^frac(x * log2(e)) : int(x * log2(e)) : y */
-	fscale			/* e^x : int(x * log2(e)) : y */
-	fst	%st(1)		/* e^x : e^x : y */
-	fxch	%st(2)		/* y : e^x : e^x */
-	fsincos			/* cos(y) : sin(y) : e^x : e^x */
-	fnstsw
-	testl	$0x400, %eax
-	jnz	7f
-	fmulp	%st, %st(3)	/* sin(y) : e^x : e^x * cos(y) */
-	fmulp	%st, %st(1)	/* e^x * sin(y) : e^x * cos(y) */
-	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-	fstpt	12(%eax)
-	fstpt	(%eax)
-	ret	$4
-
-	/* We have to reduce the argument to fsincos.  */
-	.align ALIGNARG(4)
-7:	fldt	MO(twopi)	/* 2*pi : y : e^x : e^x */
-	fxch			/* y : 2*pi : e^x : e^x */
-8:	fprem1			/* y%(2*pi) : 2*pi : e^x : e^x */
-	fnstsw
-	testl	$0x400, %eax
-	jnz	8b
-	fstp	%st(1)		/* y%(2*pi) : e^x : e^x */
-	fsincos			/* cos(y) : sin(y) : e^x : e^x */
-	fmulp	%st, %st(3)
-	fmulp	%st, %st(1)
-	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-	fstpt	12(%eax)
-	fstpt	(%eax)
-	ret	$4
-
-	/* The real part is +-inf.  We must make further differences.  */
-	.align ALIGNARG(4)
-1:	fxam			/* y : x */
-	fnstsw
-	movb	%ah, %dl
-	testb	$0x01, %ah	/* See above why 0x01 is usable here.  */
-	jne	3f
-
-
-	/* The real part is +-Inf and the imaginary part is finite.  */
-	andl	$0x245, %edx
-	cmpb	$0x40, %dl	/* Imaginary part == 0?  */
-	je	4f		/* Yes ->  */
-
-	fxch			/* x : y */
-	shrl	$5, %edx
-	fstp	%st(0)		/* y */ /* Drop the real part.  */
-	andl	$16, %edx	/* This puts the sign bit of the real part
-				   in bit 4.  So we can use it to index a
-				   small array to select 0 or Inf.  */
-	fsincos			/* cos(y) : sin(y) */
-	fnstsw
-	testl	$0x0400, %eax
-	jnz	5f
-	fldl	MOX(huge_nan_null_null,%edx,1)
-	movl	4(%esp), %edx		/* Pointer to memory for result.  */
-	fld	%st
-	fstpt	12(%edx)
-	fstpt	(%edx)
-	ftst
-	fnstsw
-	shll	$7, %eax
-	andl	$0x8000, %eax
-	orl	%eax, 8(%edx)
-	fstp	%st(0)
-	ftst
-	fnstsw
-	shll	$7, %eax
-	andl	$0x8000, %eax
-	orl	%eax, 20(%edx)
-	fstp	%st(0)
-	ret	$4
-	/* We must reduce the argument to fsincos.  */
-	.align ALIGNARG(4)
-5:	fldt	MO(twopi)
-	fxch
-6:	fprem1
-	fnstsw
-	testl	$0x400, %eax
-	jnz	6b
-	fstp	%st(1)
-	fsincos
-	fldl	MOX(huge_nan_null_null,%edx,1)
-	movl	4(%esp), %edx		/* Pointer to memory for result.  */
-	fld	%st
-	fstpt	12(%edx)
-	fstpt	(%edx)
-	ftst
-	fnstsw
-	shll	$7, %eax
-	andl	$0x8000, %eax
-	orl	%eax, 8(%edx)
-	fstp	%st(0)
-	ftst
-	fnstsw
-	shll	$7, %eax
-	andl	$0x8000, %eax
-	orl	%eax, 20(%edx)
-	fstp	%st(0)
-	ret	$4
-
-	/* The real part is +-Inf and the imaginary part is +-0.  So return
-	   +-Inf+-0i.  */
-	.align ALIGNARG(4)
-4:	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-	fstpt	12(%eax)
-	shrl	$5, %edx
-	fstp	%st(0)
-	andl	$16, %edx
-	fldl	MOX(huge_nan_null_null,%edx,1)
-	fstpt	(%eax)
-	ret	$4
-
-	/* The real part is +-Inf, the imaginary is also is not finite.  */
-	.align ALIGNARG(4)
-3:	fstp	%st(0)
-	fstp	%st(0)		/* <empty> */
-	andb	$0x45, %ah
-	andb	$0x47, %dh
-	xorb	%dh, %ah
-	jnz	30f
-	fldl	MO(infinity)	/* Raise invalid exception.  */
-	fmull	MO(zero)
-	fstp	%st(0)
-30:	movl	%edx, %eax
-	shrl	$5, %edx
-	shll	$4, %eax
-	andl	$16, %edx
-	andl	$32, %eax
-	orl	%eax, %edx
-	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-
-	fldl	MOX(huge_nan_null_null,%edx,1)
-	fldl	MOX(huge_nan_null_null+8,%edx,1)
-	fxch
-	fstpt	(%eax)
-	fstpt	12(%eax)
-	ret	$4
-
-	/* The real part is NaN.  */
-	.align ALIGNARG(4)
-20:	fldl	MO(infinity)		/* Raise invalid exception.  */
-	fmull	MO(zero)
-	fstp	%st(0)
-2:	fstp	%st(0)
-	fstp	%st(0)
-	movl	4(%esp), %eax		/* Pointer to memory for result.  */
-	fldl	MO(huge_nan_null_null+8)
-	fld	%st(0)
-	fstpt	(%eax)
-	fstpt	12(%eax)
-	ret	$4
-
-END(__cexpl)
-weak_alias (__cexpl, cexpl)
diff --git a/sysdeps/x86_64/fpu/libm-test-ulps b/sysdeps/x86_64/fpu/libm-test-ulps
index fef6ea1a8a..e5eb8f9379 100644
--- a/sysdeps/x86_64/fpu/libm-test-ulps
+++ b/sysdeps/x86_64/fpu/libm-test-ulps
@@ -482,6 +482,9 @@ float: 1
 ifloat: 1
 
 # cexp
+Test "Real part of: cexp (-10000 + 0x1p16383 i) == 1.045876464564882298442774542991176546722e-4343 + 4.421154026488516836023811173959413420548e-4344 i":
+ildouble: 1
+ldouble: 1
 Test "Imaginary part of: cexp (-2.0 - 3.0 i) == -0.13398091492954261346140525546115575 - 0.019098516261135196432576240858800925 i":
 float: 1
 ifloat: 1
@@ -491,6 +494,19 @@ ifloat: 1
 Test "Imaginary part of: cexp (0.75 + 1.25 i) == 0.667537446429131586942201977015932112 + 2.00900045494094876258347228145863909 i":
 ildouble: 1
 ldouble: 1
+Test "Real part of: cexp (50 + 0x1p127 i) == 4.053997150228616856622417636046265337193e21 + 3.232070315463388524466674772633810238819e21 i":
+double: 2
+float: 1
+idouble: 2
+ifloat: 1
+Test "Imaginary part of: cexp (50 + 0x1p127 i) == 4.053997150228616856622417636046265337193e21 + 3.232070315463388524466674772633810238819e21 i":
+double: 1
+idouble: 1
+ildouble: 1
+ldouble: 1
+Test "Real part of: cexp (500 + 0x1p1023 i) == -1.159886268932754433233243794561351783426e217 + 7.904017694554466595359379965081774849708e216 i":
+double: 1
+idouble: 1
 
 # clog
 Test "Imaginary part of: clog (-2 - 3 i) == 1.2824746787307683680267437207826593 - 2.1587989303424641704769327722648368 i":
@@ -2090,11 +2106,17 @@ ildouble: 1
 ldouble: 1
 
 Function: Real part of "cexp":
+double: 2
 float: 1
+idouble: 2
 ifloat: 1
+ildouble: 1
+ldouble: 1
 
 Function: Imaginary part of "cexp":
+double: 1
 float: 1
+idouble: 1
 ifloat: 1
 ildouble: 1
 ldouble: 1