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