Re: RFA: patch to fix scaling in ef_hypot.c

["J. Johnston" <jjohnstn at redhat dot com>]

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