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