StdLib/LibC/gdtoa/strtod.c

   1 /* $NetBSD: strtod.c,v 1.4.14.1 2008/04/08 21:10:55 jdc Exp $ */
   2
   3 /****************************************************************
   4
   5 The author of this software is David M. Gay.
   6
   7 Copyright (C) 1998-2001 by Lucent Technologies
   8 All Rights Reserved
   9
  10 Permission to use, copy, modify, and distribute this software and
  11 its documentation for any purpose and without fee is hereby
  12 granted, provided that the above copyright notice appear in all
  13 copies and that both that the copyright notice and this
  14 permission notice and warranty disclaimer appear in supporting
  15 documentation, and that the name of Lucent or any of its entities
  16 not be used in advertising or publicity pertaining to
  17 distribution of the software without specific, written prior
  18 permission.
  19
  20 LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE,
  21 INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS.
  22 IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY
  23 SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES
  24 WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER
  25 IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION,
  26 ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF
  27 THIS SOFTWARE.
  28
  29 ****************************************************************/
  30
  31 /* Please send bug reports to David M. Gay (dmg at acm dot org,
  32  * with " at " changed at "@" and " dot " changed to ".").  */
  33 #include  <LibConfig.h>
  34
  35 #include "gdtoaimp.h"
  36 #ifndef NO_FENV_H
  37 #include <fenv.h>
  38 #endif
  39
  40 #ifdef USE_LOCALE
  41 #include "locale.h"
  42 #endif
  43
  44 #ifdef IEEE_Arith
  45 #ifndef NO_IEEE_Scale
  46 #define Avoid_Underflow
  47 #undef tinytens
  48 /* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */
  49 /* flag unnecessarily.  It leads to a song and dance at the end of strtod. */
  50 static CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128,
  51     9007199254740992.e-256
  52     };
  53 #endif
  54 #endif
  55
  56 #ifdef Honor_FLT_ROUNDS
  57 #define Rounding rounding
  58 #undef Check_FLT_ROUNDS
  59 #define Check_FLT_ROUNDS
  60 #else
  61 #define Rounding Flt_Rounds
  62 #endif
  63
  64 //#ifndef __HAVE_LONG_DOUBLE
  65 //__strong_alias(_strtold, strtod)
  66 //__weak_alias(strtold, _strtold)
  67 //#endif
  68
  69 #if defined(_MSC_VER)           /* Handle Microsoft VC++ compiler specifics. */
  70 // Disable: warning C4700: uninitialized local variable 'xx' used
  71 #pragma warning ( disable : 4700 )
  72 #endif  /* defined(_MSC_VER) */
  73
  74 double
  75 strtod(CONST char *s00, char **se)
  76 {
  77 #ifdef Avoid_Underflow
  78   int scale;
  79   #endif
  80   int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign,
  81       e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign;
  82   CONST char *s, *s0, *s1;
  83   double aadj, aadj1, adj, rv, rv0;
  84   Long L;
  85   ULong y, z;
  86   Bigint *bb = NULL, *bb1, *bd0;
  87   Bigint *bd = NULL, *bs = NULL, *delta = NULL; /* pacify gcc */
  88 #ifdef SET_INEXACT
  89   int inexact, oldinexact;
  90 #endif
  91 #ifdef Honor_FLT_ROUNDS
  92   int rounding;
  93 #endif
  94
  95   sign = nz0 = nz = decpt = 0;
  96   dval(rv) = 0.;
  97   for(s = s00;;s++) {
  98     switch(*s) {
  99       case '-':
 100         sign = 1;
 101         /* FALLTHROUGH */
 102       case '+':
 103         if (*++s)
 104           goto break2;
 105         /* FALLTHROUGH */
 106       case 0:
 107         goto ret0;
 108       case '\t':
 109       case '\n':
 110       case '\v':
 111       case '\f':
 112       case '\r':
 113       case ' ':
 114         continue;
 115       default:
 116         goto break2;
 117     }
 118   }
 119  break2:
 120   if (*s == '0') {
 121 #ifndef NO_HEX_FP
 122     {
 123     static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI };
 124     Long expt;
 125     ULong bits[2];
 126     switch(s[1]) {
 127       case 'x':
 128       case 'X':
 129       {
 130 #if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD)
 131       FPI fpi1 = fpi;
 132       switch(fegetround()) {
 133         case FE_TOWARDZERO: fpi1.rounding = 0; break;
 134         case FE_UPWARD: fpi1.rounding = 2; break;
 135         case FE_DOWNWARD: fpi1.rounding = 3;
 136         }
 137 #else
 138 #endif
 139       switch((i = gethex(&s, &fpi, &expt, &bb, sign)) & STRTOG_Retmask) {
 140         case STRTOG_NoNumber:
 141         s = s00;
 142         sign = 0;
 143         /* FALLTHROUGH */
 144         case STRTOG_Zero:
 145         break;
 146         default:
 147         if (bb) {
 148           copybits(bits, fpi.nbits, bb);
 149           Bfree(bb);
 150           }
 151         ULtod((/* LINTED */(U*)&rv)->L, bits, expt, i);
 152         }}
 153       goto ret;
 154       }
 155     }
 156 #endif
 157     nz0 = 1;
 158     while(*++s == '0') ;
 159     if (!*s)
 160       goto ret;
 161   }
 162   s0 = s;
 163   y = z = 0;
 164   for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
 165     if (nd < 9)
 166       y = 10*y + c - '0';
 167     else if (nd < 16)
 168       z = 10*z + c - '0';
 169   nd0 = nd;
 170 #ifdef USE_LOCALE
 171   if (c == *localeconv()->decimal_point)
 172 #else
 173   if (c == '.')
 174 #endif
 175     {
 176     decpt = 1;
 177     c = *++s;
 178     if (!nd) {
 179       for(; c == '0'; c = *++s)
 180         nz++;
 181       if (c > '0' && c <= '9') {
 182         s0 = s;
 183         nf += nz;
 184         nz = 0;
 185         goto have_dig;
 186         }
 187       goto dig_done;
 188       }
 189     for(; c >= '0' && c <= '9'; c = *++s) {
 190  have_dig:
 191       nz++;
 192       if (c -= '0') {
 193         nf += nz;
 194         for(i = 1; i < nz; i++)
 195           if (nd++ < 9)
 196             y *= 10;
 197           else if (nd <= DBL_DIG + 1)
 198             z *= 10;
 199         if (nd++ < 9)
 200           y = 10*y + c;
 201         else if (nd <= DBL_DIG + 1)
 202           z = 10*z + c;
 203         nz = 0;
 204         }
 205       }
 206     }
 207  dig_done:
 208   e = 0;
 209   if (c == 'e' || c == 'E') {
 210     if (!nd && !nz && !nz0) {
 211       goto ret0;
 212       }
 213     s00 = s;
 214     esign = 0;
 215     switch(c = *++s) {
 216       case '-':
 217         esign = 1;
 218         /* FALLTHROUGH */
 219       case '+':
 220         c = *++s;
 221       }
 222     if (c >= '0' && c <= '9') {
 223       while(c == '0')
 224         c = *++s;
 225       if (c > '0' && c <= '9') {
 226         L = c - '0';
 227         s1 = s;
 228         while((c = *++s) >= '0' && c <= '9')
 229           L = 10*L + c - '0';
 230         if (s - s1 > 8 || L > 19999)
 231           /* Avoid confusion from exponents
 232            * so large that e might overflow.
 233            */
 234           e = 19999; /* safe for 16 bit ints */
 235         else
 236           e = (int)L;
 237         if (esign)
 238           e = -e;
 239         }
 240       else
 241         e = 0;
 242       }
 243     else
 244       s = s00;
 245     }
 246   if (!nd) {
 247     if (!nz && !nz0) {
 248 #ifdef INFNAN_CHECK
 249       /* Check for Nan and Infinity */
 250 #ifndef No_Hex_NaN
 251       ULong bits[2];
 252       static FPI fpinan = /* only 52 explicit bits */
 253         { 52, 1-1023-53+1, 2046-1023-53+1, 1, SI };
 254 #endif  // No_Hex_NaN
 255       if (!decpt)
 256        switch(c) {
 257         case 'i':
 258         case 'I':
 259         if (match(&s,"nf")) {
 260           --s;
 261           if (!match(&s,"inity"))
 262             ++s;
 263           word0(rv) = 0x7ff00000;
 264           word1(rv) = 0;
 265           goto ret;
 266           }
 267         break;
 268         case 'n':
 269         case 'N':
 270         if (match(&s, "an")) {
 271 #ifndef No_Hex_NaN
 272           if (*s == '(' /*)*/
 273            && hexnan(&s, &fpinan, bits)
 274               == STRTOG_NaNbits) {
 275             word0(rv) = (UINT32)(0x7ff00000U | bits[1]);
 276             word1(rv) = (UINT32)bits[0];
 277             }
 278           else {
 279 #endif
 280             word0(rv) = NAN_WORD0;
 281             word1(rv) = NAN_WORD1;
 282 #ifndef No_Hex_NaN
 283             }
 284 #endif
 285           goto ret;
 286           }
 287         }
 288 #endif /* INFNAN_CHECK */
 289  ret0:
 290       s = s00;
 291       sign = 0;
 292       }
 293     goto ret;
 294     }
 295   e1 = e -= nf;
 296
 297   /* Now we have nd0 digits, starting at s0, followed by a
 298    * decimal point, followed by nd-nd0 digits.  The number we're
 299    * after is the integer represented by those digits times
 300    * 10**e */
 301
 302   if (!nd0)
 303     nd0 = nd;
 304   k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
 305   dval(rv) = (double)y;
 306   if (k > 9) {
 307 #ifdef SET_INEXACT
 308     if (k > DBL_DIG)
 309       oldinexact = get_inexact();
 310 #endif
 311     dval(rv) = tens[k - 9] * dval(rv) + z;
 312     }
 313   bd0 = 0;
 314   if (nd <= DBL_DIG
 315 #ifndef RND_PRODQUOT
 316 #ifndef Honor_FLT_ROUNDS
 317     && Flt_Rounds == 1
 318 #endif
 319 #endif
 320       ) {
 321     if (!e)
 322       goto ret;
 323     if (e > 0) {
 324       if (e <= Ten_pmax) {
 325 #ifdef VAX
 326         goto vax_ovfl_check;
 327 #else
 328 #ifdef Honor_FLT_ROUNDS
 329         /* round correctly FLT_ROUNDS = 2 or 3 */
 330         if (sign) {
 331           rv = -rv;
 332           sign = 0;
 333           }
 334 #endif
 335         /* rv = */ rounded_product(dval(rv), tens[e]);
 336         goto ret;
 337 #endif
 338         }
 339       i = DBL_DIG - nd;
 340       if (e <= Ten_pmax + i) {
 341         /* A fancier test would sometimes let us do
 342          * this for larger i values.
 343          */
 344 #ifdef Honor_FLT_ROUNDS
 345         /* round correctly FLT_ROUNDS = 2 or 3 */
 346         if (sign) {
 347           rv = -rv;
 348           sign = 0;
 349           }
 350 #endif
 351         e -= i;
 352         dval(rv) *= tens[i];
 353 #ifdef VAX
 354         /* VAX exponent range is so narrow we must
 355          * worry about overflow here...
 356          */
 357  vax_ovfl_check:
 358         word0(rv) -= P*Exp_msk1;
 359         /* rv = */ rounded_product(dval(rv), tens[e]);
 360         if ((word0(rv) & Exp_mask)
 361          > Exp_msk1*(DBL_MAX_EXP+Bias-1-P))
 362           goto ovfl;
 363         word0(rv) += P*Exp_msk1;
 364 #else
 365         /* rv = */ rounded_product(dval(rv), tens[e]);
 366 #endif
 367         goto ret;
 368         }
 369       }
 370 #ifndef Inaccurate_Divide
 371     else if (e >= -Ten_pmax) {
 372 #ifdef Honor_FLT_ROUNDS
 373       /* round correctly FLT_ROUNDS = 2 or 3 */
 374       if (sign) {
 375         rv = -rv;
 376         sign = 0;
 377         }
 378 #endif
 379       /* rv = */ rounded_quotient(dval(rv), tens[-e]);
 380       goto ret;
 381       }
 382 #endif
 383     }
 384   e1 += nd - k;
 385
 386 #ifdef IEEE_Arith
 387 #ifdef SET_INEXACT
 388   inexact = 1;
 389   if (k <= DBL_DIG)
 390     oldinexact = get_inexact();
 391 #endif
 392 #ifdef Avoid_Underflow
 393   scale = 0;
 394 #endif
 395 #ifdef Honor_FLT_ROUNDS
 396   if ((rounding = Flt_Rounds) >= 2) {
 397     if (sign)
 398       rounding = rounding == 2 ? 0 : 2;
 399     else
 400       if (rounding != 2)
 401         rounding = 0;
 402     }
 403 #endif
 404 #endif /*IEEE_Arith*/
 405
 406   /* Get starting approximation = rv * 10**e1 */
 407
 408   if (e1 > 0) {
 409     if ( (i = e1 & 15) !=0)
 410       dval(rv) *= tens[i];
 411     if (e1 &= ~15) {
 412       if (e1 > DBL_MAX_10_EXP) {
 413  ovfl:
 414 #ifndef NO_ERRNO
 415         errno = ERANGE;
 416 #endif
 417         /* Can't trust HUGE_VAL */
 418 #ifdef IEEE_Arith
 419 #ifdef Honor_FLT_ROUNDS
 420         switch(rounding) {
 421           case 0: /* toward 0 */
 422           case 3: /* toward -infinity */
 423           word0(rv) = Big0;
 424           word1(rv) = Big1;
 425           break;
 426           default:
 427           word0(rv) = Exp_mask;
 428           word1(rv) = 0;
 429           }
 430 #else /*Honor_FLT_ROUNDS*/
 431         word0(rv) = Exp_mask;
 432         word1(rv) = 0;
 433 #endif /*Honor_FLT_ROUNDS*/
 434 #ifdef SET_INEXACT
 435         /* set overflow bit */
 436         dval(rv0) = 1e300;
 437         dval(rv0) *= dval(rv0);
 438 #endif
 439 #else /*IEEE_Arith*/
 440         word0(rv) = Big0;
 441         word1(rv) = Big1;
 442 #endif /*IEEE_Arith*/
 443         if (bd0)
 444           goto retfree;
 445         goto ret;
 446         }
 447       e1 = (unsigned int)e1 >> 4;
 448       for(j = 0; e1 > 1; j++, e1 = (unsigned int)e1 >> 1)
 449         if (e1 & 1)
 450           dval(rv) *= bigtens[j];
 451     /* The last multiplication could overflow. */
 452       word0(rv) -= P*Exp_msk1;
 453       dval(rv) *= bigtens[j];
 454       if ((z = word0(rv) & Exp_mask)
 455        > Exp_msk1*(DBL_MAX_EXP+Bias-P))
 456         goto ovfl;
 457       if (z > Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) {
 458         /* set to largest number */
 459         /* (Can't trust DBL_MAX) */
 460         word0(rv) = Big0;
 461         word1(rv) = Big1;
 462         }
 463       else
 464         word0(rv) += P*Exp_msk1;
 465       }
 466     }
 467   else if (e1 < 0) {
 468     e1 = -e1;
 469     if ( (i = e1 & 15) !=0)
 470       dval(rv) /= tens[i];
 471     if (e1 >>= 4) {
 472       if (e1 >= 1 << n_bigtens)
 473         goto undfl;
 474 #ifdef Avoid_Underflow
 475       if (e1 & Scale_Bit)
 476         scale = 2*P;
 477       for(j = 0; e1 > 0; j++, e1 = (unsigned int)e1 >> 1)
 478         if (e1 & 1)
 479           dval(rv) *= tinytens[j];
 480       if (scale && (j = 2*P + 1 - ((word0(rv) & Exp_mask)
 481             >> Exp_shift)) > 0) {
 482         /* scaled rv is denormal; zap j low bits */
 483         if (j >= 32) {
 484           word1(rv) = 0;
 485           if (j >= 53)
 486            word0(rv) = (P+2)*Exp_msk1;
 487           else
 488            word0(rv) &= 0xffffffff << (j-32);
 489           }
 490         else
 491           word1(rv) &= 0xffffffff << j;
 492         }
 493 #else
 494       for(j = 0; e1 > 1; j++, e1 >>= 1)
 495         if (e1 & 1)
 496           dval(rv) *= tinytens[j];
 497       /* The last multiplication could underflow. */
 498       dval(rv0) = dval(rv);
 499       dval(rv) *= tinytens[j];
 500       if (!dval(rv)) {
 501         dval(rv) = 2.*dval(rv0);
 502         dval(rv) *= tinytens[j];
 503 #endif
 504         if (!dval(rv)) {
 505  undfl:
 506           dval(rv) = 0.;
 507 #ifndef NO_ERRNO
 508           errno = ERANGE;
 509 #endif
 510           if (bd0)
 511             goto retfree;
 512           goto ret;
 513           }
 514 #ifndef Avoid_Underflow
 515         word0(rv) = Tiny0;
 516         word1(rv) = Tiny1;
 517         /* The refinement below will clean
 518          * this approximation up.
 519          */
 520         }
 521 #endif
 522       }
 523     }
 524
 525   /* Now the hard part -- adjusting rv to the correct value.*/
 526
 527   /* Put digits into bd: true value = bd * 10^e */
 528
 529   bd0 = s2b(s0, nd0, nd, y);
 530   if (bd0 == NULL)
 531     goto ovfl;
 532
 533   for(;;) {
 534     bd = Balloc(bd0->k);
 535     if (bd == NULL)
 536       goto ovfl;
 537     Bcopy(bd, bd0);
 538     bb = d2b(dval(rv), &bbe, &bbbits);  /* rv = bb * 2^bbe */
 539     if (bb == NULL)
 540       goto ovfl;
 541     bs = i2b(1);
 542     if (bs == NULL)
 543       goto ovfl;
 544
 545     if (e >= 0) {
 546       bb2 = bb5 = 0;
 547       bd2 = bd5 = e;
 548       }
 549     else {
 550       bb2 = bb5 = -e;
 551       bd2 = bd5 = 0;
 552       }
 553     if (bbe >= 0)
 554       bb2 += bbe;
 555     else
 556       bd2 -= bbe;
 557     bs2 = bb2;
 558 #ifdef Honor_FLT_ROUNDS
 559     if (rounding != 1)
 560       bs2++;
 561 #endif
 562 #ifdef Avoid_Underflow
 563     j = bbe - scale;
 564     i = j + bbbits - 1; /* logb(rv) */
 565     if (i < Emin) /* denormal */
 566       j += P - Emin;
 567     else
 568       j = P + 1 - bbbits;
 569 #else /*Avoid_Underflow*/
 570 #ifdef Sudden_Underflow
 571 #ifdef IBM
 572     j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
 573 #else
 574     j = P + 1 - bbbits;
 575 #endif
 576 #else /*Sudden_Underflow*/
 577     j = bbe;
 578     i = j + bbbits - 1; /* logb(rv) */
 579     if (i < Emin) /* denormal */
 580       j += P - Emin;
 581     else
 582       j = P + 1 - bbbits;
 583 #endif /*Sudden_Underflow*/
 584 #endif /*Avoid_Underflow*/
 585     bb2 += j;
 586     bd2 += j;
 587 #ifdef Avoid_Underflow
 588     bd2 += scale;
 589 #endif
 590     i = bb2 < bd2 ? bb2 : bd2;
 591     if (i > bs2)
 592       i = bs2;
 593     if (i > 0) {
 594       bb2 -= i;
 595       bd2 -= i;
 596       bs2 -= i;
 597       }
 598     if (bb5 > 0) {
 599       bs = pow5mult(bs, bb5);
 600       if (bs == NULL)
 601         goto ovfl;
 602       bb1 = mult(bs, bb);
 603       if (bb1 == NULL)
 604         goto ovfl;
 605       Bfree(bb);
 606       bb = bb1;
 607       }
 608     if (bb2 > 0) {
 609       bb = lshift(bb, bb2);
 610       if (bb == NULL)
 611         goto ovfl;
 612       }
 613     if (bd5 > 0) {
 614       bd = pow5mult(bd, bd5);
 615       if (bd == NULL)
 616         goto ovfl;
 617       }
 618     if (bd2 > 0) {
 619       bd = lshift(bd, bd2);
 620       if (bd == NULL)
 621         goto ovfl;
 622       }
 623     if (bs2 > 0) {
 624       bs = lshift(bs, bs2);
 625       if (bs == NULL)
 626         goto ovfl;
 627       }
 628     delta = diff(bb, bd);
 629     if (delta == NULL)
 630       goto ovfl;
 631     dsign = delta->sign;
 632     delta->sign = 0;
 633     i = cmp(delta, bs);
 634 #ifdef Honor_FLT_ROUNDS
 635     if (rounding != 1) {
 636       if (i < 0) {
 637         /* Error is less than an ulp */
 638         if (!delta->x[0] && delta->wds <= 1) {
 639           /* exact */
 640 #ifdef SET_INEXACT
 641           inexact = 0;
 642 #endif
 643           break;
 644           }
 645         if (rounding) {
 646           if (dsign) {
 647             adj = 1.;
 648             goto apply_adj;
 649             }
 650           }
 651         else if (!dsign) {
 652           adj = -1.;
 653           if (!word1(rv)
 654            && !(word0(rv) & Frac_mask)) {
 655             y = word0(rv) & Exp_mask;
 656 #ifdef Avoid_Underflow
 657             if (!scale || y > 2*P*Exp_msk1)
 658 #else
 659             if (y)
 660 #endif
 661               {
 662               delta = lshift(delta,Log2P);
 663               if (cmp(delta, bs) <= 0)
 664               adj = -0.5;
 665               }
 666             }
 667  apply_adj:
 668 #ifdef Avoid_Underflow
 669           if (scale && (y = word0(rv) & Exp_mask)
 670             <= 2*P*Exp_msk1)
 671             word0(adj) += (2*P+1)*Exp_msk1 - y;
 672 #else
 673 #ifdef Sudden_Underflow
 674           if ((word0(rv) & Exp_mask) <=
 675               P*Exp_msk1) {
 676             word0(rv) += P*Exp_msk1;
 677             dval(rv) += adj*ulp(dval(rv));
 678             word0(rv) -= P*Exp_msk1;
 679             }
 680           else
 681 #endif /*Sudden_Underflow*/
 682 #endif /*Avoid_Underflow*/
 683           dval(rv) += adj*ulp(dval(rv));
 684           }
 685         break;
 686         }
 687       adj = ratio(delta, bs);
 688       if (adj < 1.)
 689         adj = 1.;
 690       if (adj <= 0x7ffffffe) {
 691         /* adj = rounding ? ceil(adj) : floor(adj); */
 692         y = adj;
 693         if (y != adj) {
 694           if (!((rounding>>1) ^ dsign))
 695             y++;
 696           adj = y;
 697           }
 698         }
 699 #ifdef Avoid_Underflow
 700       if (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
 701         word0(adj) += (2*P+1)*Exp_msk1 - y;
 702 #else
 703 #ifdef Sudden_Underflow
 704       if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
 705         word0(rv) += P*Exp_msk1;
 706         adj *= ulp(dval(rv));
 707         if (dsign)
 708           dval(rv) += adj;
 709         else
 710           dval(rv) -= adj;
 711         word0(rv) -= P*Exp_msk1;
 712         goto cont;
 713         }
 714 #endif /*Sudden_Underflow*/
 715 #endif /*Avoid_Underflow*/
 716       adj *= ulp(dval(rv));
 717       if (dsign)
 718         dval(rv) += adj;
 719       else
 720         dval(rv) -= adj;
 721       goto cont;
 722       }
 723 #endif /*Honor_FLT_ROUNDS*/
 724
 725     if (i < 0) {
 726       /* Error is less than half an ulp -- check for
 727        * special case of mantissa a power of two.
 728        */
 729       if (dsign || word1(rv) || word0(rv) & Bndry_mask
 730 #ifdef IEEE_Arith
 731 #ifdef Avoid_Underflow
 732        || (word0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1
 733 #else
 734        || (word0(rv) & Exp_mask) <= Exp_msk1
 735 #endif
 736 #endif
 737         ) {
 738 #ifdef SET_INEXACT
 739         if (!delta->x[0] && delta->wds <= 1)
 740           inexact = 0;
 741 #endif
 742         break;
 743         }
 744       if (!delta->x[0] && delta->wds <= 1) {
 745         /* exact result */
 746 #ifdef SET_INEXACT
 747         inexact = 0;
 748 #endif
 749         break;
 750         }
 751       delta = lshift(delta,Log2P);
 752       if (cmp(delta, bs) > 0)
 753         goto drop_down;
 754       break;
 755       }
 756     if (i == 0) {
 757       /* exactly half-way between */
 758       if (dsign) {
 759         if ((word0(rv) & Bndry_mask1) == Bndry_mask1
 760          &&  word1(rv) == (
 761 #ifdef Avoid_Underflow
 762       (scale && (y = word0(rv) & Exp_mask) <= 2*P*Exp_msk1)
 763     ? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) :
 764 #endif
 765                0xffffffff)) {
 766           /*boundary case -- increment exponent*/
 767           word0(rv) = (word0(rv) & Exp_mask)
 768             + Exp_msk1
 769 #ifdef IBM
 770             | Exp_msk1 >> 4
 771 #endif
 772             ;
 773           word1(rv) = 0;
 774 #ifdef Avoid_Underflow
 775           dsign = 0;
 776 #endif
 777           break;
 778           }
 779         }
 780       else if (!(word0(rv) & Bndry_mask) && !word1(rv)) {
 781  drop_down:
 782         /* boundary case -- decrement exponent */
 783 #ifdef Sudden_Underflow /*{{*/
 784         L = word0(rv) & Exp_mask;
 785 #ifdef IBM
 786         if (L <  Exp_msk1)
 787 #else
 788 #ifdef Avoid_Underflow
 789         if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1))
 790 #else
 791         if (L <= Exp_msk1)
 792 #endif /*Avoid_Underflow*/
 793 #endif /*IBM*/
 794           goto undfl;
 795         L -= Exp_msk1;
 796 #else /*Sudden_Underflow}{*/
 797 #ifdef Avoid_Underflow
 798         if (scale) {
 799           L = word0(rv) & Exp_mask;
 800           if (L <= (2*P+1)*Exp_msk1) {
 801             if (L > (P+2)*Exp_msk1)
 802               /* round even ==> */
 803               /* accept rv */
 804               break;
 805             /* rv = smallest denormal */
 806             goto undfl;
 807             }
 808           }
 809 #endif /*Avoid_Underflow*/
 810         L = (word0(rv) & Exp_mask) - Exp_msk1;
 811 #endif /*Sudden_Underflow}*/
 812         word0(rv) = (UINT32)(L | Bndry_mask1);
 813         word1(rv) = 0xffffffffU;
 814 #ifdef IBM
 815         goto cont;
 816 #else
 817         break;
 818 #endif
 819         }
 820 #ifndef ROUND_BIASED
 821       if (!(word1(rv) & LSB))
 822         break;
 823 #endif
 824       if (dsign)
 825         dval(rv) += ulp(dval(rv));
 826 #ifndef ROUND_BIASED
 827       else {
 828         dval(rv) -= ulp(dval(rv));
 829 #ifndef Sudden_Underflow
 830         if (!dval(rv))
 831           goto undfl;
 832 #endif
 833         }
 834 #ifdef Avoid_Underflow
 835       dsign = 1 - dsign;
 836 #endif
 837 #endif
 838       break;
 839       }
 840     if ((aadj = ratio(delta, bs)) <= 2.) {
 841       if (dsign)
 842         aadj = aadj1 = 1.;
 843       else if (word1(rv) || word0(rv) & Bndry_mask) {
 844 #ifndef Sudden_Underflow
 845         if (word1(rv) == Tiny1 && !word0(rv))
 846           goto undfl;
 847 #endif
 848         aadj = 1.;
 849         aadj1 = -1.;
 850         }
 851       else {
 852         /* special case -- power of FLT_RADIX to be */
 853         /* rounded down... */
 854
 855         if (aadj < 2./FLT_RADIX)
 856           aadj = 1./FLT_RADIX;
 857         else
 858           aadj *= 0.5;
 859         aadj1 = -aadj;
 860         }
 861       }
 862     else {
 863       aadj *= 0.5;
 864       aadj1 = dsign ? aadj : -aadj;
 865 #ifdef Check_FLT_ROUNDS
 866       switch(Rounding) {
 867         case 2: /* towards +infinity */
 868           aadj1 -= 0.5;
 869           break;
 870         case 0: /* towards 0 */
 871         case 3: /* towards -infinity */
 872           aadj1 += 0.5;
 873         }
 874 #else
 875       if (Flt_Rounds == 0)
 876         aadj1 += 0.5;
 877 #endif /*Check_FLT_ROUNDS*/
 878       }
 879     y = word0(rv) & Exp_mask;
 880
 881     /* Check for overflow */
 882
 883     if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) {
 884       dval(rv0) = dval(rv);
 885       word0(rv) -= P*Exp_msk1;
 886       adj = aadj1 * ulp(dval(rv));
 887       dval(rv) += adj;
 888       if ((word0(rv) & Exp_mask) >=
 889           Exp_msk1*(DBL_MAX_EXP+Bias-P)) {
 890         if (word0(rv0) == Big0 && word1(rv0) == Big1)
 891           goto ovfl;
 892         word0(rv) = Big0;
 893         word1(rv) = Big1;
 894         goto cont;
 895         }
 896       else
 897         word0(rv) += P*Exp_msk1;
 898       }
 899     else {
 900 #ifdef Avoid_Underflow
 901       if (scale && y <= 2*P*Exp_msk1) {
 902         if (aadj <= 0x7fffffff) {
 903           if ((z = (uint32_t)aadj) == 0)
 904             z = 1;
 905           aadj = (double)z;
 906           aadj1 = dsign ? aadj : -aadj;
 907           }
 908         word0(aadj1) += (UINT32)((2*P+1)*Exp_msk1 - y);
 909         }
 910       adj = aadj1 * ulp(dval(rv));
 911       dval(rv) += adj;
 912 #else
 913 #ifdef Sudden_Underflow
 914       if ((word0(rv) & Exp_mask) <= P*Exp_msk1) {
 915         dval(rv0) = dval(rv);
 916         word0(rv) += P*Exp_msk1;
 917         adj = aadj1 * ulp(dval(rv));
 918         dval(rv) += adj;
 919 #ifdef IBM
 920         if ((word0(rv) & Exp_mask) <  P*Exp_msk1)
 921 #else
 922         if ((word0(rv) & Exp_mask) <= P*Exp_msk1)
 923 #endif
 924           {
 925           if (word0(rv0) == Tiny0
 926            && word1(rv0) == Tiny1)
 927             goto undfl;
 928           word0(rv) = Tiny0;
 929           word1(rv) = Tiny1;
 930           goto cont;
 931           }
 932         else
 933           word0(rv) -= P*Exp_msk1;
 934         }
 935       else {
 936         adj = aadj1 * ulp(dval(rv));
 937         dval(rv) += adj;
 938         }
 939 #else /*Sudden_Underflow*/
 940       /* Compute adj so that the IEEE rounding rules will
 941        * correctly round rv + adj in some half-way cases.
 942        * If rv * ulp(rv) is denormalized (i.e.,
 943        * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid
 944        * trouble from bits lost to denormalization;
 945        * example: 1.2e-307 .
 946        */
 947       if (y <= (P-1)*Exp_msk1 && aadj > 1.) {
 948         aadj1 = (double)(int)(aadj + 0.5);
 949         if (!dsign)
 950           aadj1 = -aadj1;
 951         }
 952       adj = aadj1 * ulp(dval(rv));
 953       dval(rv) += adj;
 954 #endif /*Sudden_Underflow*/
 955 #endif /*Avoid_Underflow*/
 956       }
 957     z = word0(rv) & Exp_mask;
 958 #ifndef SET_INEXACT
 959 #ifdef Avoid_Underflow
 960     if (!scale)
 961 #endif
 962     if (y == z) {
 963       /* Can we stop now? */
 964       L = (Long)aadj;
 965       aadj -= L;
 966       /* The tolerances below are conservative. */
 967       if (dsign || word1(rv) || word0(rv) & Bndry_mask) {
 968         if (aadj < .4999999 || aadj > .5000001)
 969           break;
 970         }
 971       else if (aadj < .4999999/FLT_RADIX)
 972         break;
 973       }
 974 #endif
 975  cont:
 976     Bfree(bb);
 977     Bfree(bd);
 978     Bfree(bs);
 979     Bfree(delta);
 980     }
 981 #ifdef SET_INEXACT
 982   if (inexact) {
 983     if (!oldinexact) {
 984       word0(rv0) = Exp_1 + (70 << Exp_shift);
 985       word1(rv0) = 0;
 986       dval(rv0) += 1.;
 987       }
 988     }
 989   else if (!oldinexact)
 990     clear_inexact();
 991 #endif
 992 #ifdef Avoid_Underflow
 993   if (scale) {
 994     word0(rv0) = Exp_1 - 2*P*Exp_msk1;
 995     word1(rv0) = 0;
 996     dval(rv) *= dval(rv0);
 997 #ifndef NO_ERRNO
 998     /* try to avoid the bug of testing an 8087 register value */
 999     if (word0(rv) == 0 && word1(rv) == 0)
1000       errno = ERANGE;
1001 #endif
1002     }
1003 #endif /* Avoid_Underflow */
1004 #ifdef SET_INEXACT
1005   if (inexact && !(word0(rv) & Exp_mask)) {
1006     /* set underflow bit */
1007     dval(rv0) = 1e-300;
1008     dval(rv0) *= dval(rv0);
1009     }
1010 #endif
1011  retfree:
1012   Bfree(bb);
1013   Bfree(bd);
1014   Bfree(bs);
1015   Bfree(bd0);
1016   Bfree(delta);
1017  ret:
1018   if (se)
1019     *se = __UNCONST(s);
1020   return sign ? -dval(rv) : dval(rv);
1021 }
1022