lib/libc/gdtoa/strtod.c

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