contrib/gdtoa/strtod.c

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