mirror of
https://sourceware.org/git/glibc.git
synced 2024-11-23 13:30:06 +00:00
6ace393821
Similar to various other bugs in this area, pow functions can fail to raise the underflow exception when the result is tiny and inexact but one or more low bits of the intermediate result that is scaled down (or, in the i386 case, converted from a wider evaluation format) are zero. This patch forces the exception in a similar way to previous fixes, thereby concluding the fixes for known bugs with missing underflow exceptions currently filed in Bugzilla. Tested for x86_64, x86, mips64 and powerpc. [BZ #18825] * sysdeps/i386/fpu/i386-math-asm.h (FLT_NARROW_EVAL_UFLOW_NONNAN): New macro. (DBL_NARROW_EVAL_UFLOW_NONNAN): Likewise. (LDBL_CHECK_FORCE_UFLOW_NONNAN): Likewise. * sysdeps/i386/fpu/e_pow.S: Use DEFINE_DBL_MIN. (__ieee754_pow): Use DBL_NARROW_EVAL_UFLOW_NONNAN instead of DBL_NARROW_EVAL, reloading the PIC register as needed. * sysdeps/i386/fpu/e_powf.S: Use DEFINE_FLT_MIN. (__ieee754_powf): Use FLT_NARROW_EVAL_UFLOW_NONNAN instead of FLT_NARROW_EVAL. Use separate return path for case when first argument is NaN. * sysdeps/i386/fpu/e_powl.S: Include <i386-math-asm.h>. Use DEFINE_LDBL_MIN. (__ieee754_powl): Use LDBL_CHECK_FORCE_UFLOW_NONNAN, reloading the PIC register. * sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Use math_check_force_underflow_nonneg. * sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Force underflow for subnormal result. * sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Use math_check_force_underflow_nonneg. * sysdeps/x86/fpu/powl_helper.c (__powl_helper): Use math_check_force_underflow. * sysdeps/x86_64/fpu/x86_64-math-asm.h (LDBL_CHECK_FORCE_UFLOW_NONNAN): New macro. * sysdeps/x86_64/fpu/e_powl.S: Include <x86_64-math-asm.h>. Use DEFINE_LDBL_MIN. (__ieee754_powl): Use LDBL_CHECK_FORCE_UFLOW_NONNAN. * math/auto-libm-test-in: Add more tests of pow. * math/auto-libm-test-out: Regenerated.
415 lines
8.8 KiB
ArmAsm
415 lines
8.8 KiB
ArmAsm
/* ix87 specific implementation of pow function.
|
||
Copyright (C) 1996-2015 Free Software Foundation, Inc.
|
||
This file is part of the GNU C Library.
|
||
Contributed by Ulrich Drepper <drepper@cygnus.com>, 1996.
|
||
|
||
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 <machine/asm.h>
|
||
#include <x86_64-math-asm.h>
|
||
|
||
.section .rodata.cst8,"aM",@progbits,8
|
||
|
||
.p2align 3
|
||
.type one,@object
|
||
one: .double 1.0
|
||
ASM_SIZE_DIRECTIVE(one)
|
||
.type p3,@object
|
||
p3: .byte 0, 0, 0, 0, 0, 0, 0x20, 0x40
|
||
ASM_SIZE_DIRECTIVE(p3)
|
||
.type p63,@object
|
||
p63: .byte 0, 0, 0, 0, 0, 0, 0xe0, 0x43
|
||
ASM_SIZE_DIRECTIVE(p63)
|
||
.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
|
||
.type infinity,@object
|
||
inf_zero:
|
||
infinity:
|
||
.byte 0, 0, 0, 0, 0, 0, 0xf0, 0x7f
|
||
ASM_SIZE_DIRECTIVE(infinity)
|
||
.type zero,@object
|
||
zero: .double 0.0
|
||
ASM_SIZE_DIRECTIVE(zero)
|
||
.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)
|
||
#else
|
||
# define MO(op) op
|
||
#endif
|
||
|
||
.text
|
||
ENTRY(__ieee754_powl)
|
||
fldt 24(%rsp) // y
|
||
fxam
|
||
|
||
|
||
fnstsw
|
||
movb %ah, %dl
|
||
andb $0x45, %ah
|
||
cmpb $0x40, %ah // is y == 0 ?
|
||
je 11f
|
||
|
||
cmpb $0x05, %ah // is y == <EFBFBD>inf ?
|
||
je 12f
|
||
|
||
cmpb $0x01, %ah // is y == NaN ?
|
||
je 30f
|
||
|
||
fldt 8(%rsp) // x : y
|
||
|
||
fxam
|
||
fnstsw
|
||
movb %ah, %dh
|
||
andb $0x45, %ah
|
||
cmpb $0x40, %ah
|
||
je 20f // x is <EFBFBD>0
|
||
|
||
cmpb $0x05, %ah
|
||
je 15f // x is <EFBFBD>inf
|
||
|
||
cmpb $0x01, %ah
|
||
je 31f // x is NaN
|
||
|
||
fxch // y : x
|
||
|
||
/* fistpll raises invalid exception for |y| >= 1L<<63. */
|
||
fldl MO(p63) // 1L<<63 : y : x
|
||
fld %st(1) // y : 1L<<63 : y : x
|
||
fabs // |y| : 1L<<63 : y : x
|
||
fcomip %st(1), %st // 1L<<63 : y : x
|
||
fstp %st(0) // y : x
|
||
jnc 2f
|
||
|
||
/* First see whether `y' is a natural number. In this case we
|
||
can use a more precise algorithm. */
|
||
fld %st // y : y : x
|
||
fistpll -8(%rsp) // y : x
|
||
fildll -8(%rsp) // int(y) : y : x
|
||
fucomip %st(1),%st // y : x
|
||
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 < 8), use the algorithm for real exponent to avoid
|
||
accumulation of errors. */
|
||
fldl MO(p3) // 8 : y : x
|
||
fld %st(1) // y : 8 : y : x
|
||
fabs // |y| : 8 : y : x
|
||
fcomip %st(1), %st // 8 : y : x
|
||
fstp %st(0) // y : x
|
||
jnc 3f
|
||
mov -8(%rsp),%eax
|
||
mov -4(%rsp),%edx
|
||
orl $0, %edx
|
||
fstp %st(0) // x
|
||
jns 4f // y >= 0, jump
|
||
fdivrl MO(one) // 1/x (now referred to as x)
|
||
negl %eax
|
||
adcl $0, %edx
|
||
negl %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: 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 <20>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)
|
||
ret
|
||
|
||
.align ALIGNARG(4)
|
||
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. */
|
||
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,<EFBFBD>0) = 1
|
||
.align ALIGNARG(4)
|
||
11: fstp %st(0) // pop y
|
||
fldl MO(one)
|
||
ret
|
||
|
||
// y == <EFBFBD>inf
|
||
.align ALIGNARG(4)
|
||
12: fstp %st(0) // pop y
|
||
fldl MO(one) // 1
|
||
fldt 8(%rsp) // x : 1
|
||
fabs // abs(x) : 1
|
||
fucompp // < 1, == 1, or > 1
|
||
fnstsw
|
||
andb $0x45, %ah
|
||
cmpb $0x45, %ah
|
||
je 13f // jump if x is NaN
|
||
|
||
cmpb $0x40, %ah
|
||
je 14f // jump if |x| == 1
|
||
|
||
shlb $1, %ah
|
||
xorb %ah, %dl
|
||
andl $2, %edx
|
||
#ifdef PIC
|
||
lea inf_zero(%rip),%rcx
|
||
fldl (%rcx, %rdx, 4)
|
||
#else
|
||
fldl inf_zero(,%rdx, 4)
|
||
#endif
|
||
ret
|
||
|
||
.align ALIGNARG(4)
|
||
14: fldl MO(one)
|
||
ret
|
||
|
||
.align ALIGNARG(4)
|
||
13: fldt 8(%rsp) // load x == NaN
|
||
ret
|
||
|
||
.align ALIGNARG(4)
|
||
// x is <EFBFBD>inf
|
||
15: fstp %st(0) // y
|
||
testb $2, %dh
|
||
jz 16f // jump if x == +inf
|
||
|
||
// fistpll raises invalid exception for |y| >= 1L<<63, but y
|
||
// may be odd unless we know |y| >= 1L<<64.
|
||
fldl MO(p64) // 1L<<64 : y
|
||
fld %st(1) // y : 1L<<64 : y
|
||
fabs // |y| : 1L<<64 : y
|
||
fcomip %st(1), %st // 1L<<64 : y
|
||
fstp %st(0) // y
|
||
jnc 16f
|
||
fldl MO(p63) // p63 : y
|
||
fxch // y : p63
|
||
fprem // y%p63 : p63
|
||
fstp %st(1) // y%p63
|
||
|
||
// We must find out whether y is an odd integer.
|
||
fld %st // y : y
|
||
fistpll -8(%rsp) // y
|
||
fildll -8(%rsp) // int(y) : y
|
||
fucomip %st(1),%st
|
||
ffreep %st // <empty>
|
||
jne 17f
|
||
|
||
// OK, the value is an integer, but is it odd?
|
||
mov -8(%rsp), %eax
|
||
mov -4(%rsp), %edx
|
||
andb $1, %al
|
||
jz 18f // jump if not odd
|
||
// It's an odd integer.
|
||
shrl $31, %edx
|
||
#ifdef PIC
|
||
lea minf_mzero(%rip),%rcx
|
||
fldl (%rcx, %rdx, 8)
|
||
#else
|
||
fldl minf_mzero(,%rdx, 8)
|
||
#endif
|
||
ret
|
||
|
||
.align ALIGNARG(4)
|
||
16: fcompl MO(zero)
|
||
fnstsw
|
||
shrl $5, %eax
|
||
andl $8, %eax
|
||
#ifdef PIC
|
||
lea inf_zero(%rip),%rcx
|
||
fldl (%rcx, %rax, 1)
|
||
#else
|
||
fldl inf_zero(,%rax, 1)
|
||
#endif
|
||
ret
|
||
|
||
.align ALIGNARG(4)
|
||
17: shll $30, %edx // sign bit for y in right position
|
||
18: shrl $31, %edx
|
||
#ifdef PIC
|
||
lea inf_zero(%rip),%rcx
|
||
fldl (%rcx, %rdx, 8)
|
||
#else
|
||
fldl inf_zero(,%rdx, 8)
|
||
#endif
|
||
ret
|
||
|
||
.align ALIGNARG(4)
|
||
// x is <EFBFBD>0
|
||
20: fstp %st(0) // y
|
||
testb $2, %dl
|
||
jz 21f // y > 0
|
||
|
||
// x is <EFBFBD>0 and y is < 0. We must find out whether y is an odd integer.
|
||
testb $2, %dh
|
||
jz 25f
|
||
|
||
// fistpll raises invalid exception for |y| >= 1L<<63, but y
|
||
// may be odd unless we know |y| >= 1L<<64.
|
||
fldl MO(p64) // 1L<<64 : y
|
||
fld %st(1) // y : 1L<<64 : y
|
||
fabs // |y| : 1L<<64 : y
|
||
fcomip %st(1), %st // 1L<<64 : y
|
||
fstp %st(0) // y
|
||
jnc 25f
|
||
fldl MO(p63) // p63 : y
|
||
fxch // y : p63
|
||
fprem // y%p63 : p63
|
||
fstp %st(1) // y%p63
|
||
|
||
fld %st // y : y
|
||
fistpll -8(%rsp) // y
|
||
fildll -8(%rsp) // int(y) : y
|
||
fucomip %st(1),%st
|
||
ffreep %st // <empty>
|
||
jne 26f
|
||
|
||
// OK, the value is an integer, but is it odd?
|
||
mov -8(%rsp),%eax
|
||
mov -4(%rsp),%edx
|
||
andb $1, %al
|
||
jz 27f // jump if not odd
|
||
// It's an odd integer.
|
||
// Raise divide-by-zero exception and get minus infinity value.
|
||
fldl MO(one)
|
||
fdivl MO(zero)
|
||
fchs
|
||
ret
|
||
|
||
25: fstp %st(0)
|
||
26:
|
||
27: // Raise divide-by-zero exception and get infinity value.
|
||
fldl MO(one)
|
||
fdivl MO(zero)
|
||
ret
|
||
|
||
.align ALIGNARG(4)
|
||
// x is <EFBFBD>0 and y is > 0. We must find out whether y is an odd integer.
|
||
21: testb $2, %dh
|
||
jz 22f
|
||
|
||
// fistpll raises invalid exception for |y| >= 1L<<63, but y
|
||
// may be odd unless we know |y| >= 1L<<64.
|
||
fldl MO(p64) // 1L<<64 : y
|
||
fxch // y : 1L<<64
|
||
fcomi %st(1), %st // y : 1L<<64
|
||
fstp %st(1) // y
|
||
jnc 22f
|
||
fldl MO(p63) // p63 : y
|
||
fxch // y : p63
|
||
fprem // y%p63 : p63
|
||
fstp %st(1) // y%p63
|
||
|
||
fld %st // y : y
|
||
fistpll -8(%rsp) // y
|
||
fildll -8(%rsp) // int(y) : y
|
||
fucomip %st(1),%st
|
||
ffreep %st // <empty>
|
||
jne 23f
|
||
|
||
// OK, the value is an integer, but is it odd?
|
||
mov -8(%rsp),%eax
|
||
mov -4(%rsp),%edx
|
||
andb $1, %al
|
||
jz 24f // jump if not odd
|
||
// It's an odd integer.
|
||
fldl MO(mzero)
|
||
ret
|
||
|
||
22: fstp %st(0)
|
||
23:
|
||
24: fldl MO(zero)
|
||
ret
|
||
|
||
END(__ieee754_powl)
|
||
strong_alias (__ieee754_powl, __powl_finite)
|