Skip to content

Commit

Permalink
Fix log1p missing underflows (bug 16339).
Browse files Browse the repository at this point in the history 
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.

	[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.
  • Loading branch information
Joseph Myers committed May 14, 2015
1 parent 95b07fb commit 0b7a5f9
Show file tree
Hide file tree
Showing 9 changed files with 202 additions and 109 deletions.
22 changes: 22 additions & 0 deletions ChangeLog
View file
Original file line number Diff line number Diff line change
@@ -1,3 +1,25 @@
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.

2015-05-14 Jakub Bogusz <qboosh@pld-linux.org>
Adhemerval Zanella <adhemerval.zanella@linaro.org>

Expand Down
2 changes: 1 addition & 1 deletion NEWS
View file
Original file line number Diff line number Diff line change
Expand Up @@ -9,7 +9,7 @@ Version 2.22

* The following bugs are resolved with this release:

4719, 6792, 13064, 14094, 14841, 14906, 15319, 15467, 15790, 15969,
4719, 6792, 13064, 14094, 14841, 14906, 15319, 15467, 15790, 15969, 16339,
16351, 16512, 16560, 16783, 16850, 17090, 17195, 17269, 17523, 17542,
17569, 17588, 17596, 17620, 17621, 17628, 17631, 17692, 17711, 17715,
17776, 17779, 17792, 17836, 17912, 17916, 17930, 17932, 17944, 17949,
Expand Down
9 changes: 4 additions & 5 deletions math/auto-libm-test-in
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down
200 changes: 100 additions & 100 deletions math/auto-libm-test-out
View file

Large diffs are not rendered by default.

27 changes: 26 additions & 1 deletion sysdeps/i386/fpu/s_log1p.S
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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)
Expand Down
27 changes: 26 additions & 1 deletion sysdeps/i386/fpu/s_log1pf.S
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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.
Expand Down Expand Up @@ -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)
Expand Down
10 changes: 9 additions & 1 deletion sysdeps/ieee754/dbl-64/s_log1p.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -78,6 +78,7 @@
* See HP-15C Advanced Functions Handbook, p.193.
*/

#include <float.h>
#include <math.h>
#include <math_private.h>

Expand Down Expand Up @@ -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;
}
Expand Down
8 changes: 8 additions & 0 deletions sysdeps/ieee754/flt-32/s_log1pf.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -13,6 +13,7 @@
* ====================================================
*/

#include <float.h>
#include <math.h>
#include <math_private.h>

Expand Down Expand Up @@ -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;
}
Expand Down
6 changes: 6 additions & 0 deletions sysdeps/ieee754/ldbl-128/s_log1pl.c
View file
Original file line number Diff line number Diff line change
Expand Up @@ -53,6 +53,7 @@
<http://www.gnu.org/licenses/>. */


#include <float.h>
#include <math.h>
#include <math_private.h>

Expand Down Expand Up @@ -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;
}
Expand Down

0 comments on commit 0b7a5f9

Please sign in to comment.