glibc/sysdeps/x86_64/fpu/s_log1pl.S
Joseph Myers 1a84c3d6d4 Fix log1pl (LDBL_MAX) in FE_UPWARD mode (bug 16564).
Bug 16564 is spurious overflow of log1pl (LDBL_MAX) in FE_UPWARD mode,
resulting from log1pl adding 1 to its argument (for arguments not
close to 0), which overflows in that mode.  This patch fixes this by
avoiding adding 1 to large arguments (precisely what counts as large
depends on the floating-point format).

Tested x86_64 and x86, and spot-checked log1pl tests on mips64 and
powerpc64.

	[BZ #16564]
	* sysdeps/i386/fpu/s_log1pl.S (__log1pl): Do not add 1 to positive
	arguments with exponent 65 or above.
	* sysdeps/ieee754/ldbl-128/s_log1pl.c (__log1pl): Do not add 1 to
	arguments 0x1p113L or above.
	* sysdeps/ieee754/ldbl-128ibm/s_log1pl.c (__log1pl): Do not add 1
	to arguments 0x1p107L or above.
	* sysdeps/x86_64/fpu/s_log1pl.S (__log1pl): Do not add 1 to
	positive arguments with exponent 65 or above.
	* math/auto-libm-test-in: Add more tests of log1p.
	* math/auto-libm-test-out: Regenerated.
2014-05-14 12:38:56 +00:00

75 lines
1.3 KiB
ArmAsm

/*
* Written by J.T. Conklin <jtc@netbsd.org>.
* Public domain.
*
* Adapted for `long double' by Ulrich Drepper <drepper@cygnus.com>.
* Adapted for x86-64 by Andreas Jaeger <aj@suse.de>.
*/
#include <machine/asm.h>
RCSID("$NetBSD: s_log1p.S,v 1.7 1995/05/09 00:10:58 jtc Exp $")
.section .rodata
.align ALIGNARG(4)
/* The fyl2xp1 can only be used for values in
-1 + sqrt(2) / 2 <= x <= 1 - sqrt(2) / 2
0.29 is a safe value.
*/
limit: .tfloat 0.29
/* Please note: we use a double value here. Since 1.0 has
an exact representation this does not effect the accuracy
but it helps to optimize the code. */
one: .double 1.0
/*
* Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29,
* otherwise fyl2x with the needed extra computation.
*/
#ifdef PIC
#define MO(op) op##(%rip)
#else
#define MO(op) op
#endif
.text
ENTRY(__log1pl)
fldln2
fldt 8(%rsp)
fxam
fnstsw
fld %st
testb $1, %ah
jnz 3f // in case x is NaN or ±Inf
4:
fabs
fldt MO(limit)
fcompp
fnstsw
andb $1,%ah
jz 2f
movzwl 8+8(%rsp), %eax
xorb $0x80, %ah
cmpl $0xc040, %eax
jae 5f
faddl MO(one)
5: fyl2x
ret
2: fyl2xp1
ret
3: testb $4, %ah
jnz 4b // in case x is ±Inf
fstp %st(1)
fstp %st(1)
ret
END (__log1pl)
weak_alias (__log1pl, log1pl)