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