23  *  P[0] x^n  +  P[1] x^(n-1)  +  ...  +  P[n]
42  *  x^n  +  P[0] x^(n-1)  +  P[1] x^(n-2)  +  ...  +  P[n]
51 	n -= 1;  in __p1evll()
92  * of the approximation interval the logarithm equal 1/32
93  * and its relative error is about 1 lsb = 1.1e-19.  Hence
94  * the predicted relative error in the result is 2.3e-21 y .
99  *    IEEE     +-1000       40000      2.8e-18      3.7e-19
103  *    IEEE     0,8700       60000      6.5e-18      1.0e-18
111  * pow underflow   x**y < 1/MAXNUM       0.0
126 /* log(1+x) =  x - .5x^2 + x^3 *  P(z)/Q(z)
127  * on the domain  2^(-1/32) - 1  <=  x  <=  2^(1/32) - 1
130  8.3319510773868690346226E-4L,
131  4.9000050881978028599627E-1L,
146  9.7857206208770013448287E-1L,
147  9.5760328069857364691013E-1L,
148  9.3708381705514995065011E-1L,
149  9.1700404320467123175367E-1L,
150  8.9735453750155359320742E-1L,
151  8.7812608018664974155474E-1L,
152  8.5930964906123895780165E-1L,
153  8.4089641525371454301892E-1L,
154  8.2287773907698242225554E-1L,
155  8.0524516597462715409607E-1L,
156  7.8799042255394324325455E-1L,
157  7.7110541270397041179298E-1L,
158  7.5458221379671136985669E-1L,
159  7.3841307296974965571198E-1L,
160  7.2259040348852331001267E-1L,
161  7.0710678118654752438189E-1L,
162  6.9195494098191597746178E-1L,
163  6.7712777346844636413344E-1L,
164  6.6261832157987064729696E-1L,
165  6.4841977732550483296079E-1L,
166  6.3452547859586661129850E-1L,
167  6.2092890603674202431705E-1L,
168  6.0762367999023443907803E-1L,
169  5.9460355750136053334378E-1L,
170  5.8186242938878875689693E-1L,
171  5.6939431737834582684856E-1L,
172  5.5719337129794626814472E-1L,
173  5.4525386633262882960438E-1L,
174  5.3357020033841180906486E-1L,
175  5.2213689121370692017331E-1L,
176  5.1094857432705833910408E-1L,
177  5.0000000000000000000000E-1L,
181  2.6176170809902549338711E-20L,
182 -1.0126791927256478897086E-20L,
183  1.3438228172316276937655E-21L,
184  1.2207982955417546912101E-20L,
185 -6.3084814358060867200133E-21L,
186  1.3164426894366316434230E-20L,
187 -1.8527916071632873716786E-20L,
188  1.8950325588932570796551E-20L,
189  1.5564775779538780478155E-20L,
190  6.0859793637556860974380E-21L,
191 -2.0208749253662532228949E-20L,
192  1.4966292219224761844552E-20L,
193  3.3540909728056476875639E-21L,
194 -8.6987564101742849540743E-22L,
195 -1.2327176863327626135542E-20L,
199 /* 2^x = 1 + x P(x),
200  * on the interval -1/32 <= x <= 0
203  1.5089970579127659901157E-5L,
204  1.5402715328927013076125E-4L,
205  1.3333556028915671091390E-3L,
206  9.6181291046036762031786E-3L,
207  5.5504108664798463044015E-2L,
208  2.4022650695910062854352E-1L,
209  6.9314718055994530931447E-1L,
217 /* log2(e) - 1 */
232 static const long double LOGE2L = 6.9314718055994530941723E-1L;
250 long e;  in powl()  local
299 /* Set iyflg to 1 if y is an integer.  */  in powl()
302 	iyflg = 1;  in powl()
312 		yoddint = 1;  in powl()
332 nflg = 0;	/* flag = 1 if x<0 raised to integer power */  in powl()
358 		nflg = 1;  in powl()
381 e = i;  in powl()
384 i = 1;  in powl()
393 if( x >= douba(1) )  in powl()
394 	i = -1;  in powl()
395 i += 1;  in powl()
403  * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a  in powl()
410 /* rational approximation for log(1+v):  in powl()
412  * log(1+v)  =  v  -  v**2/2  +  v**3 P(v) / Q(v)  in powl()
416 w = w - ldexpl( z, -1 );   /*  w - 0.5 * z  */  in powl()
419  * multiply by log2(e) = 1 + LOG2EA  in powl()
429 w += e;  in powl()
435  * and small part yb less than 1/NXT  in powl()
462 e = w;  in powl()
467 	e += 1;  in powl()
471 /* Now the product y * log2(x)  =  Hb + e/NXT.  in powl()
476 z = Hb * __polevll( Hb, R, 6 );  /*    z  =  2**Hb - 1    */  in powl()
478 /* Express e/NXT as an integer plus a negative number of (1/NXT)ths.  in powl()
481 if( e < 0 )  in powl()
484 	i = 1;  in powl()
485 i = e/NXT + i;  in powl()
486 e = NXT*i - e;  in powl()
487 w = douba( e );  in powl()
488 z = w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */  in powl()
498 	w = ldexpl( y, -1 );  in powl()
500 	w = ldexpl( w, 1 );  in powl()
509 /* Find a multiple of 1/NXT that is within 1/NXT of x. */
551  *    IEEE     .001,1000  -1022,1023  50000       4.3e-17     7.8e-18
552  *    IEEE        1,2     -1022,1023  20000       3.9e-17     7.6e-18
553  *    IEEE     .99,1.01     0,8700    10000       3.6e-16     7.2e-17
564 int n, e, sign, asign, lx;  in powil()  local
582 	asign = -1;  in powil()
591 	sign = -1;  in powil()
596 	sign = 1;  in powil()
605 e = (lx - 1)*n;  in powil()
606 if( (e == 0) || (e > 64) || (e < -64) )  in powil()
608 	s = (s - 7.0710678118654752e-1L) / (s +  7.0710678118654752e-1L);  in powil()
613 	s = LOGE2L * e;  in powil()
632 if( n & 1 )  in powil()
642 n >>= 1;  in powil()
646 	if( n & 1 )	/* if that bit is set, then include in product */  in powil()
648 	n >>= 1;  in powil()