From f6987f5aa4c2ef7270bcb0a828808bd3db53520b Mon Sep 17 00:00:00 2001 From: Joseph Myers Date: Thu, 24 Sep 2015 23:43:57 +0000 Subject: Fix hypot missing underflows (bug 18803). Similar to various other bugs in this area, hypot 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. Note that this issue cannot arise for implementations of hypotf using double (or wider) for intermediate evaluation (if hypotf should underflow, that means the double square root is being computed of some number of the form N*2^-298, for 0 < N < 2^46, which is exactly represented as a double, and whatever the rounding mode such a square root cannot have a mantissa with all zeroes after the initial 23 bits). Thus no changes are made to hypotf implementations in this patch, only to hypot and hypotl. Tested for x86_64, x86, mips64 and powerpc. [BZ #18803] * sysdeps/i386/fpu/e_hypot.S: Use DEFINE_DBL_MIN. (MO): New macro. (__ieee754_hypot) [PIC]: Load PIC register. (__ieee754_hypot): Use DBL_NARROW_EVAL_UFLOW_NONNEG instead of DBL_NARROW_EVAL. * sysdeps/ieee754/dbl-64/e_hypot.c (__ieee754_hypot): Use math_check_force_underflow_nonneg in case where result might be tiny. * sysdeps/ieee754/ldbl-128/e_hypotl.c (__ieee754_hypotl): Likewise. * sysdeps/ieee754/ldbl-128ibm/e_hypotl.c (__ieee754_hypotl): Likewise. * sysdeps/ieee754/ldbl-96/e_hypotl.c (__ieee754_hypotl): Likewise. * sysdeps/powerpc/fpu/e_hypot.c (__ieee754_hypot): Likewise. * math/auto-libm-test-in: Add more tests of hypot. * math/auto-libm-test-out: Regenerated. --- sysdeps/i386/fpu/e_hypot.S | 13 ++++++++++++- sysdeps/ieee754/dbl-64/e_hypot.c | 4 +++- sysdeps/ieee754/ldbl-128/e_hypotl.c | 4 +++- sysdeps/ieee754/ldbl-128ibm/e_hypotl.c | 6 +++++- sysdeps/ieee754/ldbl-96/e_hypotl.c | 4 +++- sysdeps/powerpc/fpu/e_hypot.c | 4 +++- 6 files changed, 29 insertions(+), 6 deletions(-) (limited to 'sysdeps') diff --git a/sysdeps/i386/fpu/e_hypot.S b/sysdeps/i386/fpu/e_hypot.S index 5323fde01e..e69ea096ab 100644 --- a/sysdeps/i386/fpu/e_hypot.S +++ b/sysdeps/i386/fpu/e_hypot.S @@ -20,8 +20,19 @@ #include #include +DEFINE_DBL_MIN + +#ifdef PIC +# define MO(op) op##@GOTOFF(%edx) +#else +# define MO(op) op +#endif + .text ENTRY(__ieee754_hypot) +#ifdef PIC + LOAD_PIC_REG (dx) +#endif fldl 4(%esp) // x fxam fnstsw @@ -38,7 +49,7 @@ ENTRY(__ieee754_hypot) fmul %st(0) // x * x : y * y faddp // x * x + y * y fsqrt - DBL_NARROW_EVAL + DBL_NARROW_EVAL_UFLOW_NONNEG 2: ret // We have to test whether any of the parameters is Inf. diff --git a/sysdeps/ieee754/dbl-64/e_hypot.c b/sysdeps/ieee754/dbl-64/e_hypot.c index 5cbfcbeb48..f142c450a2 100644 --- a/sysdeps/ieee754/dbl-64/e_hypot.c +++ b/sysdeps/ieee754/dbl-64/e_hypot.c @@ -149,7 +149,9 @@ __ieee754_hypot (double x, double y) t1 = 1.0; GET_HIGH_WORD (high, t1); SET_HIGH_WORD (t1, high + (k << 20)); - return t1 * w; + w *= t1; + math_check_force_underflow_nonneg (w); + return w; } else return w; diff --git a/sysdeps/ieee754/ldbl-128/e_hypotl.c b/sysdeps/ieee754/ldbl-128/e_hypotl.c index 01444cfb4e..80e5e38c72 100644 --- a/sysdeps/ieee754/ldbl-128/e_hypotl.c +++ b/sysdeps/ieee754/ldbl-128/e_hypotl.c @@ -130,7 +130,9 @@ __ieee754_hypotl(long double x, long double y) t1 = 1.0L; GET_LDOUBLE_MSW64(high,t1); SET_LDOUBLE_MSW64(t1,high+(k<<48)); - return t1*w; + w *= t1; + math_check_force_underflow_nonneg (w); + return w; } else return w; } strong_alias (__ieee754_hypotl, __hypotl_finite) diff --git a/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c b/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c index 3b07a47b40..c68dac03b0 100644 --- a/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c +++ b/sysdeps/ieee754/ldbl-128ibm/e_hypotl.c @@ -125,7 +125,11 @@ __ieee754_hypotl(long double x, long double y) w = __ieee754_sqrtl(a1*b1-(w*(-w)-(a1*b2+a2*b))); } if(k!=0) - return w*kld; + { + w *= kld; + math_check_force_underflow_nonneg (w); + return w; + } else return w; } diff --git a/sysdeps/ieee754/ldbl-96/e_hypotl.c b/sysdeps/ieee754/ldbl-96/e_hypotl.c index d3152f91e5..ee3a07055b 100644 --- a/sysdeps/ieee754/ldbl-96/e_hypotl.c +++ b/sysdeps/ieee754/ldbl-96/e_hypotl.c @@ -132,7 +132,9 @@ long double __ieee754_hypotl(long double x, long double y) t1 = 1.0; GET_LDOUBLE_EXP(exp,t1); SET_LDOUBLE_EXP(t1,exp+k); - return t1*w; + w *= t1; + math_check_force_underflow_nonneg (w); + return w; } else return w; } strong_alias (__ieee754_hypotl, __hypotl_finite) diff --git a/sysdeps/powerpc/fpu/e_hypot.c b/sysdeps/powerpc/fpu/e_hypot.c index 98e4ae635f..a3b9d3d371 100644 --- a/sysdeps/powerpc/fpu/e_hypot.c +++ b/sysdeps/powerpc/fpu/e_hypot.c @@ -116,7 +116,9 @@ __ieee754_hypot (double x, double y) { x *= two1022; y *= two1022; - return __ieee754_sqrt (x * x + y * y) / two1022; + double ret = __ieee754_sqrt (x * x + y * y) / two1022; + math_check_force_underflow_nonneg (ret); + return ret; } else { -- cgit 1.4.1