RFA: patch to fix scaling in ef_hypot.c
J. Johnston
jjohnstn@redhat.com
Fri Feb 22 12:13:00 GMT 2002
Richard Sandiford wrote:
>
> 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?
>
Yes. Thanks.
-- Jeff J.
>
> * 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);
> }
>
> --------------------------------------------------------------------------------
> Name: hypot-test.c
> hypot-test.c Type: Plain Text (text/plain)
> Description: Test program
More information about the Newlib
mailing list