Remove inaccurate x86 cexp implementations (bug 13883).

This commit is contained in:
Joseph Myers 2012-03-21 15:28:05 +00:00
parent a458e7fe38
commit 1a4ac776eb
8 changed files with 114 additions and 787 deletions

View file

@ -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

2
NEWS
View file

@ -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:

View file

@ -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);
}

View file

@ -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

View file

@ -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)

View file

@ -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)

View file

@ -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)

View file

@ -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