contrib/gdtoa/strtod.c

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