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