From 1a4ac776ebc9bb07287f59f63e473db531318dff Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Wed, 21 Mar 2012 15:28:05 +0000 Subject: Remove inaccurate x86 cexp implementations (bug 13883). --- sysdeps/i386/fpu/libm-test-ulps | 86 ++++++++++---- sysdeps/i386/fpu/s_cexp.S | 253 --------------------------------------- sysdeps/i386/fpu/s_cexpf.S | 257 ---------------------------------------- sysdeps/i386/fpu/s_cexpl.S | 256 --------------------------------------- 4 files changed, 66 insertions(+), 786 deletions(-) delete mode 100644 sysdeps/i386/fpu/s_cexp.S delete mode 100644 sysdeps/i386/fpu/s_cexpf.S delete mode 100644 sysdeps/i386/fpu/s_cexpl.S (limited to 'sysdeps/i386') 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 , 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 - . */ - -#include - - .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) /* */ - 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 , 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 - . */ - -#include - - .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) /* */ - 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 , 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 - . */ - -#include - - .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) /* */ - 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) -- cgit 1.4.1