newlib/libc/stdlib/strtod.c

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