Re: lib/49240: Some corner cases for pow() fail to work correctly

[Havard Eidnes <he%NetBSD.org@localhost>]

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

From: Havard Eidnes <he%NetBSD.org@localhost>
To: gnats-bugs%NetBSD.org@localhost
Cc: 
Subject: Re: lib/49240: Some corner cases for pow() fail to work correctly
Date: Wed, 01 Oct 2014 12:26:54 +0200 (CEST)

 > I'll work out a minimal diff for our e_pow.c.
 
 This follows below.  With that and the test correction for C99
 behaviour, the tests now all pass when I test locally.
 
 The changes rolled into this are:
 
  * Change yy1 to y1; no need for yy1.
  * Add a case for 1**y = 1, even if y is NaN (or +-Inf)
  * Use (x+0.0)+(y+0.0) instead of x+y when mixing NaNs.  There's
    a long explanation of why this is done in FreeBSD's commit log
    for e_pow.c at
    http://svnweb.freebsd.org/base/head/lib/msun/src/e_pow.c?view=log
  * Add a case for -1**+-Inf, since C99 and IEEE 754 revised 2008(?)
    reportedly state that this should be 1
  * Add a fix to get correct sign for some over- and underflow
    cases.  This seems to come from fdlib 5.3, FreeBSD revision
    129956. (Ref. above SVN log.)
  * Add a fix to get well-defined behaviour when right-shifting a
    possibly negative value by adding a cast.
  * Fix comment detailing the expected result for the special cases
 
 
 --- /usr/src/lib/libm/src/e_pow.c	2013-03-27 19:12:07.000000000 +0100
 +++ e_pow.c	2014-10-01 11:41:10.000000000 +0200
 @@ -29,13 +29,13 @@
   * Special cases:
   *	1.  (anything) ** 0  is 1
   *	2.  (anything) ** 1  is itself
 - *	3.  (anything) ** NAN is NAN
 + *	3.  (anything) ** NAN is NAN except 1 ** NAN = 1
   *	4.  NAN ** (anything except 0) is NAN
   *	5.  +-(|x| > 1) **  +INF is +INF
   *	6.  +-(|x| > 1) **  -INF is +0
   *	7.  +-(|x| < 1) **  +INF is +0
   *	8.  +-(|x| < 1) **  -INF is +INF
 - *	9.  +-1         ** +-INF is NAN
 + *	9.  +-1         ** +-INF is 1
   *	10. +0 ** (+anything except 0, NAN)               is +0
   *	11. -0 ** (+anything except 0, NAN, odd integer)  is +0
   *	12. +0 ** (-anything except 0, NAN)               is +INF
 @@ -98,10 +98,10 @@
  ivln2_l  =  1.92596299112661746887e-08; /* 0x3E54AE0B, 0xF85DDF44 =1/ln2 tail*/
  
  double
 -__ieee754_pow(double x, double y)
 +pow(double x, double y)
  {
  	double z,ax,z_h,z_l,p_h,p_l;
 -	double yy1,t1,t2,r,s,t,u,v,w;
 +	double y1,t1,t2,r,s,t,u,v,w;
  	int32_t i,j,k,yisint,n;
  	int32_t hx,hy,ix,iy;
  	u_int32_t lx,ly;
 @@ -113,10 +113,13 @@
      /* y==zero: x**0 = 1 */
  	if((iy|ly)==0) return one;
  
 -    /* +-NaN return x+y */
 +    /* x==1: 1**y = 1, even if y is NaN */
 +	if (hx==0x3ff00000 && lx == 0) return one;
 +
 +    /* y!=zero: result is NaN if either arg is NaN */
  	if(ix > 0x7ff00000 || ((ix==0x7ff00000)&&(lx!=0)) ||
  	   iy > 0x7ff00000 || ((iy==0x7ff00000)&&(ly!=0)))
 -		return x+y;
 +		return (x+0.0)+(y+0.0);
  
      /* determine if y is an odd int when x < 0
       * yisint = 0	... y is not an integer
 @@ -142,7 +145,7 @@
  	if(ly==0) {
  	    if (iy==0x7ff00000) {	/* y is +-inf */
  	        if(((ix-0x3ff00000)|lx)==0)
 -		    return  y - y;	/* inf**+-1 is NaN */
 +		    return one;		/* (-1)**+-inf is 1 */
  	        else if (ix >= 0x3ff00000)/* (|x|>1)**+-inf = inf,0 */
  		    return (hy>=0)? y: zero;
  	        else			/* (|x|<1)**-,+inf = inf,0 */
 @@ -174,7 +177,11 @@
  	    }
  	}
  
 +    /* CYGNUS LOCAL + fdlibm-5.3 fix: This used to be
  	n = (hx>>31)+1;
 +	but ANSI C says a right shift of a signed negative quantity is
 +	implementation defined.  */
 +	n = ((u_int32_t)hx>>31)-1;
  
      /* (x<0)**(non-int) is NaN */
  	if((n|yisint)==0) return (x-x)/(x-x);
 @@ -250,11 +257,11 @@
  	    t2 = z_l-(((t1-t)-dp_h[k])-z_h);
  	}
  
 -    /* split up y into yy1+y2 and compute (yy1+y2)*(t1+t2) */
 -	yy1  = y;
 -	SET_LOW_WORD(yy1,0);
 -	p_l = (y-yy1)*t1+y*t2;
 -	p_h = yy1*t1;
 +    /* split up y into y1+y2 and compute (y1+y2)*(t1+t2) */
 +	y1  = y;
 +	SET_LOW_WORD(y1,0);
 +	p_l = (y-y1)*t1+y*t2;
 +	p_h = y1*t1;
  	z = p_l+p_h;
  	EXTRACT_WORDS(j,i,z);
  	if (j>=0x40900000) {				/* z >= 1024 */