contrib/gdtoa/strtod.c

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