This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix log1p missing underflows (bug 16339) [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Thu, 14 May 2015 23:39:02 +0000
- Subject: Fix log1p missing underflows (bug 16339) [committed]
- Authentication-results: sourceware.org; auth=none
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 --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;
}
--
Joseph S. Myers
joseph@codesourcery.com