glibc/sysdeps/x86_64/fpu/e_log10l.S
Joseph Myers b72592e75f Fix log10 (1) in round-downward mode (bug 16977).
As with various other issues of this kind, bug 16977 is log10 (1)
wrongly returning -0 rather than +0 in round-downward mode because of
an implementation effectively in terms of log1p (x - 1).  This patch
fixes the issue in the same way used for log.

Tested x86_64 and x86 and ulps updated accordingly.  Also tested for
mips64 to confirm a fix was needed for ldbl-128 and to validate that
fix (also applied to ldbl-128ibm since that version of logl is
essentially the same as the ldbl-128 one).

	[BZ #16977]
	* sysdeps/i386/fpu/e_log10.S (__ieee754_log10): Take absolute
	value when x - 1 is zero.
	* sysdeps/i386/fpu/e_log10f.S (__ieee754_log10f): Likewise.
	* sysdeps/i386/fpu/e_log10l.S (__ieee754_log10l): Likewise.
	* sysdeps/ieee754/ldbl-128/e_log10l.c (__ieee754_log10l): Return
	0.0L for an argument of 1.0L.
	* sysdeps/ieee754/ldbl-128ibm/e_log10l.c (__ieee754_log10l):
	Likewise.
	* sysdeps/x86_64/fpu/e_log10l.S (__ieee754_log10l): Take absolute
	value when x - 1 is zero.
	* math/libm-test.inc (log10_test): Use ALL_RM_TEST.
	* sysdeps/i386/fpu/libm-test-ulps: Update.
	* sysdeps/x86_64/fpu/libm-test-ulps: Likewise.
2014-05-23 12:07:50 +00:00

86 lines
1.8 KiB
ArmAsm

/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
*
* Changed to use fyl2xp1 for values near 1, <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
*/
#include <machine/asm.h>
.section .rodata.cst8,"aM",@progbits,8
.p2align 3
.type one,@object
one: .double 1.0
ASM_SIZE_DIRECTIVE(one)
/* It is not important that this constant is precise. It is only
a value which is known to be on the safe side for using the
fyl2xp1 instruction. */
.type limit,@object
limit: .double 0.29
ASM_SIZE_DIRECTIVE(limit)
#ifdef PIC
# define MO(op) op##(%rip)
#else
# define MO(op) op
#endif
.text
ENTRY(__ieee754_log10l)
fldlg2 // log10(2)
fldt 8(%rsp) // x : log10(2)
fxam
fnstsw
fld %st // x : x : log10(2)
testb $1, %ah
jnz 3f // in case x is NaN or ±Inf
4: fsubl MO(one) // x-1 : x : log10(2)
fld %st // x-1 : x-1 : x : log10(2)
fabs // |x-1| : x-1 : x : log10(2)
fcompl MO(limit) // x-1 : x : log10(2)
fnstsw // x-1 : x : log10(2)
andb $0x45, %ah
jz 2f
fxam
fnstsw
andb $0x45, %ah
cmpb $0x40, %ah
jne 5f
fabs // log10(1) is +0 in all rounding modes.
5: fstp %st(1) // x-1 : log10(2)
fyl2xp1 // log10(x)
ret
2: fstp %st(0) // x : log10(2)
fyl2x // log10(x)
ret
3: testb $4, %ah
jnz 4b // in case x is ±Inf
fstp %st(1)
fstp %st(1)
ret
END(__ieee754_log10l)
ENTRY(__log10l_finite)
fldlg2 // log10(2)
fldt 8(%rsp) // x : log10(2)
fld %st // x : x : log10(2)
4: fsubl MO(one) // x-1 : x : log10(2)
fld %st // x-1 : x-1 : x : log10(2)
fabs // |x-1| : x-1 : x : log10(2)
fcompl MO(limit) // x-1 : x : log10(2)
fnstsw // x-1 : x : log10(2)
andb $0x45, %ah
jz 2b
fstp %st(1) // x-1 : log10(2)
fyl2xp1 // log10(x)
ret
END(__log10l_finite)