This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

Fix hypot missing underflows (bug 18803) [committed]


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.  Committed.

(auto-libm-test-out diffs omitted below.)

2015-09-24  Joseph Myers  <joseph@codesourcery.com>

	[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.

diff --git a/math/auto-libm-test-in b/math/auto-libm-test-in
index 8dea2b4..2403216 100644
--- a/math/auto-libm-test-in
+++ b/math/auto-libm-test-in
@@ -1955,6 +1955,13 @@ hypot 0x0.ffffffp-16382 0x0.fp-16383 no-test-inline
 hypot 0x0.ffffffp-16382 0x0.fp-16386 no-test-inline
 hypot 0 min_subnorm no-test-inline
 
+hypot 0x1p-127 0x1p-149
+hypot 0x1p-1023 0x1p-1074
+hypot 0x1p-970 0x1p-1074
+hypot 0x1p-16383 0x1p-16445 no-test-inline
+hypot 0x1p-16384 0x1p-16446 no-test-inline
+hypot 0x1p-16383 0x1p-16494 no-test-inline
+
 hypot -0x1.fa7deap+0 0x1.a761bab383ac8p+0
 
 j0 -1.0
diff --git a/sysdeps/i386/fpu/e_hypot.S b/sysdeps/i386/fpu/e_hypot.S
index 5323fde..e69ea09 100644
--- a/sysdeps/i386/fpu/e_hypot.S
+++ b/sysdeps/i386/fpu/e_hypot.S
@@ -20,8 +20,19 @@
 #include <sysdep.h>
 #include <i386-math-asm.h>
 
+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 5cbfcbe..f142c45 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 01444cf..80e5e38 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 3b07a47..c68dac0 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 d3152f9..ee3a070 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 98e4ae6..a3b9d3d 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
 	{

-- 
Joseph S. Myers
joseph@codesourcery.com


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]