RFA: patch to fix scaling in ef_hypot.c

[Richard Sandiford]

This patch fixes a scaling bug in ef_hypot.c.  After checking that
the arguments are close enough together for a sqrt to be needed,
it uses the following code to rescale them:

	if(ha > 0x58800000L) {	/* a>2**50 */
	   if(!FLT_UWORD_IS_FINITE(ha)) {	/* Inf or NaN */
	       w = a+b;			/* for sNaN */
	       if(FLT_UWORD_IS_INFINITE(ha)) w = a;
	       if(FLT_UWORD_IS_INFINITE(hb)) w = b;
	       return w;
	   }
	   /* scale a and b by 2**-60 */
	   ha -= 0x5d800000L; hb -= 0x5d800000L;	k += 60;
	   SET_FLOAT_WORD(a,ha);
	   SET_FLOAT_WORD(b,hb);
	}

(similar code for b<2**50).

There's two problems here.

 (a) Scaling by 2**-60 isn't enough if both arguments are near the
     extreme end of the range, such as in hypotf (2**125, 2**125).
     We end up multiplying two numbers O(2**65), which overflows.

 (b) More importantly, the code adjusts "ha" and "hb" by the float
     representation of 2**60 (with biased exponent) rather than
     simply subtracting 60 from the exponent field.  The exponents
     of the new numbers are far too small.

The patch fixes (a) by using an adjustment of 80 and changes the "ha"
and "hb" adjustments to match.

I used the attached program as a test case.  It checks the results of
hypotf() against the double version, hypot(), assuming it to be fairly
trustworthy.  I compared the program's output before and after the
patch on a mips64-elf target.  The patch seemed to be a strict improvement:
some things that failed before now pass, and nothing now fails that
passed before.

OK to apply?

Note: there seems to be problems with the handling of denormals,
but this patch doesn't make it any worse.

Richard

	* libm/math/ef_hypot.c: Increasing scale factor to 80.

Index: libm/math/ef_hypot.c
===================================================================
RCS file: /cvs/cvsfiles/devo/newlib/libm/math/ef_hypot.c,v
retrieving revision 1.5
diff -c -p -d -r1.5 ef_hypot.c
*** libm/math/ef_hypot.c	2001/04/04 13:36:56	1.5
--- libm/math/ef_hypot.c	2002/02/20 16:38:43
***************
*** 41,48 ****
  	       if(FLT_UWORD_IS_INFINITE(hb)) w = b;
  	       return w;
  	   }
! 	   /* scale a and b by 2**-60 */
! 	   ha -= 0x5d800000L; hb -= 0x5d800000L;	k += 60;
  	   SET_FLOAT_WORD(a,ha);
  	   SET_FLOAT_WORD(b,hb);
  	}
--- 41,48 ----
  	       if(FLT_UWORD_IS_INFINITE(hb)) w = b;
  	       return w;
  	   }
! 	   /* scale a and b by 2**-80 */
! 	   ha -= 0x28000000L; hb -= 0x28000000L;	k += 80;
  	   SET_FLOAT_WORD(a,ha);
  	   SET_FLOAT_WORD(b,hb);
  	}
***************
*** 54,63 ****
  		b *= t1;
  		a *= t1;
  		k -= 126;
! 	    } else {		/* scale a and b by 2^60 */
! 	        ha += 0x5d800000; 	/* a *= 2^60 */
! 		hb += 0x5d800000;	/* b *= 2^60 */
! 		k -= 60;
  		SET_FLOAT_WORD(a,ha);
  		SET_FLOAT_WORD(b,hb);
  	    }
--- 54,63 ----
  		b *= t1;
  		a *= t1;
  		k -= 126;
! 	    } else {		/* scale a and b by 2^80 */
! 	        ha += 0x28000000; 	/* a *= 2^80 */
! 		hb += 0x28000000;	/* b *= 2^80 */
! 		k -= 80;
  		SET_FLOAT_WORD(a,ha);
  		SET_FLOAT_WORD(b,hb);
  	    }

-------------- next part --------------
An embedded and charset-unspecified text was scrubbed...
Name: hypot-test.c
URL: <http://sourceware.org/pipermail/newlib/attachments/20020220/8d0d04a1/attachment.c>