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