Re: lib/50698: hypotf() of small numbers yields infinity

[Rin Okuyama <okuyama%flex.phys.tohoku.ac.jp@localhost>]

The following reply was made to PR lib/50698; it has been noted by GNATS.

From: Rin Okuyama <okuyama%flex.phys.tohoku.ac.jp@localhost>
To: gnats-bugs%NetBSD.org@localhost
Cc: 
Subject: Re: lib/50698: hypotf() of small numbers yields infinity
Date: Sun, 24 Jan 2016 19:26:06 +0900

 This is because some magic numbers in e_hypotf.c are wrong. In FreeBSD,
 this bug was fixed by commit #23397:
 
   https://svnweb.freebsd.org/base/head/lib/msun/src/e_hypotf.c#rev23397
 
 I attached the minimal patch below. With this patch, I obtain (on amd64)
 
   % ./test-hypotf 1e-18 1e-18
   1.41421363268178e-18
   %
 
 Other changes have been made for libm in FreeBSD. IMO, we should merge
 reasonable ones.
 
 --- src/lib/libm/src/e_hypotf.c.orig	2016-01-24 18:40:38.426417738 +0900
 +++ src/lib/libm/src/e_hypotf.c	2016-01-24 19:07:19.114100741 +0900
 @@ -43,22 +43,22 @@
  	       if(hb == 0x7f800000) w = b;
  	       return w;
  	   }
 -	   /* scale a and b by 2**-60 */
 -	   ha -= 0x5d800000; hb -= 0x5d800000;	k += 60;
 +	   /* scale a and b by 2**-68 */
 +	   ha -= 0x22000000; hb -= 0x22000000;	k += 68;
  	   SET_FLOAT_WORD(a,ha);
  	   SET_FLOAT_WORD(b,hb);
  	}
  	if(hb < 0x26800000) {	/* b < 2**-50 */
  	    if(hb <= 0x007fffff) {	/* subnormal b or 0 */
  	        if(hb==0) return a;
 -		SET_FLOAT_WORD(t1,0x3f000000);	/* t1=2^126 */
 +		SET_FLOAT_WORD(t1,0x7e800000);	/* t1=2^126 */
  		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;
 +	    } else {		/* scale a and b by 2^68 */
 +	        ha += 0x22000000; 	/* a *= 2^68 */
 +		hb += 0x22000000;	/* b *= 2^68 */
 +		k -= 68;
  		SET_FLOAT_WORD(a,ha);
  		SET_FLOAT_WORD(b,hb);
  	    }