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