From 6ace39382120dd0668b2168591611f71b04ac8f3 Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Fri, 25 Sep 2015 22:29:10 +0000 Subject: Fix pow missing underflows (bug 18825). 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. [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 . 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 . 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. --- sysdeps/i386/fpu/e_pow.S | 8 +++++-- sysdeps/i386/fpu/e_powf.S | 13 +++++++--- sysdeps/i386/fpu/e_powl.S | 6 +++++ sysdeps/i386/fpu/i386-math-asm.h | 46 ++++++++++++++++++++++++++++++++++++ sysdeps/ieee754/dbl-64/e_pow.c | 2 ++ sysdeps/ieee754/flt-32/e_powf.c | 7 +++++- sysdeps/ieee754/ldbl-128/e_powl.c | 6 ++++- sysdeps/ieee754/ldbl-128ibm/e_powl.c | 1 + sysdeps/x86/fpu/powl_helper.c | 1 + sysdeps/x86_64/fpu/e_powl.S | 3 +++ sysdeps/x86_64/fpu/x86_64-math-asm.h | 13 ++++++++++ 11 files changed, 99 insertions(+), 7 deletions(-) (limited to 'sysdeps') diff --git a/sysdeps/i386/fpu/e_pow.S b/sysdeps/i386/fpu/e_pow.S index 2903e13e59..8203c9625a 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 d929bdc9f1..6608d3e9d1 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 bdfb5cf635..76ca8213fd 100644 --- a/sysdeps/i386/fpu/e_powl.S +++ b/sysdeps/i386/fpu/e_powl.S @@ -18,6 +18,7 @@ . */ #include +#include .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 5d8d847cc3..4bac3bf5ab 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 3c027fe7f5..3b0bbe3103 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 8e8d918456..18042960c2 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 f531385232..7f3037fb51 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 c942f2f249..90340e890e 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 f2774b7059..7f46c7327a 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 358abb8dcb..a44c311bb6 100644 --- a/sysdeps/x86_64/fpu/e_powl.S +++ b/sysdeps/x86_64/fpu/e_powl.S @@ -18,6 +18,7 @@ . */ #include +#include .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 07eecf7079..b5482518d7 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); \ -- cgit 1.4.1