Fix log1pl (LDBL_MAX) in FE_UPWARD mode (bug 16564)
Commit Message
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
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
@@ -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
@@ -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
@@ -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)
@@ -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)
@@ -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