Remove inaccurate x86 cexp implementations (bug 13883).

8 files changed, 114 insertions, 787 deletions

diff --git a/ChangeLog b/ChangeLog
index 75dbb36..ebffe4f 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,13 @@
+2012-03-21  Joseph Myers  <joseph@codesourcery.com>
+
+	[BZ #13883]
+	* sysdeps/i386/fpu/s_cexp.S: Remove.
+	* sysdeps/i386/fpu/s_cexpf.S: Likewise.
+	* sysdeps/i386/fpu/s_cexpl.S: Likewise.
+	* math/libm-test.inc (cexp_test): Add more tests.
+	* sysdeps/i386/fpu/libm-test-ulps: Update.
+	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
+
 2012-03-21  Allan McRae  <allan@archlinux.org>
 
 	* timezone/Makefile: Do not install iso3166.tab and zone.tab
diff --git a/NEWS b/NEWS
index a8a3e57..5fc4444 100644
--- a/NEWS
+++ b/NEWS
@@ -16,7 +16,7 @@ Version 2.16
   13526, 13527, 13528, 13529, 13530, 13531, 13532, 13533, 13547, 13551,
   13552, 13553, 13555, 13559, 13566, 13583, 13618, 13637, 13656, 13658,
   13673, 13695, 13704, 13706, 13726, 13738, 13786, 13792, 13806, 13840,
-  13841, 13844, 13846, 13851, 13852, 13854, 13871
+  13841, 13844, 13846, 13851, 13852, 13854, 13871, 13883
 
 * ISO C11 support:
 
diff --git a/math/libm-test.inc b/math/libm-test.inc
index afc6f55..05a000e 100644
--- a/math/libm-test.inc
+++ b/math/libm-test.inc
@@ -1894,6 +1894,21 @@ cexp_test (void)
   TEST_c_c (cexp, 0.75L, 1.25L, 0.667537446429131586942201977015932112L, 2.00900045494094876258347228145863909L);
   TEST_c_c (cexp, -2.0, -3.0, -0.13398091492954261346140525546115575L, -0.019098516261135196432576240858800925L);
 
+  TEST_c_c (cexp, 0, 0x1p65, 0.99888622066058013610642172179340364209972L, -0.047183876212354673805106149805700013943218L);
+  TEST_c_c (cexp, 0, -0x1p65, 0.99888622066058013610642172179340364209972L, 0.047183876212354673805106149805700013943218L);
+  TEST_c_c (cexp, 50, 0x1p127, 4.053997150228616856622417636046265337193e21L, 3.232070315463388524466674772633810238819e21L);
+
+#ifndef TEST_FLOAT
+  TEST_c_c (cexp, 0, 1e22, 0.5232147853951389454975944733847094921409L, -0.8522008497671888017727058937530293682618L);
+  TEST_c_c (cexp, 0, 0x1p1023, -0.826369834614147994500785680811743734805L, 0.5631277798508840134529434079444683477104L);
+  TEST_c_c (cexp, 500, 0x1p1023, -1.159886268932754433233243794561351783426e217L, 7.904017694554466595359379965081774849708e216L);
+#endif
+
+#if defined TEST_LDOUBLE && LDBL_MAX_EXP >= 16384
+  TEST_c_c (cexp, 0, 0x1p16383L, 0.9210843909921906206874509522505756251609L, 0.3893629985894208126948115852610595405563L);
+  TEST_c_c (cexp, -10000, 0x1p16383L, 1.045876464564882298442774542991176546722e-4343L, 4.421154026488516836023811173959413420548e-4344L);
+#endif
+
   END (cexp, complex);
 }
 
diff --git a/sysdeps/i386/fpu/libm-test-ulps b/sysdeps/i386/fpu/libm-test-ulps
index 1d87514..4d61635 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 e5fdb7d..0000000
--- 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 6ed66e6..0000000
--- 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 ab02a17..0000000
--- 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 fef6ea1..e5eb8f9 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