This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Remove inaccurate x86 cexp implementations (bug 13883)
- From: "Joseph S. Myers" <joseph at codesourcery dot com>
- To: libc-alpha at sourceware dot org
- Date: Wed, 21 Mar 2012 14:52:04 +0000 (UTC)
- Subject: Remove inaccurate x86 cexp implementations (bug 13883)
This patch fixes bug 13883, cexp inaccurate for large imaginary parts
on x86, the same way other functions using fsincos and related
instructions were fixed: remove the x86 versions and use the generic C
versions instead, which have similar accuracy for small inputs and
much greater accuracy as the range reduction of the x86 versions gets
worse. The testcases added are based on those for sincos.
Tested x86 and x86_64 and ulps updated accordingly.
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.
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
--
Joseph S. Myers
joseph@codesourcery.com