lib/libc/gdtoa/strtod.c

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