1
0
mirror of https://sourceware.org/git/glibc.git synced 2025-10-24 13:33:08 +03:00

Fix exp2 spurious underflows (bug 16560).

This patch fixes the remaining part of bug 16560, spurious underflows
from exp2 of arguments close to 0 (when the result is close to 1, so
should not underflow), by just using 1+x instead of a more complicated
calculation when the argument is sufficiently small.

Tested for x86_64, x86 and mips64.

	[BZ #16560]
	* math/e_exp2l.c [LDBL_MANT_DIG == 106] (LDBL_EPSILON): Undefine
	and redefine.
	(__ieee754_exp2l): Do not multiply small fractional parts by
	M_LN2l.
	* sysdeps/i386/fpu/e_exp2l.S (__ieee754_exp2l): Just add 1 to
	small argument.
	* sysdeps/ieee754/dbl-64/e_exp2.c (__ieee754_exp2): Likewise.
	* sysdeps/ieee754/flt-32/e_exp2f.c (__ieee754_exp2f): Likewise.
	* sysdeps/x86_64/fpu/e_exp2l.S (__ieee754_exp2l): Likewise.
	* math/auto-libm-test-in: Add more tests of exp2.
	* math/auto-libm-test-out: Regenerated.
This commit is contained in:
Joseph Myers
2015-02-12 19:02:45 +00:00
parent 98408b95b1
commit 03d95bd483
9 changed files with 733 additions and 3 deletions

View File

@@ -18,7 +18,15 @@ ENTRY(__ieee754_exp2l)
andb %ah, %dh
cmpb $0x05, %dh
je 1f /* Is +-Inf, jump. */
fld %st
movzwl 4+8(%esp), %eax
andl $0x7fff, %eax
cmpl $0x3fbe, %eax
jge 3f
/* Argument's exponent below -65, result rounds to 1. */
fld1
faddp
ret
3: fld %st
frndint /* int(x) */
fsubr %st,%st(1) /* fract(x) */
fxch

View File

@@ -61,6 +61,9 @@ __ieee754_exp2 (double x)
double rx, x22, result;
union ieee754_double ex2_u, scale_u;
if (fabs (x) < DBL_EPSILON / 4.0)
return 1.0 + x;
{
SET_RESTORE_ROUND_NOEX (FE_TONEAREST);

View File

@@ -54,6 +54,9 @@ __ieee754_exp2f (float x)
float rx, x22, result;
union ieee754_float ex2_u, scale_u;
if (fabsf (x) < FLT_EPSILON / 4.0f)
return 1.0f + x;
{
SET_RESTORE_ROUND_NOEXF (FE_TONEAREST);

View File

@@ -19,7 +19,15 @@ ENTRY(__ieee754_exp2l)
andb %ah, %dh
cmpb $0x05, %dh
je 1f /* Is +-Inf, jump. */
fld %st
movzwl 8+8(%rsp), %eax
andl $0x7fff, %eax
cmpl $0x3fbe, %eax
jge 3f
/* Argument's exponent below -65, result rounds to 1. */
fld1
faddp
ret
3: fld %st
frndint /* int(x) */
fsubr %st,%st(1) /* fract(x) */
fxch