lib/libm/pow.c

   1 /*                                                      pow.c
   2  *
   3  *      Power function
   4  *
   5  *
   6  *
   7  * SYNOPSIS:
   8  *
   9  * double x, y, z, pow();
  10  *
  11  * z = pow( x, y );
  12  *
  13  *
  14  *
  15  * DESCRIPTION:
  16  *
  17  * Computes x raised to the yth power.  Analytically,
  18  *
  19  *      x**y  =  exp( y log(x) ).
  20  *
  21  * Following Cody and Waite, this program uses a lookup table
  22  * of 2**-i/16 and pseudo extended precision arithmetic to
  23  * obtain an extra three bits of accuracy in both the logarithm
  24  * and the exponential.
  25  *
  26  *
  27  *
  28  * ACCURACY:
  29  *
  30  *                      Relative error:
  31  * arithmetic   domain     # trials      peak         rms
  32  *    IEEE     -26,26       30000      4.2e-16      7.7e-17
  33  *    DEC      -26,26       60000      4.8e-17      9.1e-18
  34  * 1/26 < x < 26, with log(x) uniformly distributed.
  35  * -26 < y < 26, y uniformly distributed.
  36  *    IEEE     0,8700       30000      1.5e-14      2.1e-15
  37  * 0.99 < x < 1.01, 0 < y < 8700, uniformly distributed.
  38  *
  39  *
  40  * ERROR MESSAGES:
  41  *
  42  *   message         condition      value returned
  43  * pow overflow     x**y > MAXNUM      INFINITY
  44  * pow underflow   x**y < 1/MAXNUM       0.0
  45  * pow domain      x<0 and y noninteger  0.0
  46  *
  47  */
  48 \f
  49 /*
  50 Cephes Math Library Release 2.8:  June, 2000
  51 Copyright 1984, 1995, 2000 by Stephen L. Moshier
  52 */
  53
  54
  55 #include "mconf.h"
  56 static char fname[] = {"pow"};
  57
  58 #define SQRTH 0.70710678118654752440
  59
  60 #ifdef UNK
  61 static double P[] = {
  62   4.97778295871696322025E-1,
  63   3.73336776063286838734E0,
  64   7.69994162726912503298E0,
  65   4.66651806774358464979E0
  66 };
  67 static double Q[] = {
  68 /* 1.00000000000000000000E0, */
  69   9.33340916416696166113E0,
  70   2.79999886606328401649E1,
  71   3.35994905342304405431E1,
  72   1.39995542032307539578E1
  73 };
  74 /* 2^(-i/16), IEEE precision */
  75 static double A[] = {
  76   1.00000000000000000000E0,
  77   9.57603280698573700036E-1,
  78   9.17004043204671215328E-1,
  79   8.78126080186649726755E-1,
  80   8.40896415253714502036E-1,
  81   8.05245165974627141736E-1,
  82   7.71105412703970372057E-1,
  83   7.38413072969749673113E-1,
  84   7.07106781186547572737E-1,
  85   6.77127773468446325644E-1,
  86   6.48419777325504820276E-1,
  87   6.20928906036742001007E-1,
  88   5.94603557501360513449E-1,
  89   5.69394317378345782288E-1,
  90   5.45253866332628844837E-1,
  91   5.22136891213706877402E-1,
  92   5.00000000000000000000E-1
  93 };
  94 static double B[] = {
  95  0.00000000000000000000E0,
  96  1.64155361212281360176E-17,
  97  4.09950501029074826006E-17,
  98  3.97491740484881042808E-17,
  99 -4.83364665672645672553E-17,
 100  1.26912513974441574796E-17,
 101  1.99100761573282305549E-17,
 102 -1.52339103990623557348E-17,
 103  0.00000000000000000000E0
 104 };
 105 static double R[] = {
 106  1.49664108433729301083E-5,
 107  1.54010762792771901396E-4,
 108  1.33335476964097721140E-3,
 109  9.61812908476554225149E-3,
 110  5.55041086645832347466E-2,
 111  2.40226506959099779976E-1,
 112  6.93147180559945308821E-1
 113 };
 114
 115 #define douba(k) A[k]
 116 #define doubb(k) B[k]
 117 #define MEXP 16383.0
 118 #ifdef DENORMAL
 119 #define MNEXP -17183.0
 120 #else
 121 #define MNEXP -16383.0
 122 #endif
 123 #endif
 124
 125 #ifdef DEC
 126 static unsigned short P[] = {
 127 0037776,0156313,0175332,0163602,
 128 0040556,0167577,0052366,0174245,
 129 0040766,0062753,0175707,0055564,
 130 0040625,0052035,0131344,0155636,
 131 };
 132 static unsigned short Q[] = {
 133 /*0040200,0000000,0000000,0000000,*/
 134 0041025,0052644,0154404,0105155,
 135 0041337,0177772,0007016,0047646,
 136 0041406,0062740,0154273,0020020,
 137 0041137,0177054,0106127,0044555,
 138 };
 139 static unsigned short A[] = {
 140 0040200,0000000,0000000,0000000,
 141 0040165,0022575,0012444,0103314,
 142 0040152,0140306,0163735,0022071,
 143 0040140,0146336,0166052,0112341,
 144 0040127,0042374,0145326,0116553,
 145 0040116,0022214,0012437,0102201,
 146 0040105,0063452,0010525,0003333,
 147 0040075,0004243,0117530,0006067,
 148 0040065,0002363,0031771,0157145,
 149 0040055,0054076,0165102,0120513,
 150 0040045,0177326,0124661,0050471,
 151 0040036,0172462,0060221,0120422,
 152 0040030,0033760,0050615,0134251,
 153 0040021,0141723,0071653,0010703,
 154 0040013,0112701,0161752,0105727,
 155 0040005,0125303,0063714,0044173,
 156 0040000,0000000,0000000,0000000
 157 };
 158 static unsigned short B[] = {
 159 0000000,0000000,0000000,0000000,
 160 0021473,0040265,0153315,0140671,
 161 0121074,0062627,0042146,0176454,
 162 0121413,0003524,0136332,0066212,
 163 0121767,0046404,0166231,0012553,
 164 0121257,0015024,0002357,0043574,
 165 0021736,0106532,0043060,0056206,
 166 0121310,0020334,0165705,0035326,
 167 0000000,0000000,0000000,0000000
 168 };
 169
 170 static unsigned short R[] = {
 171 0034173,0014076,0137624,0115771,
 172 0035041,0076763,0003744,0111311,
 173 0035656,0141766,0041127,0074351,
 174 0036435,0112533,0073611,0116664,
 175 0037143,0054106,0134040,0152223,
 176 0037565,0176757,0176026,0025551,
 177 0040061,0071027,0173721,0147572
 178 };
 179
 180 /*
 181 static double R[] = {
 182 0.14928852680595608186e-4,
 183 0.15400290440989764601e-3,
 184 0.13333541313585784703e-2,
 185 0.96181290595172416964e-2,
 186 0.55504108664085595326e-1,
 187 0.24022650695909537056e0,
 188 0.69314718055994529629e0
 189 };
 190 */
 191 #define douba(k) (*(double *)&A[(k)<<2])
 192 #define doubb(k) (*(double *)&B[(k)<<2])
 193 #define MEXP 2031.0
 194 #define MNEXP -2031.0
 195 #endif
 196
 197 #ifdef IBMPC
 198 static unsigned short P[] = {
 199 0x5cf0,0x7f5b,0xdb99,0x3fdf,
 200 0xdf15,0xea9e,0xddef,0x400d,
 201 0xeb6f,0x7f78,0xccbd,0x401e,
 202 0x9b74,0xb65c,0xaa83,0x4012,
 203 };
 204 static unsigned short Q[] = {
 205 /*0x0000,0x0000,0x0000,0x3ff0,*/
 206 0x914e,0x9b20,0xaab4,0x4022,
 207 0xc9f5,0x41c1,0xffff,0x403b,
 208 0x6402,0x1b17,0xccbc,0x4040,
 209 0xe92e,0x918a,0xffc5,0x402b,
 210 };
 211 static unsigned short A[] = {
 212 0x0000,0x0000,0x0000,0x3ff0,
 213 0x90da,0xa2a4,0xa4af,0x3fee,
 214 0xa487,0xdcfb,0x5818,0x3fed,
 215 0x529c,0xdd85,0x199b,0x3fec,
 216 0xd3ad,0x995a,0xe89f,0x3fea,
 217 0xf090,0x82a3,0xc491,0x3fe9,
 218 0xa0db,0x422a,0xace5,0x3fe8,
 219 0x0187,0x73eb,0xa114,0x3fe7,
 220 0x3bcd,0x667f,0xa09e,0x3fe6,
 221 0x5429,0xdd48,0xab07,0x3fe5,
 222 0x2a27,0xd536,0xbfda,0x3fe4,
 223 0x3422,0x4c12,0xdea6,0x3fe3,
 224 0xb715,0x0a31,0x06fe,0x3fe3,
 225 0x6238,0x6e75,0x387a,0x3fe2,
 226 0x517b,0x3c7d,0x72b8,0x3fe1,
 227 0x890f,0x6cf9,0xb558,0x3fe0,
 228 0x0000,0x0000,0x0000,0x3fe0
 229 };
 230 static unsigned short B[] = {
 231 0x0000,0x0000,0x0000,0x0000,
 232 0x3707,0xd75b,0xed02,0x3c72,
 233 0xcc81,0x345d,0xa1cd,0x3c87,
 234 0x4b27,0x5686,0xe9f1,0x3c86,
 235 0x6456,0x13b2,0xdd34,0xbc8b,
 236 0x42e2,0xafec,0x4397,0x3c6d,
 237 0x82e4,0xd231,0xf46a,0x3c76,
 238 0x8a76,0xb9d7,0x9041,0xbc71,
 239 0x0000,0x0000,0x0000,0x0000
 240 };
 241 static unsigned short R[] = {
 242 0x937f,0xd7f2,0x6307,0x3eef,
 243 0x9259,0x60fc,0x2fbe,0x3f24,
 244 0xef1d,0xc84a,0xd87e,0x3f55,
 245 0x33b7,0x6ef1,0xb2ab,0x3f83,
 246 0x1a92,0xd704,0x6b08,0x3fac,
 247 0xc56d,0xff82,0xbfbd,0x3fce,
 248 0x39ef,0xfefa,0x2e42,0x3fe6
 249 };
 250
 251 #define douba(k) (*(double *)&A[(k)<<2])
 252 #define doubb(k) (*(double *)&B[(k)<<2])
 253 #define MEXP 16383.0
 254 #ifdef DENORMAL
 255 #define MNEXP -17183.0
 256 #else
 257 #define MNEXP -16383.0
 258 #endif
 259 #endif
 260
 261 #ifdef MIEEE
 262 static unsigned short P[] = {
 263 0x3fdf,0xdb99,0x7f5b,0x5cf0,
 264 0x400d,0xddef,0xea9e,0xdf15,
 265 0x401e,0xccbd,0x7f78,0xeb6f,
 266 0x4012,0xaa83,0xb65c,0x9b74
 267 };
 268 static unsigned short Q[] = {
 269 0x4022,0xaab4,0x9b20,0x914e,
 270 0x403b,0xffff,0x41c1,0xc9f5,
 271 0x4040,0xccbc,0x1b17,0x6402,
 272 0x402b,0xffc5,0x918a,0xe92e
 273 };
 274 static unsigned short A[] = {
 275 0x3ff0,0x0000,0x0000,0x0000,
 276 0x3fee,0xa4af,0xa2a4,0x90da,
 277 0x3fed,0x5818,0xdcfb,0xa487,
 278 0x3fec,0x199b,0xdd85,0x529c,
 279 0x3fea,0xe89f,0x995a,0xd3ad,
 280 0x3fe9,0xc491,0x82a3,0xf090,
 281 0x3fe8,0xace5,0x422a,0xa0db,
 282 0x3fe7,0xa114,0x73eb,0x0187,
 283 0x3fe6,0xa09e,0x667f,0x3bcd,
 284 0x3fe5,0xab07,0xdd48,0x5429,
 285 0x3fe4,0xbfda,0xd536,0x2a27,
 286 0x3fe3,0xdea6,0x4c12,0x3422,
 287 0x3fe3,0x06fe,0x0a31,0xb715,
 288 0x3fe2,0x387a,0x6e75,0x6238,
 289 0x3fe1,0x72b8,0x3c7d,0x517b,
 290 0x3fe0,0xb558,0x6cf9,0x890f,
 291 0x3fe0,0x0000,0x0000,0x0000
 292 };
 293 static unsigned short B[] = {
 294 0x0000,0x0000,0x0000,0x0000,
 295 0x3c72,0xed02,0xd75b,0x3707,
 296 0x3c87,0xa1cd,0x345d,0xcc81,
 297 0x3c86,0xe9f1,0x5686,0x4b27,
 298 0xbc8b,0xdd34,0x13b2,0x6456,
 299 0x3c6d,0x4397,0xafec,0x42e2,
 300 0x3c76,0xf46a,0xd231,0x82e4,
 301 0xbc71,0x9041,0xb9d7,0x8a76,
 302 0x0000,0x0000,0x0000,0x0000
 303 };
 304 static unsigned short R[] = {
 305 0x3eef,0x6307,0xd7f2,0x937f,
 306 0x3f24,0x2fbe,0x60fc,0x9259,
 307 0x3f55,0xd87e,0xc84a,0xef1d,
 308 0x3f83,0xb2ab,0x6ef1,0x33b7,
 309 0x3fac,0x6b08,0xd704,0x1a92,
 310 0x3fce,0xbfbd,0xff82,0xc56d,
 311 0x3fe6,0x2e42,0xfefa,0x39ef
 312 };
 313
 314 #define douba(k) (*(double *)&A[(k)<<2])
 315 #define doubb(k) (*(double *)&B[(k)<<2])
 316 #define MEXP 16383.0
 317 #ifdef DENORMAL
 318 #define MNEXP -17183.0
 319 #else
 320 #define MNEXP -16383.0
 321 #endif
 322 #endif
 323
 324 /* log2(e) - 1 */
 325 #define LOG2EA 0.44269504088896340736
 326
 327 #define F W
 328 #define Fa Wa
 329 #define Fb Wb
 330 #define G W
 331 #define Ga Wa
 332 #define Gb u
 333 #define H W
 334 #define Ha Wb
 335 #define Hb Wb
 336
 337 #ifdef ANSIPROT
 338 extern double floor ( double );
 339 extern double fabs ( double );
 340 extern double frexp ( double, int * );
 341 extern double ldexp ( double, int );
 342 extern double polevl ( double, void *, int );
 343 extern double p1evl ( double, void *, int );
 344 extern double powi ( double, int );
 345 extern int signbit ( double );
 346 extern int isnan ( double );
 347 extern int isfinite ( double );
 348 static double reduc ( double );
 349 #else
 350 double floor(), fabs(), frexp(), ldexp();
 351 double polevl(), p1evl(), powi();
 352 int signbit(), isnan(), isfinite();
 353 static double reduc();
 354 #endif
 355 extern double MAXNUM;
 356 extern double INFINITY;
 357 extern double NAN;
 358 extern double NEGZERO;
 359
 360 double pow( x, y )
 361 double x, y;
 362 {
 363 double w, z, W, Wa, Wb, ya, yb, u;
 364 /* double F, Fa, Fb, G, Ga, Gb, H, Ha, Hb */
 365 double aw, ay, wy;
 366 int e, i, nflg, iyflg, yoddint;
 367
 368 if( y == 0.0 )
 369         return( 1.0 );
 370 #ifdef NANS
 371 if( isnan(x) )
 372         return( x );
 373 if( isnan(y) )
 374         return( y );
 375 #endif
 376 if( y == 1.0 )
 377         return( x );
 378
 379
 380 #ifdef INFINITIES
 381 if( !isfinite(y) && (x == 1.0 || x == -1.0) )
 382         {
 383         mtherr( "pow", DOMAIN );
 384 #ifdef NANS
 385         return( NAN );
 386 #else
 387         return( INFINITY );
 388 #endif
 389         }
 390 #endif
 391
 392 if( x == 1.0 )
 393         return( 1.0 );
 394
 395 if( y >= MAXNUM )
 396         {
 397 #ifdef INFINITIES
 398         if( x > 1.0 )
 399                 return( INFINITY );
 400 #else
 401         if( x > 1.0 )
 402                 return( MAXNUM );
 403 #endif
 404         if( x > 0.0 && x < 1.0 )
 405                 return( 0.0);
 406         if( x < -1.0 )
 407                 {
 408 #ifdef INFINITIES
 409                 return( INFINITY );
 410 #else
 411                 return( MAXNUM );
 412 #endif
 413                 }
 414         if( x > -1.0 && x < 0.0 )
 415                 return( 0.0 );
 416         }
 417 if( y <= -MAXNUM )
 418         {
 419         if( x > 1.0 )
 420                 return( 0.0 );
 421 #ifdef INFINITIES
 422         if( x > 0.0 && x < 1.0 )
 423                 return( INFINITY );
 424 #else
 425         if( x > 0.0 && x < 1.0 )
 426                 return( MAXNUM );
 427 #endif
 428         if( x < -1.0 )
 429                 return( 0.0 );
 430 #ifdef INFINITIES
 431         if( x > -1.0 && x < 0.0 )
 432                 return( INFINITY );
 433 #else
 434         if( x > -1.0 && x < 0.0 )
 435                 return( MAXNUM );
 436 #endif
 437         }
 438 if( x >= MAXNUM )
 439         {
 440 #if INFINITIES
 441         if( y > 0.0 )
 442                 return( INFINITY );
 443 #else
 444         if( y > 0.0 )
 445                 return( MAXNUM );
 446 #endif
 447         return(0.0);
 448         }
 449 /* Set iyflg to 1 if y is an integer.  */
 450 iyflg = 0;
 451 w = floor(y);
 452 if( w == y )
 453         iyflg = 1;
 454
 455 /* Test for odd integer y.  */
 456 yoddint = 0;
 457 if( iyflg )
 458         {
 459         ya = fabs(y);
 460         ya = floor(0.5 * ya);
 461         yb = 0.5 * fabs(w);
 462         if( ya != yb )
 463                 yoddint = 1;
 464         }
 465
 466 if( x <= -MAXNUM )
 467         {
 468         if( y > 0.0 )
 469                 {
 470 #ifdef INFINITIES
 471                 if( yoddint )
 472                         return( -INFINITY );
 473                 return( INFINITY );
 474 #else
 475                 if( yoddint )
 476                         return( -MAXNUM );
 477                 return( MAXNUM );
 478 #endif
 479                 }
 480         if( y < 0.0 )
 481                 {
 482 #ifdef MINUSZERO
 483                 if( yoddint )
 484                         return( NEGZERO );
 485 #endif
 486                 return( 0.0 );
 487                 }
 488         }
 489
 490 nflg = 0;       /* flag = 1 if x<0 raised to integer power */
 491 if( x <= 0.0 )
 492         {
 493         if( x == 0.0 )
 494                 {
 495                 if( y < 0.0 )
 496                         {
 497 #ifdef MINUSZERO
 498                         if( signbit(x) && yoddint )
 499                                 return( -INFINITY );
 500 #endif
 501 #ifdef INFINITIES
 502                         return( INFINITY );
 503 #else
 504                         return( MAXNUM );
 505 #endif
 506                         }
 507                 if( y > 0.0 )
 508                         {
 509 #ifdef MINUSZERO
 510                         if( signbit(x) && yoddint )
 511                                 return( NEGZERO );
 512 #endif
 513                         return( 0.0 );
 514                         }
 515                 return( 1.0 );
 516                 }
 517         else
 518                 {
 519                 if( iyflg == 0 )
 520                         { /* noninteger power of negative number */
 521                         mtherr( fname, DOMAIN );
 522 #ifdef NANS
 523                         return(NAN);
 524 #else
 525                         return(0.0);
 526 #endif
 527                         }
 528                 nflg = 1;
 529                 }
 530         }
 531
 532 /* Integer power of an integer.  */
 533
 534 if( iyflg )
 535         {
 536         i = w;
 537         w = floor(x);
 538         if( (w == x) && (fabs(y) < 32768.0) )
 539                 {
 540                 w = powi( x, (int) y );
 541                 return( w );
 542                 }
 543         }
 544
 545 if( nflg )
 546         x = fabs(x);
 547
 548 /* For results close to 1, use a series expansion.  */
 549 w = x - 1.0;
 550 aw = fabs(w);
 551 ay = fabs(y);
 552 wy = w * y;
 553 ya = fabs(wy);
 554 if((aw <= 1.0e-3 && ay <= 1.0)
 555    || (ya <= 1.0e-3 && ay >= 1.0))
 556         {
 557         z = (((((w*(y-5.)/720. + 1./120.)*w*(y-4.) + 1./24.)*w*(y-3.)
 558                 + 1./6.)*w*(y-2.) + 0.5)*w*(y-1.) )*wy + wy + 1.;
 559         goto done;
 560         }
 561 /* These are probably too much trouble.  */
 562 #if 0
 563 w = y * log(x);
 564 if (aw > 1.0e-3 && fabs(w) < 1.0e-3)
 565   {
 566     z = ((((((
 567     w/7. + 1.)*w/6. + 1.)*w/5. + 1.)*w/4. + 1.)*w/3. + 1.)*w/2. + 1.)*w + 1.;
 568     goto done;
 569   }
 570
 571 if(ya <= 1.0e-3 && aw <= 1.0e-4)
 572   {
 573     z = (((((
 574              wy*1./720.
 575              + (-w*1./48. + 1./120.) )*wy
 576             + ((w*17./144. - 1./12.)*w + 1./24.) )*wy
 577            + (((-w*5./16. + 7./24.)*w - 1./4.)*w + 1./6.) )*wy
 578           + ((((w*137./360. - 5./12.)*w + 11./24.)*w - 1./2.)*w + 1./2.) )*wy
 579          + (((((-w*1./6. + 1./5.)*w - 1./4)*w + 1./3.)*w -1./2.)*w ) )*wy
 580            + wy + 1.0;
 581     goto done;
 582   }
 583 #endif
 584
 585 /* separate significand from exponent */
 586 x = frexp( x, &e );
 587
 588 #if 0
 589 /* For debugging, check for gross overflow. */
 590 if( (e * y)  > (MEXP + 1024) )
 591         goto overflow;
 592 #endif
 593
 594 /* Find significand of x in antilog table A[]. */
 595 i = 1;
 596 if( x <= douba(9) )
 597         i = 9;
 598 if( x <= douba(i+4) )
 599         i += 4;
 600 if( x <= douba(i+2) )
 601         i += 2;
 602 if( x >= douba(1) )
 603         i = -1;
 604 i += 1;
 605
 606
 607 /* Find (x - A[i])/A[i]
 608  * in order to compute log(x/A[i]):
 609  *
 610  * log(x) = log( a x/a ) = log(a) + log(x/a)
 611  *
 612  * log(x/a) = log(1+v),  v = x/a - 1 = (x-a)/a
 613  */
 614 x -= douba(i);
 615 x -= doubb(i/2);
 616 x /= douba(i);
 617
 618
 619 /* rational approximation for log(1+v):
 620  *
 621  * log(1+v)  =  v  -  v**2/2  +  v**3 P(v) / Q(v)
 622  */
 623 z = x*x;
 624 w = x * ( z * polevl( x, P, 3 ) / p1evl( x, Q, 4 ) );
 625 w = w - ldexp( z, -1 );   /*  w - 0.5 * z  */
 626
 627 /* Convert to base 2 logarithm:
 628  * multiply by log2(e)
 629  */
 630 w = w + LOG2EA * w;
 631 /* Note x was not yet added in
 632  * to above rational approximation,
 633  * so do it now, while multiplying
 634  * by log2(e).
 635  */
 636 z = w + LOG2EA * x;
 637 z = z + x;
 638
 639 /* Compute exponent term of the base 2 logarithm. */
 640 w = -i;
 641 w = ldexp( w, -4 );     /* divide by 16 */
 642 w += e;
 643 /* Now base 2 log of x is w + z. */
 644
 645 /* Multiply base 2 log by y, in extended precision. */
 646
 647 /* separate y into large part ya
 648  * and small part yb less than 1/16
 649  */
 650 ya = reduc(y);
 651 yb = y - ya;
 652
 653
 654 F = z * y  +  w * yb;
 655 Fa = reduc(F);
 656 Fb = F - Fa;
 657
 658 G = Fa + w * ya;
 659 Ga = reduc(G);
 660 Gb = G - Ga;
 661
 662 H = Fb + Gb;
 663 Ha = reduc(H);
 664 w = ldexp( Ga+Ha, 4 );
 665
 666 /* Test the power of 2 for overflow */
 667 if( w > MEXP )
 668         {
 669 #ifndef INFINITIES
 670         mtherr( fname, OVERFLOW );
 671 #endif
 672 #ifdef INFINITIES
 673         if( nflg && yoddint )
 674           return( -INFINITY );
 675         return( INFINITY );
 676 #else
 677         if( nflg && yoddint )
 678           return( -MAXNUM );
 679         return( MAXNUM );
 680 #endif
 681         }
 682
 683 if( w < (MNEXP - 1) )
 684         {
 685 #ifndef DENORMAL
 686         mtherr( fname, UNDERFLOW );
 687 #endif
 688 #ifdef MINUSZERO
 689         if( nflg && yoddint )
 690           return( NEGZERO );
 691 #endif
 692         return( 0.0 );
 693         }
 694
 695 e = w;
 696 Hb = H - Ha;
 697
 698 if( Hb > 0.0 )
 699         {
 700         e += 1;
 701         Hb -= 0.0625;
 702         }
 703
 704 /* Now the product y * log2(x)  =  Hb + e/16.0.
 705  *
 706  * Compute base 2 exponential of Hb,
 707  * where -0.0625 <= Hb <= 0.
 708  */
 709 z = Hb * polevl( Hb, R, 6 );  /*    z  =  2**Hb - 1    */
 710
 711 /* Express e/16 as an integer plus a negative number of 16ths.
 712  * Find lookup table entry for the fractional power of 2.
 713  */
 714 if( e < 0 )
 715         i = 0;
 716 else
 717         i = 1;
 718 i = e/16 + i;
 719 e = 16*i - e;
 720 w = douba( e );
 721 z = w + w * z;      /*    2**-e * ( 1 + (2**Hb-1) )    */
 722 z = ldexp( z, i );  /* multiply by integer power of 2 */
 723
 724 done:
 725
 726 /* Negate if odd integer power of negative number */
 727 if( nflg && yoddint )
 728         {
 729 #ifdef MINUSZERO
 730         if( z == 0.0 )
 731                 z = NEGZERO;
 732         else
 733 #endif
 734                 z = -z;
 735         }
 736 return( z );
 737 }
 738
 739
 740 /* Find a multiple of 1/16 that is within 1/16 of x. */
 741 static double reduc(x)
 742 double x;
 743 {
 744 double t;
 745
 746 t = ldexp( x, 4 );
 747 t = floor( t );
 748 t = ldexp( t, -4 );
 749 return(t);
 750 }