Fix log1pl (LDBL_MAX) in FE_UPWARD mode (bug 16564)
diff mbox

Message ID Pine.LNX.4.64.1405122017540.25715@digraph.polyomino.org.uk
State New
Headers show

Commit Message

Joseph Myers May 12, 2014, 8:18 p.m. UTC
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.

(auto-libm-test-out diffs omitted below.)

2014-05-12  Joseph Myers  <joseph@codesourcery.com>

	[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.

Comments

Andreas Jaeger May 14, 2014, 10:24 a.m. UTC | #1
On 05/12/2014 04:18 PM, Joseph S. Myers wrote:
> 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).

Looks fine, thanks,
Andreas

Patch
diff mbox

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 420f657..85f3c4f 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1417,6 +1417,14 @@  log1p min missing-underflow
 log1p min_subnorm missing-underflow
 log1p -min missing-underflow
 log1p -min_subnorm missing-underflow
+log1p 0x1p10
+log1p 0x1p20
+log1p 0x1p30
+log1p 0x1p50
+log1p 0x1p60
+log1p 0x1p100
+log1p 0x1p1000
+log1p max
 
 log2 1
 log2 e
diff --git a/sysdeps/i386/fpu/s_log1pl.S b/sysdeps/i386/fpu/s_log1pl.S
index 93c07b6..d2d5d3b 100644
--- a/sysdeps/i386/fpu/s_log1pl.S
+++ b/sysdeps/i386/fpu/s_log1pl.S
@@ -53,12 +53,17 @@  ENTRY(__log1pl)
 	sahf
 	jnc	2f
 
+	movzwl	4+8(%esp), %eax
+	xorb	$0x80, %ah
+	cmpl	$0xc040, %eax
+	jae	5f
+
 #ifdef PIC
 	faddl	one@GOTOFF(%edx)
 #else
 	faddl	one
 #endif
-	fyl2x
+5:	fyl2x
 	ret
 
 2:	fyl2xp1
diff --git a/sysdeps/ieee754/ldbl-128/s_log1pl.c b/sysdeps/ieee754/ldbl-128/s_log1pl.c
index d8d89f0..4a30af6 100644
--- a/sysdeps/ieee754/ldbl-128/s_log1pl.c
+++ b/sysdeps/ieee754/ldbl-128/s_log1pl.c
@@ -144,7 +144,10 @@  __log1pl (long double xm1)
 	return xm1;
     }
 
-  x = xm1 + 1.0L;
+  if (xm1 >= 0x1p113L)
+    x = xm1;
+  else
+    x = xm1 + 1.0L;
 
   /* log1p(-1) = -inf */
   if (x <= 0.0L)
diff --git a/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c b/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c
index a346383..e4bb6e8 100644
--- a/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/s_log1pl.c
@@ -140,7 +140,10 @@  __log1pl (long double xm1)
   if (((hx & 0x7fffffff) | lx) == 0)
     return xm1;
 
-  x = xm1 + 1.0L;
+  if (xm1 >= 0x1p107L)
+    x = xm1;
+  else
+    x = xm1 + 1.0L;
 
   /* log1p(-1) = -inf */
   if (x <= 0.0L)
diff --git a/sysdeps/x86_64/fpu/s_log1pl.S b/sysdeps/x86_64/fpu/s_log1pl.S
index b4dbcdf..af3024a 100644
--- a/sysdeps/x86_64/fpu/s_log1pl.S
+++ b/sysdeps/x86_64/fpu/s_log1pl.S
@@ -52,8 +52,13 @@  ENTRY(__log1pl)
 	andb	$1,%ah
 	jz	2f
 
+	movzwl	8+8(%rsp), %eax
+	xorb	$0x80, %ah
+	cmpl	$0xc040, %eax
+	jae	5f
+
 	faddl	MO(one)
-	fyl2x
+5:	fyl2x
 	ret
 
 2:	fyl2xp1