StdLib/LibC/gdtoa/strtod.c

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