[PATCH, LIBM] powf is wrong with x near 1 and |y| >> 1

[Steve Kargl]

See the last few comments here:

https://github.com/JuliaMath/openlibm/issues/211

Test program:

#include <stdio.h>
#include <math.h>
  
#if defined(__i386__)
#include <ieeefp.h>
#endif

int main()
{
  float x, y, z, a;
#if defined(__i386__)
        fpsetprec(FP_PE);
#endif
  x = 0x1.ffffeep-1f;
  y = -0x1.000002p+27f;
  z = powf (x, y);
  printf ("x=%e y=%e z=%e\n", x, y, z);
  printf ("x=%a y=%a z=%a\n", x, y, z);

  a = expf(y * logf(x)); // 0x1.d53532p+103f;
  printf ("floats: a=%a a=%e\n", a, a);
  a = (float)exp(y * log(x)); // 0x1.d53532p+103f;
  printf (" dbles: a=%a a=%e\n", a, a);
  return 0;
}

Patch without too much analysis for x near 1 or |y| > 2**27.

Index: src/e_powf.c
===================================================================
--- src/e_powf.c	(revision 2342)
+++ src/e_powf.c	(working copy)
@@ -136,7 +136,7 @@ __ieee754_powf(float x, float y)
     /* |y| is huge */
 	if(iy>0x4d000000) { /* if |y| > 2**27 */
 	/* over/underflow if x is not close to one */
-	    if(ix<0x3f7ffff8) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
+	    if(ix<0x3f7ffff7) return (hy<0)? sn*huge*huge:sn*tiny*tiny;
 	    if(ix>0x3f800007) return (hy>0)? sn*huge*huge:sn*tiny*tiny;
 	/* now |1-x| is tiny <= 2**-20, suffice to compute
 	   log(x) by x-x^2/2+x^3/3-x^4/4 */

-- 
Steve