This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
Fix pow missing underflows (bug 18825) [committed]
- From: Joseph Myers <joseph at codesourcery dot com>
- To: <libc-alpha at sourceware dot org>
- Date: Fri, 25 Sep 2015 22:30:02 +0000
- Subject: Fix pow missing underflows (bug 18825) [committed]
- Authentication-results: sourceware.org; auth=none
Similar to various other bugs in this area, pow functions can fail to
raise the underflow exception when the result is tiny and inexact but
one or more low bits of the intermediate result that is scaled down
(or, in the i386 case, converted from a wider evaluation format) are
zero. This patch forces the exception in a similar way to previous
fixes, thereby concluding the fixes for known bugs with missing
underflow exceptions currently filed in Bugzilla.
Tested for x86_64, x86, mips64 and powerpc. Committed.
(auto-libm-test-out diffs omitted below.)
2015-09-25 Joseph Myers <joseph@codesourcery.com>
[BZ #18825]
* sysdeps/i386/fpu/i386-math-asm.h (FLT_NARROW_EVAL_UFLOW_NONNAN):
New macro.
(DBL_NARROW_EVAL_UFLOW_NONNAN): Likewise.
(LDBL_CHECK_FORCE_UFLOW_NONNAN): Likewise.
* sysdeps/i386/fpu/e_pow.S: Use DEFINE_DBL_MIN.
(__ieee754_pow): Use DBL_NARROW_EVAL_UFLOW_NONNAN instead of
DBL_NARROW_EVAL, reloading the PIC register as needed.
* sysdeps/i386/fpu/e_powf.S: Use DEFINE_FLT_MIN.
(__ieee754_powf): Use FLT_NARROW_EVAL_UFLOW_NONNAN instead of
FLT_NARROW_EVAL. Use separate return path for case when first
argument is NaN.
* sysdeps/i386/fpu/e_powl.S: Include <i386-math-asm.h>. Use
DEFINE_LDBL_MIN.
(__ieee754_powl): Use LDBL_CHECK_FORCE_UFLOW_NONNAN, reloading the
PIC register.
* sysdeps/ieee754/dbl-64/e_pow.c (__ieee754_pow): Use
math_check_force_underflow_nonneg.
* sysdeps/ieee754/flt-32/e_powf.c (__ieee754_powf): Force
underflow for subnormal result.
* sysdeps/ieee754/ldbl-128/e_powl.c (__ieee754_powl): Likewise.
* sysdeps/ieee754/ldbl-128ibm/e_powl.c (__ieee754_powl): Use
math_check_force_underflow_nonneg.
* sysdeps/x86/fpu/powl_helper.c (__powl_helper): Use
math_check_force_underflow.
* sysdeps/x86_64/fpu/x86_64-math-asm.h
(LDBL_CHECK_FORCE_UFLOW_NONNAN): New macro.
* sysdeps/x86_64/fpu/e_powl.S: Include <x86_64-math-asm.h>. Use
DEFINE_LDBL_MIN.
(__ieee754_powl): Use LDBL_CHECK_FORCE_UFLOW_NONNAN.
* math/auto-libm-test-in: Add more tests of pow.
* math/auto-libm-test-out: Regenerated.
diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 2403216..4b93857 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -3067,6 +3067,151 @@ pow 1.01 -min_subnorm
pow 2.0 -100000.0
+pow 2 -126.125
+pow 2 -126.25
+pow 2 -126.375
+pow 2 -126.5
+pow 2 -126.625
+pow 2 -126.75
+pow 2 -126.875
+pow 2 -969.125
+pow 2 -969.25
+pow 2 -969.375
+pow 2 -969.5
+pow 2 -969.625
+pow 2 -969.75
+pow 2 -969.875
+pow 2 -1022.125
+pow 2 -1022.25
+pow 2 -1022.375
+pow 2 -1022.5
+pow 2 -1022.625
+pow 2 -1022.75
+pow 2 -1022.875
+pow 2 -16382.125
+pow 2 -16382.25
+pow 2 -16382.375
+pow 2 -16382.5
+pow 2 -16382.625
+pow 2 -16382.75
+pow 2 -16382.875
+pow 2 -16383.125
+pow 2 -16383.25
+pow 2 -16383.375
+pow 2 -16383.5
+pow 2 -16383.625
+pow 2 -16383.75
+pow 2 -16383.875
+pow 0.5 126.125
+pow 0.5 126.25
+pow 0.5 126.375
+pow 0.5 126.5
+pow 0.5 126.625
+pow 0.5 126.75
+pow 0.5 126.875
+pow 0.5 969.125
+pow 0.5 969.25
+pow 0.5 969.375
+pow 0.5 969.5
+pow 0.5 969.625
+pow 0.5 969.75
+pow 0.5 969.875
+pow 0.5 1022.125
+pow 0.5 1022.25
+pow 0.5 1022.375
+pow 0.5 1022.5
+pow 0.5 1022.625
+pow 0.5 1022.75
+pow 0.5 1022.875
+pow 0.5 16382.125
+pow 0.5 16382.25
+pow 0.5 16382.375
+pow 0.5 16382.5
+pow 0.5 16382.625
+pow 0.5 16382.75
+pow 0.5 16382.875
+pow 0.5 16383.125
+pow 0.5 16383.25
+pow 0.5 16383.375
+pow 0.5 16383.5
+pow 0.5 16383.625
+pow 0.5 16383.75
+pow 0.5 16383.875
+pow 0x1.00001p1 -126
+pow -0x1.00002p1 -126
+pow 0x1.00003p1 -126
+pow -0x1.00004p1 -126
+pow 0x1.00005p1 -126
+pow -0x1.00006p1 -126
+pow 0x1.00007p1 -126
+pow 0x1.00001p1 -969
+pow -0x1.00002p1 -969
+pow 0x1.00003p1 -969
+pow -0x1.00004p1 -969
+pow 0x1.00005p1 -969
+pow -0x1.00006p1 -969
+pow 0x1.00007p1 -969
+pow 0x1.00001p1 -1022
+pow -0x1.00002p1 -1022
+pow 0x1.00003p1 -1022
+pow -0x1.00004p1 -1022
+pow 0x1.00005p1 -1022
+pow -0x1.00006p1 -1022
+pow 0x1.00007p1 -1022
+pow 0x1.00001p1 -16382
+pow -0x1.00002p1 -16382
+pow 0x1.00003p1 -16382
+pow -0x1.00004p1 -16382
+pow 0x1.00005p1 -16382
+pow -0x1.00006p1 -16382
+pow 0x1.00007p1 -16382
+pow 0x1.00001p1 -16383
+pow -0x1.00002p1 -16383
+pow 0x1.00003p1 -16383
+pow -0x1.00004p1 -16383
+pow 0x1.00005p1 -16383
+pow -0x1.00006p1 -16383
+pow 0x1.00007p1 -16383
+pow 0x0.ffff1p-1 126
+pow -0x0.ffff2p-1 126
+pow 0x0.ffff3p-1 126
+pow -0x0.ffff4p-1 126
+pow 0x0.ffff5p-1 126
+pow -0x0.ffff6p-1 126
+pow 0x0.ffff7p-1 126
+pow 0x0.ffff1p-1 969
+pow -0x0.ffff2p-1 969
+pow 0x0.ffff3p-1 969
+pow -0x0.ffff4p-1 969
+pow 0x0.ffff5p-1 969
+pow -0x0.ffff6p-1 969
+pow 0x0.ffff7p-1 969
+pow 0x0.ffff1p-1 1022
+pow -0x0.ffff2p-1 1022
+pow 0x0.ffff3p-1 1022
+pow -0x0.ffff4p-1 1022
+pow 0x0.ffff5p-1 1022
+pow -0x0.ffff6p-1 1022
+pow 0x0.ffff7p-1 1022
+pow 0x0.ffff1p-1 16382
+pow -0x0.ffff2p-1 16382
+pow 0x0.ffff3p-1 16382
+pow -0x0.ffff4p-1 16382
+pow 0x0.ffff5p-1 16382
+pow -0x0.ffff6p-1 16382
+pow 0x0.ffff7p-1 16382
+pow 0x0.ffff1p-1 16383
+pow -0x0.ffff2p-1 16383
+pow 0x0.ffff3p-1 16383
+pow -0x0.ffff4p-1 16383
+pow 0x0.ffff5p-1 16383
+pow -0x0.ffff6p-1 16383
+pow 0x0.ffff7p-1 16383
+pow 0x2.000b3p0 -1022
+pow 0x2.000582p0 -1022
+pow 2 -0x3.fe513p+8
+pow 2 -0x3.fe4e8p+8
+
pow 1.0625 1.125
pow 1.5 1.03125
pow 0x1.7d1a0a6f2p+681 1.5
diff --git a/sysdeps/i386/fpu/e_pow.S b/sysdeps/i386/fpu/e_pow.S
index 2903e13..8203c96 100644
--- a/sysdeps/i386/fpu/e_pow.S
+++ b/sysdeps/i386/fpu/e_pow.S
@@ -54,6 +54,7 @@ minfinity:
mzero:
.byte 0, 0, 0, 0, 0, 0, 0, 0x80
ASM_SIZE_DIRECTIVE(minf_mzero)
+DEFINE_DBL_MIN
#ifdef PIC
# define MO(op) op##@GOTOFF(%ecx)
@@ -166,7 +167,10 @@ ENTRY(__ieee754_pow)
orl %edx, %ecx
jnz 6b
fstp %st(0) // ST*x
- DBL_NARROW_EVAL
+#ifdef PIC
+ LOAD_PIC_REG (cx)
+#endif
+ DBL_NARROW_EVAL_UFLOW_NONNAN
ret
/* y is ?NAN */
@@ -259,7 +263,7 @@ ENTRY(__ieee754_pow)
cfi_adjust_cfa_offset (-8)
292: fscale // +/- 2^fract(y*log2(x))*2^int(y*log2(x)) : int(y*log2(x))
fstp %st(1) // +/- 2^fract(y*log2(x))*2^int(y*log2(x))
- DBL_NARROW_EVAL
+ DBL_NARROW_EVAL_UFLOW_NONNAN
ret
diff --git a/sysdeps/i386/fpu/e_powf.S b/sysdeps/i386/fpu/e_powf.S
index d929bdc..6608d3e 100644
--- a/sysdeps/i386/fpu/e_powf.S
+++ b/sysdeps/i386/fpu/e_powf.S
@@ -51,6 +51,7 @@ minfinity:
mzero:
.byte 0, 0, 0, 0, 0, 0, 0, 0x80
ASM_SIZE_DIRECTIVE(minf_mzero)
+DEFINE_FLT_MIN
#ifdef PIC
# define MO(op) op##@GOTOFF(%ecx)
@@ -97,7 +98,7 @@ ENTRY(__ieee754_powf)
je 15f // x is ?inf
cmpb $0x01, %ah
- je 32f // x is NaN
+ je 33f // x is NaN
fxch // y : x
@@ -149,7 +150,7 @@ ENTRY(__ieee754_powf)
testl %edx, %edx
jnz 6b
fstp %st(0) // ST*x
- FLT_NARROW_EVAL
+ FLT_NARROW_EVAL_UFLOW_NONNAN
ret
/* y is ?NAN */
@@ -198,9 +199,15 @@ ENTRY(__ieee754_powf)
32: addl $4, %esp
cfi_adjust_cfa_offset (-4)
fstp %st(1) // 2^fract(y*log2(x))*2^int(y*log2(x))
- FLT_NARROW_EVAL
+ FLT_NARROW_EVAL_UFLOW_NONNAN
ret
+ /* x is NaN. */
+ cfi_adjust_cfa_offset (4)
+33: addl $4, %esp
+ cfi_adjust_cfa_offset (-4)
+ fstp %st(1)
+ ret
// pow(x,?0) = 1
.align ALIGNARG(4)
diff --git a/sysdeps/i386/fpu/e_powl.S b/sysdeps/i386/fpu/e_powl.S
index bdfb5cf..76ca821 100644
--- a/sysdeps/i386/fpu/e_powl.S
+++ b/sysdeps/i386/fpu/e_powl.S
@@ -18,6 +18,7 @@
<http://www.gnu.org/licenses/>. */
#include <machine/asm.h>
+#include <i386-math-asm.h>
.section .rodata.cst8,"aM",@progbits,8
@@ -59,6 +60,7 @@ minfinity:
mzero:
.byte 0, 0, 0, 0, 0, 0, 0, 0x80
ASM_SIZE_DIRECTIVE(minf_mzero)
+DEFINE_LDBL_MIN
#ifdef PIC
# define MO(op) op##@GOTOFF(%ecx)
@@ -187,6 +189,10 @@ ENTRY(__ieee754_powl)
orl %edx, %ecx
jnz 6b
fstp %st(0) // ST*x
+#ifdef PIC
+ LOAD_PIC_REG (cx)
+#endif
+ LDBL_CHECK_FORCE_UFLOW_NONNAN
ret
/* y is ?NAN */
diff --git a/sysdeps/i386/fpu/i386-math-asm.h b/sysdeps/i386/fpu/i386-math-asm.h
index 5d8d847..4bac3bf 100644
--- a/sysdeps/i386/fpu/i386-math-asm.h
+++ b/sysdeps/i386/fpu/i386-math-asm.h
@@ -131,6 +131,40 @@ ldbl_min: \
addl $8, %esp; \
cfi_adjust_cfa_offset (-8);
+/* Likewise, but the non-NaN argument may be negative. */
+#define FLT_NARROW_EVAL_UFLOW_NONNAN \
+ subl $4, %esp; \
+ cfi_adjust_cfa_offset (4); \
+ fld %st(0); \
+ fabs; \
+ fcomps MO(flt_min); \
+ fnstsw; \
+ sahf; \
+ jnc 6424f; \
+ fld %st(0); \
+ fmul %st(0); \
+ fstps (%esp); \
+6424: fstps (%esp); \
+ flds (%esp); \
+ addl $4, %esp; \
+ cfi_adjust_cfa_offset (-4);
+#define DBL_NARROW_EVAL_UFLOW_NONNAN \
+ subl $8, %esp; \
+ cfi_adjust_cfa_offset (8); \
+ fld %st(0); \
+ fabs; \
+ fcompl MO(dbl_min); \
+ fnstsw; \
+ sahf; \
+ jnc 6453f; \
+ fld %st(0); \
+ fmul %st(0); \
+ fstpl (%esp); \
+6453: fstpl (%esp); \
+ fldl (%esp); \
+ addl $8, %esp; \
+ cfi_adjust_cfa_offset (-8);
+
/* Force an underflow exception if the given value is subnormal. The
relevant constant for the minimum of the type must have been
defined, the MO macro must have been defined for access to memory
@@ -251,6 +285,18 @@ ldbl_min: \
addl $8, %esp; \
cfi_adjust_cfa_offset (-8); \
6453:
+#define LDBL_CHECK_FORCE_UFLOW_NONNAN \
+ fldt MO(ldbl_min); \
+ fld %st(1); \
+ fabs; \
+ fcompp; \
+ fnstsw; \
+ sahf; \
+ jnc 6464f; \
+ fld %st(0); \
+ fmul %st(0); \
+ fstp %st(0); \
+6464:
/* Likewise, but the argument is nonnegative and not a NaN. */
#define FLT_CHECK_FORCE_UFLOW_NONNEG \
diff --git a/sysdeps/ieee754/dbl-64/e_pow.c b/sysdeps/ieee754/dbl-64/e_pow.c
index 3c027fe..3b0bbe3 100644
--- a/sysdeps/ieee754/dbl-64/e_pow.c
+++ b/sysdeps/ieee754/dbl-64/e_pow.c
@@ -119,6 +119,8 @@ __ieee754_pow (double x, double y)
retval = huge * huge;
else if (retval == 0)
retval = tiny * tiny;
+ else
+ math_check_force_underflow_nonneg (retval);
return retval;
}
diff --git a/sysdeps/ieee754/flt-32/e_powf.c b/sysdeps/ieee754/flt-32/e_powf.c
index 8e8d918..1804296 100644
--- a/sysdeps/ieee754/flt-32/e_powf.c
+++ b/sysdeps/ieee754/flt-32/e_powf.c
@@ -244,7 +244,12 @@ __ieee754_powf(float x, float y)
z = one-(r-z);
GET_FLOAT_WORD(j,z);
j += (n<<23);
- if((j>>23)<=0) z = __scalbnf(z,n); /* subnormal output */
+ if((j>>23)<=0) /* subnormal output */
+ {
+ z = __scalbnf (z, n);
+ float force_underflow = z * z;
+ math_force_eval (force_underflow);
+ }
else SET_FLOAT_WORD(z,j);
return s*z;
}
diff --git a/sysdeps/ieee754/ldbl-128/e_powl.c b/sysdeps/ieee754/ldbl-128/e_powl.c
index f531385..7f3037f 100644
--- a/sysdeps/ieee754/ldbl-128/e_powl.c
+++ b/sysdeps/ieee754/ldbl-128/e_powl.c
@@ -435,7 +435,11 @@ __ieee754_powl (long double x, long double y)
j = o.parts32.w0;
j += (n << 16);
if ((j >> 16) <= 0)
- z = __scalbnl (z, n); /* subnormal output */
+ {
+ z = __scalbnl (z, n); /* subnormal output */
+ long double force_underflow = z * z;
+ math_force_eval (force_underflow);
+ }
else
{
o.parts32.w0 = j;
diff --git a/sysdeps/ieee754/ldbl-128ibm/e_powl.c b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
index c942f2f..90340e8 100644
--- a/sysdeps/ieee754/ldbl-128ibm/e_powl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/e_powl.c
@@ -409,6 +409,7 @@ __ieee754_powl (long double x, long double y)
r = (z * t1) / (t1 - two) - (w + z * w);
z = one - (r - z);
z = __scalbnl (z, n);
+ math_check_force_underflow_nonneg (z);
return s * z;
}
strong_alias (__ieee754_powl, __powl_finite)
diff --git a/sysdeps/x86/fpu/powl_helper.c b/sysdeps/x86/fpu/powl_helper.c
index f2774b7..7f46c73 100644
--- a/sysdeps/x86/fpu/powl_helper.c
+++ b/sysdeps/x86/fpu/powl_helper.c
@@ -229,6 +229,7 @@ __powl_helper (long double x, long double y)
if (negate)
res = -res;
asm ("fscale" : "=t" (res) : "0" (res), "u" (log2_res_int));
+ math_check_force_underflow (res);
return res;
}
diff --git a/sysdeps/x86_64/fpu/e_powl.S b/sysdeps/x86_64/fpu/e_powl.S
index 358abb8..a44c311 100644
--- a/sysdeps/x86_64/fpu/e_powl.S
+++ b/sysdeps/x86_64/fpu/e_powl.S
@@ -18,6 +18,7 @@
<http://www.gnu.org/licenses/>. */
#include <machine/asm.h>
+#include <x86_64-math-asm.h>
.section .rodata.cst8,"aM",@progbits,8
@@ -59,6 +60,7 @@ minfinity:
mzero:
.byte 0, 0, 0, 0, 0, 0, 0, 0x80
ASM_SIZE_DIRECTIVE(minf_mzero)
+DEFINE_LDBL_MIN
#ifdef PIC
# define MO(op) op##(%rip)
@@ -175,6 +177,7 @@ ENTRY(__ieee754_powl)
orl %edx, %ecx
jnz 6b
fstp %st(0) // ST*x
+ LDBL_CHECK_FORCE_UFLOW_NONNAN
ret
/* y is ?NAN */
diff --git a/sysdeps/x86_64/fpu/x86_64-math-asm.h b/sysdeps/x86_64/fpu/x86_64-math-asm.h
index 07eecf7..b548251 100644
--- a/sysdeps/x86_64/fpu/x86_64-math-asm.h
+++ b/sysdeps/x86_64/fpu/x86_64-math-asm.h
@@ -46,6 +46,19 @@ ldbl_min: \
fstp %st(0); \
6464:
+/* Likewise, but the argument is not a NaN. */
+#define LDBL_CHECK_FORCE_UFLOW_NONNAN \
+ fldt MO(ldbl_min); \
+ fld %st(1); \
+ fabs; \
+ fcomip %st(1), %st(0); \
+ fstp %st(0); \
+ jnc 6464f; \
+ fld %st(0); \
+ fmul %st(0); \
+ fstp %st(0); \
+6464:
+
/* Likewise, but the argument is nonnegative and not a NaN. */
#define LDBL_CHECK_FORCE_UFLOW_NONNEG \
fldt MO(ldbl_min); \
--
Joseph S. Myers
joseph@codesourcery.com