/* ix87 specific implementation of pow function.
- Copyright (C) 1996-1999, 2001, 2004, 2007, 2011-2012
- Free Software Foundation, Inc.
+ Copyright (C) 1996-2019 Free Software Foundation, Inc.
This file is part of the GNU C Library.
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
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/>. */
+ <https://www.gnu.org/licenses/>. */
#include <machine/asm.h>
+#include <x86_64-math-asm.h>
.section .rodata.cst8,"aM",@progbits,8
.p2align 3
- ASM_TYPE_DIRECTIVE(one,@object)
+ .type one,@object
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
- ASM_TYPE_DIRECTIVE(limit,@object)
-limit: .double 0.29
- ASM_SIZE_DIRECTIVE(limit)
- ASM_TYPE_DIRECTIVE(p63,@object)
+ .type p2,@object
+p2: .byte 0, 0, 0, 0, 0, 0, 0x10, 0x40
+ ASM_SIZE_DIRECTIVE(p2)
+ .type p63,@object
p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
ASM_SIZE_DIRECTIVE(p63)
- ASM_TYPE_DIRECTIVE(p64,@object)
+ .type p64,@object
p64: .byte 0, 0, 0, 0, 0, 0, 0xf0, 0x43
ASM_SIZE_DIRECTIVE(p64)
+ .type p78,@object
+p78: .byte 0, 0, 0, 0, 0, 0, 0xd0, 0x44
+ ASM_SIZE_DIRECTIVE(p78)
+ .type pm79,@object
+pm79: .byte 0, 0, 0, 0, 0, 0, 0, 0x3b
+ ASM_SIZE_DIRECTIVE(pm79)
.section .rodata.cst16,"aM",@progbits,16
.p2align 3
- ASM_TYPE_DIRECTIVE(infinity,@object)
+ .type infinity,@object
inf_zero:
infinity:
.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
ASM_SIZE_DIRECTIVE(infinity)
- ASM_TYPE_DIRECTIVE(zero,@object)
+ .type zero,@object
zero: .double 0.0
ASM_SIZE_DIRECTIVE(zero)
- ASM_TYPE_DIRECTIVE(minf_mzero,@object)
+ .type minf_mzero,@object
minf_mzero:
minfinity:
.byte 0, 0, 0, 0, 0, 0, 0xf0, 0xff
mzero:
.byte 0, 0, 0, 0, 0, 0, 0, 0x80
ASM_SIZE_DIRECTIVE(minf_mzero)
+DEFINE_LDBL_MIN
#ifdef PIC
# define MO(op) op##(%rip)
cmpb $0x05, %ah
je 15f // x is ±inf
+ cmpb $0x01, %ah
+ je 31f // x is NaN
+
fxch // y : x
/* fistpll raises invalid exception for |y| >= 1L<<63. */
fistpll -8(%rsp) // y : x
fildll -8(%rsp) // int(y) : y : x
fucomip %st(1),%st // y : x
- jne 3f
-
- /* OK, we have an integer value for y. */
+ je 9f
+
+ // If y has absolute value at most 0x1p-79, then any finite
+ // nonzero x will result in 1. Saturate y to those bounds to
+ // avoid underflow in the calculation of y*log2(x).
+ fldl MO(pm79) // 0x1p-79 : y : x
+ fld %st(1) // y : 0x1p-79 : y : x
+ fabs // |y| : 0x1p-79 : y : x
+ fcomip %st(1), %st // 0x1p-79 : y : x
+ fstp %st(0) // y : x
+ jnc 3f
+ fstp %st(0) // pop y
+ fldl MO(pm79) // 0x1p-79 : x
+ testb $2, %dl
+ jnz 3f // y > 0
+ fchs // -0x1p-79 : x
+ jmp 3f
+
+9: /* OK, we have an integer value for y. Unless very small
+ (we use < 4), use the algorithm for real exponent to avoid
+ accumulation of errors. */
+ fldl MO(p2) // 4 : y : x
+ fld %st(1) // y : 4 : y : x
+ fabs // |y| : 4 : y : x
+ fcomip %st(1), %st // 4 : y : x
+ fstp %st(0) // y : x
+ jnc 3f
mov -8(%rsp),%eax
mov -4(%rsp),%edx
orl $0, %edx
4: fldl MO(one) // 1 : x
fxch
+ /* If y is even, take the absolute value of x. Otherwise,
+ ensure all intermediate values that might overflow have the
+ sign of x. */
+ testb $1, %al
+ jnz 6f
+ fabs
+
6: shrdl $1, %edx, %eax
jnc 5f
fxch
+ fabs
fmul %st(1) // x : ST*x
fxch
-5: fmul %st(0), %st // x*x : ST*x
+5: fld %st // x : x : ST*x
+ fabs // |x| : x : ST*x
+ fmulp // |x|*x : ST*x
shrl $1, %edx
movl %eax, %ecx
orl %edx, %ecx
jnz 6b
fstp %st(0) // ST*x
+ LDBL_CHECK_FORCE_UFLOW_NONNAN
ret
/* y is ±NAN */
30: fldt 8(%rsp) // x : y
fldl MO(one) // 1.0 : x : y
fucomip %st(1),%st // x : y
- je 31f
- fxch // y : x
-31: fstp %st(1)
+ je 32f
+31: /* At least one argument NaN, and result should be NaN. */
+ faddp
+ ret
+32: jc 31b
+ /* pow (1, NaN); check if the NaN signaling. */
+ testb $0x40, 31(%rsp)
+ jz 31b
+ fstp %st(1)
ret
.align ALIGNARG(4)
-2: // y is a large integer (absolute value at least 1L<<63), but
- // may be odd unless at least 1L<<64. So it may be necessary
- // to adjust the sign of a negative result afterwards.
- fxch // x : y
- fabs // |x| : y
- fxch // y : |x|
+2: // y is a large integer (absolute value at least 1L<<63).
+ // If y has absolute value at least 1L<<78, then any finite
+ // nonzero x will result in 0 (underflow), 1 or infinity (overflow).
+ // Saturate y to those bounds to avoid overflow in the calculation
+ // of y*log2(x).
+ fldl MO(p78) // 1L<<78 : y : x
+ fld %st(1) // y : 1L<<78 : y : x
+ fabs // |y| : 1L<<78 : y : x
+ fcomip %st(1), %st // 1L<<78 : y : x
+ fstp %st(0) // y : x
+ jc 3f
+ fstp %st(0) // pop y
+ fldl MO(p78) // 1L<<78 : x
+ testb $2, %dl
+ jz 3f // y > 0
+ fchs // -(1L<<78) : x
.align ALIGNARG(4)
3: /* y is a real number. */
- fxch // x : y
- fldl MO(one) // 1.0 : x : y
- fldl MO(limit) // 0.29 : 1.0 : x : y
- fld %st(2) // x : 0.29 : 1.0 : x : y
- fsub %st(2) // x-1 : 0.29 : 1.0 : x : y
- fabs // |x-1| : 0.29 : 1.0 : x : y
- fucompp // 1.0 : x : y
- fnstsw
- fxch // x : 1.0 : y
- test $0x4500,%eax
- jz 7f
- fsub %st(1) // x-1 : 1.0 : y
- fyl2xp1 // log2(x) : y
- jmp 8f
-
-7: fyl2x // log2(x) : y
-8: fmul %st(1) // y*log2(x) : y
- fxam
- fnstsw
- andb $0x45, %ah
- cmpb $0x05, %ah // is y*log2(x) == ±inf ?
- je 28f
- fst %st(1) // y*log2(x) : y*log2(x)
- frndint // int(y*log2(x)) : y*log2(x)
- fsubr %st, %st(1) // int(y*log2(x)) : fract(y*log2(x))
- fxch // fract(y*log2(x)) : int(y*log2(x))
- f2xm1 // 2^fract(y*log2(x))-1 : int(y*log2(x))
- faddl MO(one) // 2^fract(y*log2(x)) : int(y*log2(x))
- fscale // 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
- fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
- jmp 29f
-
-28: fstp %st(1) // y*log2(x)
- fldl MO(one) // 1 : y*log2(x)
- fscale // 2^(y*log2(x)) : y*log2(x)
- fstp %st(1) // 2^(y*log2(x))
-29: testb $2, %dh
- jz 292f
- // x is negative. If y is an odd integer, negate the result.
- fldt 24(%rsp) // y : abs(result)
- fldl MO(p64) // 1L<<64 : y : abs(result)
- fld %st(1) // y : 1L<<64 : y : abs(result)
- fabs // |y| : 1L<<64 : y : abs(result)
- fcomip %st(1), %st // 1L<<64 : y : abs(result)
- fstp %st(0) // y : abs(result)
- jnc 291f
- fldl MO(p63) // p63 : y : abs(result)
- fxch // y : p63 : abs(result)
- fprem // y%p63 : p63 : abs(result)
- fstp %st(1) // y%p63 : abs(result)
-
- // We must find out whether y is an odd integer.
- fld %st // y : y : abs(result)
- fistpll -8(%rsp) // y : abs(result)
- fildll -8(%rsp) // int(y) : y : abs(result)
- fucomip %st(1),%st // y : abs(result)
- ffreep %st // abs(result)
- jne 292f
-
- // OK, the value is an integer, but is it odd?
- mov -8(%rsp), %eax
- mov -4(%rsp), %edx
- andb $1, %al
- jz 290f // jump if not odd
- // It's an odd integer.
- fchs
-290: ret
-291: fstp %st(0) // abs(result)
-292: ret
+ subq $40, %rsp
+ cfi_adjust_cfa_offset (40)
+ fstpt 16(%rsp) // x
+ fstpt (%rsp) // <empty>
+ call HIDDEN_JUMPTARGET (__powl_helper) // <result>
+ addq $40, %rsp
+ cfi_adjust_cfa_offset (-40)
+ ret
- // pow(x,±0) = 1
+ // pow(x,±0) = 1, unless x is sNaN
.align ALIGNARG(4)
11: fstp %st(0) // pop y
+ fldt 8(%rsp) // x
+ fxam
+ fnstsw
+ andb $0x45, %ah
+ cmpb $0x01, %ah
+ je 112f // x is NaN
+111: fstp %st(0)
fldl MO(one)
ret
+112: testb $0x40, 15(%rsp)
+ jnz 111b
+ fadd %st(0)
+ ret
+
// y == ±inf
.align ALIGNARG(4)
12: fstp %st(0) // pop y
.align ALIGNARG(4)
13: fldt 8(%rsp) // load x == NaN
+ fadd %st(0)
ret
.align ALIGNARG(4)
projects / thirdparty / glibc.git / blobdiff