diff mbox

Fix log1p missing underflows (bug 16339) [committed]

Message ID alpine.DEB.2.10.1505142338420.20495@digraph.polyomino.org.uk
State New
Headers show

Commit Message

Joseph Myers May 14, 2015, 11:39 p.m. UTC 
Similar to various other bugs in this area, some log1p implementations
do not raise the underflow exception for subnormal arguments, when the
result is tiny and inexact.  This patch forces the exception in a
similar way to previous fixes.  (The ldbl-128ibm implementation
doesn't currently need any change as it already generates this
exception, albeit through code that would generate spurious exceptions
in other cases; special code for this issue will only be needed there
when fixing the spurious exceptions.)

Tested for x86_64, x86, powerpc and mips64.  Committed.

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

2015-05-14  Joseph Myers  <joseph@codesourcery.com>

	[BZ #16339]
	* sysdeps/i386/fpu/s_log1p.S (dbl_min): New object.
	(__log1p): Force underflow exception for results with small
	absolute value.
	* sysdeps/i386/fpu/s_log1pf.S (flt_min): New object.
	(__log1pf): Force underflow exception for results with small
	absolute value.
	* sysdeps/ieee754/dbl-64/s_log1p.c: Include <float.h>.
	(__log1p): Force underflow exception for results with small
	absolute value.
	* sysdeps/ieee754/flt-32/s_log1pf.c: Include <float.h>.
	(__log1pf): Force underflow exception for results with small
	absolute value.
	* sysdeps/ieee754/ldbl-128/s_log1pl.c: Include <float.h>.
	(__log1pl): Force underflow exception for results with small
	absolute value.
	* math/auto-libm-test-in: Do not allow missing underflow
	exceptions from log1p.
	* math/auto-libm-test-out: Regenerated.
diff mbox

Patch

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 883397c..7e5c7bb 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1850,11 +1850,10 @@  log1p -0
 log1p e-1
 log1p -0.25
 log1p -0.875
-# Bug 16339: underflow exception may be missing.
-log1p min missing-underflow
-log1p min_subnorm missing-underflow
-log1p -min missing-underflow
-log1p -min_subnorm missing-underflow
+log1p min
+log1p min_subnorm
+log1p -min
+log1p -min_subnorm
 log1p 0x1p10
 log1p 0x1p20
 log1p 0x1p30
diff --git a/sysdeps/i386/fpu/s_log1p.S b/sysdeps/i386/fpu/s_log1p.S
index 8624249..c2559a3 100644
--- a/sysdeps/i386/fpu/s_log1p.S
+++ b/sysdeps/i386/fpu/s_log1p.S
@@ -17,6 +17,13 @@  RCSID("$NetBSD: s_log1p.S,v 1.7 1995/05/09 00:10:58 jtc Exp $")
 limit:	.double 0.29
 one:	.double 1.0
 
+	.section .rodata.cst8,"aM",@progbits,8
+
+	.p2align 3
+	.type dbl_min,@object
+dbl_min:	.byte 0, 0, 0, 0, 0, 0, 0x10, 0
+	ASM_SIZE_DIRECTIVE(dbl_min)
+
 /*
  * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29,
  * otherwise fyl2x with the needed extra computation.
@@ -55,7 +62,25 @@  ENTRY(__log1p)
 	ret
 
 2:	fyl2xp1
-	ret
+#ifdef PIC
+	fldl	dbl_min@GOTOFF(%edx)
+#else
+	fldl	dbl_min
+#endif
+	fld	%st(1)
+	fabs
+	fucompp
+	fnstsw
+	sahf
+	jnc	1f
+	subl	$8, %esp
+	cfi_adjust_cfa_offset (8)
+	fld	%st(0)
+	fmul	%st(0)
+	fstpl	(%esp)
+	addl	$8, %esp
+	cfi_adjust_cfa_offset (-8)
+1:	ret
 
 3:	jp	4b		// in case x is ?Inf
 	fstp	%st(1)
diff --git a/sysdeps/i386/fpu/s_log1pf.S b/sysdeps/i386/fpu/s_log1pf.S
index b071e73..8fca22e 100644
--- a/sysdeps/i386/fpu/s_log1pf.S
+++ b/sysdeps/i386/fpu/s_log1pf.S
@@ -17,6 +17,13 @@  RCSID("$NetBSD: s_log1pf.S,v 1.4 1995/05/09 00:13:05 jtc Exp $")
 limit:	.float 0.29
 one:	.float 1.0
 
+	.section .rodata.cst4,"aM",@progbits,4
+
+	.p2align 2
+	.type flt_min,@object
+flt_min:	.byte 0, 0, 0x80, 0
+	ASM_SIZE_DIRECTIVE(flt_min)
+
 /*
  * Use the fyl2xp1 function when the argument is in the range -0.29 to 0.29,
  * otherwise fyl2x with the needed extra computation.
@@ -55,7 +62,25 @@  ENTRY(__log1pf)
 	ret
 
 2:	fyl2xp1
-	ret
+#ifdef PIC
+	flds	flt_min@GOTOFF(%edx)
+#else
+	flds	flt_min
+#endif
+	fld	%st(1)
+	fabs
+	fucompp
+	fnstsw
+	sahf
+	jnc	1f
+	subl	$4, %esp
+	cfi_adjust_cfa_offset (4)
+	fld	%st(0)
+	fmul	%st(0)
+	fstps	(%esp)
+	addl	$4, %esp
+	cfi_adjust_cfa_offset (-4)
+1:	ret
 
 3:	jp	4b		// in case x is ?Inf
 	fstp	%st(1)
diff --git a/sysdeps/ieee754/dbl-64/s_log1p.c b/sysdeps/ieee754/dbl-64/s_log1p.c
index 86bbfba..cff555b 100644
--- a/sysdeps/ieee754/dbl-64/s_log1p.c
+++ b/sysdeps/ieee754/dbl-64/s_log1p.c
@@ -78,6 +78,7 @@ 
  *	 See HP-15C Advanced Functions Handbook, p.193.
  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -118,7 +119,14 @@  __log1p (double x)
 	{
 	  math_force_eval (two54 + x);          /* raise inexact */
 	  if (ax < 0x3c900000)                  /* |x| < 2**-54 */
-	    return x;
+	    {
+	      if (fabs (x) < DBL_MIN)
+		{
+		  double force_underflow = x * x;
+		  math_force_eval (force_underflow);
+		}
+	      return x;
+	    }
 	  else
 	    return x - x * x * 0.5;
 	}
diff --git a/sysdeps/ieee754/flt-32/s_log1pf.c b/sysdeps/ieee754/flt-32/s_log1pf.c
index 94c33fc..83a09f1 100644
--- a/sysdeps/ieee754/flt-32/s_log1pf.c
+++ b/sysdeps/ieee754/flt-32/s_log1pf.c
@@ -13,6 +13,7 @@ 
  * ====================================================
  */
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -48,7 +49,14 @@  __log1pf(float x)
 	    if(ax<0x31000000) {			/* |x| < 2**-29 */
 		math_force_eval(two25+x);	/* raise inexact */
 		if (ax<0x24800000)		/* |x| < 2**-54 */
+		  {
+		    if (fabsf (x) < FLT_MIN)
+		      {
+			float force_underflow = x * x;
+			math_force_eval (force_underflow);
+		      }
 		    return x;
+		  }
 		else
 		    return x - x*x*(float)0.5;
 	    }
diff --git a/sysdeps/ieee754/ldbl-128/s_log1pl.c b/sysdeps/ieee754/ldbl-128/s_log1pl.c
index b70a55b..ff759bc 100644
--- a/sysdeps/ieee754/ldbl-128/s_log1pl.c
+++ b/sysdeps/ieee754/ldbl-128/s_log1pl.c
@@ -53,6 +53,7 @@ 
     <http://www.gnu.org/licenses/>.  */
 
 
+#include <float.h>
 #include <math.h>
 #include <math_private.h>
 
@@ -140,6 +141,11 @@  __log1pl (long double xm1)
 
   if ((hx & 0x7fffffff) < 0x3f8e0000)
     {
+      if (fabsl (xm1) < LDBL_MIN)
+	{
+	  long double force_underflow = xm1 * xm1;
+	  math_force_eval (force_underflow);
+	}
       if ((int) xm1 == 0)
 	return xm1;
     }