27  * long double x, y, z, powl();
29  * z = powl( x, y );
37  *      x**y  =  exp( y log(x) ).
48  * The relative error of pow(x,y) can be estimated
49  * by   y dl ln(2),   where dl is the absolute error of
53  * the predicted relative error in the result is 2.3e-21 y .
60  * -1000 < y < 1000, y uniformly distributed.
63  * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
69  * pow overflow     x**y > MAXNUM      INFINITY
70  * pow underflow   x**y < 1/MAXNUM       0.0
71  * pow domain      x<0 and y noninteger  0.0
186 powl(long double x, long double y)  in powl()  argument
193 if( !isnanl_inline(y) && y == 0.0L ) {  in powl()
195                 return x + y;  in powl()
200         if (issignalingl_inline(y))  in powl()
201                 return x + y;  in powl()
206 	return( x + y );  in powl()
207 if( isnanl_inline(y) )  in powl()
208 	return( x + y );  in powl()
210 w = floorl(y);  in powl()
211 /* Set iyflg to 1 if y is an integer.  */  in powl()
212 iyflg = (w == y);  in powl()
217 /* Test for odd integer y and negative x (including negative zero)  */  in powl()
221         ya = ldexpl(y, -1);  in powl()
226         if( y < 0 )  in powl()
230 if( y == 1.0L )  in powl()
233 if( !isfinite(y) ) {  in powl()
236         if (y < 0) {  in powl()
247 /* y >= 2**80, infinity (x > 1) or zero (x < 1) */  in powl()
248 if( y >= 0x1p80L)  in powl()
259 /* y <= -2**80, zero (x > 1) or infinity (x < 1) */  in powl()
260 if( y <= -0x1p80L)  in powl()
273 	if( y > 0.0L )  in powl()
282 	if( y > 0.0L )  in powl()
288 	if( y < 0.0L )  in powl()
301 		if( y < 0.0L )  in powl()
307 		if( y > 0.0L )  in powl()
313 		if( y == 0.0L )  in powl()
316 			return( 0.0L );  /*   0**y   */  in powl()
330 	if( (w == x) && (fabsl(y) < 32768.0L) )  in powl()
332 		w = powil( x, (int) y );  in powl()
394 /* Multiply base 2 log by y, in extended precision. */  in powl()
396 /* separate y into large part ya  in powl()
399 ya = reducl(y);  in powl()
400 yb = y - ya;  in powl()
403  * = w*ya + w*yb + z*y  in powl()
405 F = z * y  +  w * yb;  in powl()
433 /* Now the product y * log2(x)  =  Hb + e/NXT.  in powl()
460 	w = ldexpl( y, -1 );  in powl()
463 	if( w != y )  in powl()
493  * long double x, y, powil();
496  * y = powil( x, n );
526 long double ww, y;  in powil()  local
598 	y = x;  in powil()
602 	y = 1.0L;  in powil()
612 		y *= ww;  in powil()
617 	y = -y; /* odd power of negative number */  in powil()
619 	y = 1.0L/y;  in powil()
620 return(y);  in powil()