Re: Fix hypotf overflow/underflow (bug 13840)

[David Miller <davem at davemloft dot net>]

From: David Miller <davem@davemloft.net>
Date: Tue, 13 Mar 2012 15:43:25 -0700 (PDT)

> By using a cast to double all the scaling code can be removed, and all
> that's left are the checks against 0x7f800000 and zero.

Something like this.  It passes all the existing math tests at
least on sparc.  I'll try with the new tests from your patch
Joseph.

diff --git a/sysdeps/ieee754/flt-32/e_hypotf.c b/sysdeps/ieee754/flt-32/e_hypotf.c
index 97e7d34..8f499b5 100644
--- a/sysdeps/ieee754/flt-32/e_hypotf.c
+++ b/sysdeps/ieee754/flt-32/e_hypotf.c
@@ -19,62 +19,35 @@
 float
 __ieee754_hypotf(float x, float y)
 {
-	float a,b,t1,t2,y1,y2,w;
-	int32_t j,k,ha,hb;
+	double d_x, d_y;
+	int32_t ha, hb;
 
 	GET_FLOAT_WORD(ha,x);
 	ha &= 0x7fffffff;
 	GET_FLOAT_WORD(hb,y);
 	hb &= 0x7fffffff;
-	if(hb > ha) {a=y;b=x;j=ha; ha=hb;hb=j;} else {a=x;b=y;}
-	SET_FLOAT_WORD(a,ha);	/* a <- |a| */
-	SET_FLOAT_WORD(b,hb);	/* b <- |b| */
-	if((ha-hb)>0xf000000) {return a+b;} /* x/y > 2**30 */
-	k=0;
-	if(__builtin_expect(ha > 0x58800000, 0)) {	/* a>2**50 */
-	   if(ha >= 0x7f800000) {	/* Inf or NaN */
-	       w = a+b;			/* for sNaN */
-	       if(ha == 0x7f800000) w = a;
-	       if(hb == 0x7f800000) w = b;
-	       return w;
-	   }
-	   /* scale a and b by 2**-60 */
-	   ha -= 0x1e000000; hb -= 0x1e000000;	k += 60;
-	   SET_FLOAT_WORD(a,ha);
-	   SET_FLOAT_WORD(b,hb);
-	}
-	if(__builtin_expect(hb < 0x26800000, 0)) {	/* b < 2**-50 */
-	    if(hb <= 0x007fffff) {	/* subnormal b or 0 */
-		if(hb==0) return a;
-		SET_FLOAT_WORD(t1,0x7e800000);	/* t1=2^126 */
-		b *= t1;
-		a *= t1;
-		k -= 126;
-	    } else {		/* scale a and b by 2^60 */
-		ha += 0x1e000000;	/* a *= 2^60 */
-		hb += 0x1e000000;	/* b *= 2^60 */
-		k -= 60;
-		SET_FLOAT_WORD(a,ha);
-		SET_FLOAT_WORD(b,hb);
-	    }
-	}
-    /* medium size a and b */
-	w = a-b;
-	if (w>b) {
-	    SET_FLOAT_WORD(t1,ha&0xfffff000);
-	    t2 = a-t1;
-	    w  = __ieee754_sqrtf(t1*t1-(b*(-b)-t2*(a+t1)));
-	} else {
-	    a  = a+a;
-	    SET_FLOAT_WORD(y1,hb&0xfffff000);
-	    y2 = b - y1;
-	    SET_FLOAT_WORD(t1,ha+0x00800000);
-	    t2 = a - t1;
-	    w  = __ieee754_sqrtf(t1*y1-(w*(-w)-(t1*y2+t2*b)));
-	}
-	if(k!=0) {
-	    SET_FLOAT_WORD(t1,0x3f800000+(k<<23));
-	    return t1*w;
-	} else return w;
+	if (ha == 0x7f800000)
+	  {
+	    if (x == y)
+	      return fabsf(y);
+	    return fabsf(x);
+	  }
+	else if (hb == 0x7f800000)
+	  {
+	    if (x == y)
+	      return fabsf(x);
+	    return fabsf(y);
+	  }
+	else if (ha > 0x7f800000 || hb > 0x7f800000)
+	  return fabsf(x) * fabsf(y);
+	else if (ha == 0)
+	  return fabsf(y);
+	else if (hb == 0)
+	  return fabsf(x);
+
+	d_x = (double) x;
+	d_y = (double) y;
+
+	return (float) __ieee754_sqrt(d_x * d_x + d_y * d_y);
 }
 strong_alias (__ieee754_hypotf, __hypotf_finite)