1 /* 2 FUNCTION 3 <<strtod>>, <<strtodf>>---string to double or float 4 5 INDEX 6 strtod 7 INDEX 8 _strtod_r 9 INDEX 10 strtodf 11 12 ANSI_SYNOPSIS 13 #include <stdlib.h> 14 double strtod(const char *<[str]>, char **<[tail]>); 15 float strtodf(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 strtodf(<[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 <<strtodf>> 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) 1991 by AT&T. 77 * 78 * Permission to use, copy, modify, and distribute this software for any 79 * purpose without fee is hereby granted, provided that this entire notice 80 * is included in all copies of any software which is or includes a copy 81 * or modification of this software and in all copies of the supporting 82 * documentation for such software. 83 * 84 * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED 85 * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY 86 * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY 87 * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. 88 * 89 ***************************************************************/ 90 91 /* Please send bug reports to 92 David M. Gay 93 AT&T Bell Laboratories, Room 2C-463 94 600 Mountain Avenue 95 Murray Hill, NJ 07974-2070 96 U.S.A. 97 dmg@research.att.com or research!dmg 98 */ 99 100 #include <string.h> 101 #include <float.h> 102 #include <errno.h> 103 #include "mprec.h" 104 105 double 106 _DEFUN (_strtod_r, (ptr, s00, se), 107 struct _Jv_reent *ptr _AND 108 _CONST char *s00 _AND 109 char **se) 110 { 111 int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, e1, esign, i, j, 112 k, nd, nd0, nf, nz, nz0, sign; 113 int digits = 0; /* Number of digits found in fraction part. */ 114 long e; 115 _CONST char *s, *s0, *s1; 116 double aadj, aadj1, adj; 117 long L; 118 unsigned long y, z; 119 union double_union rv, rv0; 120 121 _Jv_Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; 122 sign = nz0 = nz = 0; 123 rv.d = 0.; 124 for (s = s00;; s++) 125 switch (*s) 126 { 127 case '-': 128 sign = 1; 129 /* no break */ 130 case '+': 131 if (*++s) 132 goto break2; 133 /* no break */ 134 case 0: 135 s = s00; 136 goto ret; 137 case '\t': 138 case '\n': 139 case '\v': 140 case '\f': 141 case '\r': 142 case ' ': 143 continue; 144 default: 145 goto break2; 146 } 147 break2: 148 if (*s == '0') 149 { 150 digits++; 151 nz0 = 1; 152 while (*++s == '0') 153 digits++; 154 if (!*s) 155 goto ret; 156 } 157 s0 = s; 158 y = z = 0; 159 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) 160 { 161 digits++; 162 if (nd < 9) 163 y = 10 * y + c - '0'; 164 else if (nd < 16) 165 z = 10 * z + c - '0'; 166 } 167 nd0 = nd; 168 if (c == '.') 169 { 170 c = *++s; 171 if (!nd) 172 { 173 for (; c == '0'; c = *++s) 174 { 175 digits++; 176 nz++; 177 } 178 if (c > '0' && c <= '9') 179 { 180 digits++; 181 s0 = s; 182 nf += nz; 183 nz = 0; 184 goto have_dig; 185 } 186 goto dig_done; 187 } 188 for (; c >= '0' && c <= '9'; c = *++s) 189 { 190 digits++; 191 have_dig: 192 nz++; 193 if (c -= '0') 194 { 195 nf += nz; 196 for (i = 1; i < nz; i++) 197 if (nd++ < 9) 198 y *= 10; 199 else if (nd <= DBL_DIG + 1) 200 z *= 10; 201 if (nd++ < 9) 202 y = 10 * y + c; 203 else if (nd <= DBL_DIG + 1) 204 z = 10 * z + c; 205 nz = 0; 206 } 207 } 208 } 209 dig_done: 210 e = 0; 211 if (c == 'e' || c == 'E') 212 { 213 if (!nd && !nz && !nz0) 214 { 215 s = s00; 216 goto ret; 217 } 218 s00 = s; 219 esign = 0; 220 switch (c = *++s) 221 { 222 case '-': 223 esign = 1; 224 case '+': 225 c = *++s; 226 } 227 if (c >= '0' && c <= '9') 228 { 229 while (c == '0') 230 c = *++s; 231 if (c > '0' && c <= '9') 232 { 233 e = c - '0'; 234 s1 = s; 235 while ((c = *++s) >= '0' && c <= '9') 236 e = 10 * e + c - '0'; 237 if (s - s1 > 8) 238 /* Avoid confusion from exponents 239 * so large that e might overflow. 240 */ 241 e = 9999999L; 242 if (esign) 243 e = -e; 244 } 245 } 246 else 247 { 248 /* No exponent after an 'E' : that's an error. */ 249 ptr->_errno = EINVAL; 250 e = 0; 251 s = s00; 252 goto ret; 253 } 254 } 255 if (!nd) 256 { 257 if (!nz && !nz0) 258 s = s00; 259 goto ret; 260 } 261 e1 = e -= nf; 262 263 /* Now we have nd0 digits, starting at s0, followed by a 264 * decimal point, followed by nd-nd0 digits. The number we're 265 * after is the integer represented by those digits times 266 * 10**e */ 267 268 if (!nd0) 269 nd0 = nd; 270 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; 271 rv.d = y; 272 if (k > 9) 273 rv.d = tens[k - 9] * rv.d + z; 274 bd0 = 0; 275 if (nd <= DBL_DIG 276 #ifndef RND_PRODQUOT 277 && FLT_ROUNDS == 1 278 #endif 279 ) 280 { 281 if (!e) 282 goto ret; 283 if (e > 0) 284 { 285 if (e <= Ten_pmax) 286 { 287 #ifdef VAX 288 goto vax_ovfl_check; 289 #else 290 /* rv.d = */ rounded_product (rv.d, tens[e]); 291 goto ret; 292 #endif 293 } 294 i = DBL_DIG - nd; 295 if (e <= Ten_pmax + i) 296 { 297 /* A fancier test would sometimes let us do 298 * this for larger i values. 299 */ 300 e -= i; 301 rv.d *= tens[i]; 302 #ifdef VAX 303 /* VAX exponent range is so narrow we must 304 * worry about overflow here... 305 */ 306 vax_ovfl_check: 307 word0 (rv) -= P * Exp_msk1; 308 /* rv.d = */ rounded_product (rv.d, tens[e]); 309 if ((word0 (rv) & Exp_mask) 310 > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) 311 goto ovfl; 312 word0 (rv) += P * Exp_msk1; 313 #else 314 /* rv.d = */ rounded_product (rv.d, tens[e]); 315 #endif 316 goto ret; 317 } 318 } 319 #ifndef Inaccurate_Divide 320 else if (e >= -Ten_pmax) 321 { 322 /* rv.d = */ rounded_quotient (rv.d, tens[-e]); 323 goto ret; 324 } 325 #endif 326 } 327 e1 += nd - k; 328 329 /* Get starting approximation = rv.d * 10**e1 */ 330 331 if (e1 > 0) 332 { 333 if ((i = e1 & 15)) 334 rv.d *= tens[i]; 335 336 if (e1 &= ~15) 337 { 338 if (e1 > DBL_MAX_10_EXP) 339 { 340 ovfl: 341 ptr->_errno = ERANGE; 342 343 /* Force result to IEEE infinity. */ 344 word0 (rv) = Exp_mask; 345 word1 (rv) = 0; 346 347 if (bd0) 348 goto retfree; 349 goto ret; 350 } 351 if (e1 >>= 4) 352 { 353 for (j = 0; e1 > 1; j++, e1 >>= 1) 354 if (e1 & 1) 355 rv.d *= bigtens[j]; 356 /* The last multiplication could overflow. */ 357 word0 (rv) -= P * Exp_msk1; 358 rv.d *= bigtens[j]; 359 if ((z = word0 (rv) & Exp_mask) 360 > Exp_msk1 * (DBL_MAX_EXP + Bias - P)) 361 goto ovfl; 362 if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) 363 { 364 /* set to largest number */ 365 /* (Can't trust DBL_MAX) */ 366 word0 (rv) = Big0; 367 #ifndef _DOUBLE_IS_32BITS 368 word1 (rv) = Big1; 369 #endif 370 } 371 else 372 word0 (rv) += P * Exp_msk1; 373 } 374 375 } 376 } 377 else if (e1 < 0) 378 { 379 e1 = -e1; 380 if ((i = e1 & 15)) 381 rv.d /= tens[i]; 382 if (e1 &= ~15) 383 { 384 e1 >>= 4; 385 if (e1 >= 1 << n_bigtens) 386 goto undfl; 387 for (j = 0; e1 > 1; j++, e1 >>= 1) 388 if (e1 & 1) 389 rv.d *= tinytens[j]; 390 /* The last multiplication could underflow. */ 391 rv0.d = rv.d; 392 rv.d *= tinytens[j]; 393 if (!rv.d) 394 { 395 rv.d = 2. * rv0.d; 396 rv.d *= tinytens[j]; 397 if (!rv.d) 398 { 399 undfl: 400 rv.d = 0.; 401 ptr->_errno = ERANGE; 402 if (bd0) 403 goto retfree; 404 goto ret; 405 } 406 #ifndef _DOUBLE_IS_32BITS 407 word0 (rv) = Tiny0; 408 word1 (rv) = Tiny1; 409 #else 410 word0 (rv) = Tiny1; 411 #endif 412 /* The refinement below will clean 413 * this approximation up. 414 */ 415 } 416 } 417 } 418 419 /* Now the hard part -- adjusting rv to the correct value.*/ 420 421 /* Put digits into bd: true value = bd * 10^e */ 422 423 bd0 = s2b (ptr, s0, nd0, nd, y); 424 425 for (;;) 426 { 427 bd = Balloc (ptr, bd0->_k); 428 Bcopy (bd, bd0); 429 bb = d2b (ptr, rv.d, &bbe, &bbbits); /* rv.d = bb * 2^bbe */ 430 bs = i2b (ptr, 1); 431 432 if (e >= 0) 433 { 434 bb2 = bb5 = 0; 435 bd2 = bd5 = e; 436 } 437 else 438 { 439 bb2 = bb5 = -e; 440 bd2 = bd5 = 0; 441 } 442 if (bbe >= 0) 443 bb2 += bbe; 444 else 445 bd2 -= bbe; 446 bs2 = bb2; 447 #ifdef Sudden_Underflow 448 #ifdef IBM 449 j = 1 + 4 * P - 3 - bbbits + ((bbe + bbbits - 1) & 3); 450 #else 451 j = P + 1 - bbbits; 452 #endif 453 #else 454 i = bbe + bbbits - 1; /* logb(rv.d) */ 455 if (i < Emin) /* denormal */ 456 j = bbe + (P - Emin); 457 else 458 j = P + 1 - bbbits; 459 #endif 460 bb2 += j; 461 bd2 += j; 462 i = bb2 < bd2 ? bb2 : bd2; 463 if (i > bs2) 464 i = bs2; 465 if (i > 0) 466 { 467 bb2 -= i; 468 bd2 -= i; 469 bs2 -= i; 470 } 471 if (bb5 > 0) 472 { 473 bs = pow5mult (ptr, bs, bb5); 474 bb1 = mult (ptr, bs, bb); 475 Bfree (ptr, bb); 476 bb = bb1; 477 } 478 if (bb2 > 0) 479 bb = lshift (ptr, bb, bb2); 480 if (bd5 > 0) 481 bd = pow5mult (ptr, bd, bd5); 482 if (bd2 > 0) 483 bd = lshift (ptr, bd, bd2); 484 if (bs2 > 0) 485 bs = lshift (ptr, bs, bs2); 486 delta = diff (ptr, bb, bd); 487 dsign = delta->_sign; 488 delta->_sign = 0; 489 i = cmp (delta, bs); 490 if (i < 0) 491 { 492 /* Error is less than half an ulp -- check for 493 * special case of mantissa a power of two. 494 */ 495 if (dsign || word1 (rv) || word0 (rv) & Bndry_mask) 496 break; 497 delta = lshift (ptr, delta, Log2P); 498 if (cmp (delta, bs) > 0) 499 goto drop_down; 500 break; 501 } 502 if (i == 0) 503 { 504 /* exactly half-way between */ 505 if (dsign) 506 { 507 if ((word0 (rv) & Bndry_mask1) == Bndry_mask1 508 && word1 (rv) == 0xffffffff) 509 { 510 /*boundary case -- increment exponent*/ 511 word0 (rv) = (word0 (rv) & Exp_mask) 512 + Exp_msk1 513 #ifdef IBM 514 | Exp_msk1 >> 4 515 #endif 516 ; 517 #ifndef _DOUBLE_IS_32BITS 518 word1 (rv) = 0; 519 #endif 520 break; 521 } 522 } 523 else if (!(word0 (rv) & Bndry_mask) && !word1 (rv)) 524 { 525 drop_down: 526 /* boundary case -- decrement exponent */ 527 #ifdef Sudden_Underflow 528 L = word0 (rv) & Exp_mask; 529 #ifdef IBM 530 if (L < Exp_msk1) 531 #else 532 if (L <= Exp_msk1) 533 #endif 534 goto undfl; 535 L -= Exp_msk1; 536 #else 537 L = (word0 (rv) & Exp_mask) - Exp_msk1; 538 #endif 539 word0 (rv) = L | Bndry_mask1; 540 #ifndef _DOUBLE_IS_32BITS 541 word1 (rv) = 0xffffffff; 542 #endif 543 #ifdef IBM 544 goto cont; 545 #else 546 break; 547 #endif 548 } 549 #ifndef ROUND_BIASED 550 if (!(word1 (rv) & LSB)) 551 break; 552 #endif 553 if (dsign) 554 rv.d += ulp (rv.d); 555 #ifndef ROUND_BIASED 556 else 557 { 558 rv.d -= ulp (rv.d); 559 #ifndef Sudden_Underflow 560 if (!rv.d) 561 goto undfl; 562 #endif 563 } 564 #endif 565 break; 566 } 567 if ((aadj = ratio (delta, bs)) <= 2.) 568 { 569 if (dsign) 570 aadj = aadj1 = 1.; 571 else if (word1 (rv) || word0 (rv) & Bndry_mask) 572 { 573 #ifndef Sudden_Underflow 574 if (word1 (rv) == Tiny1 && !word0 (rv)) 575 goto undfl; 576 #endif 577 aadj = 1.; 578 aadj1 = -1.; 579 } 580 else 581 { 582 /* special case -- power of FLT_RADIX to be */ 583 /* rounded down... */ 584 585 if (aadj < 2. / FLT_RADIX) 586 aadj = 1. / FLT_RADIX; 587 else 588 aadj *= 0.5; 589 aadj1 = -aadj; 590 } 591 } 592 else 593 { 594 aadj *= 0.5; 595 aadj1 = dsign ? aadj : -aadj; 596 #ifdef Check_FLT_ROUNDS 597 switch (FLT_ROUNDS) 598 { 599 case 2: /* towards +infinity */ 600 aadj1 -= 0.5; 601 break; 602 case 0: /* towards 0 */ 603 case 3: /* towards -infinity */ 604 aadj1 += 0.5; 605 } 606 #else 607 if (FLT_ROUNDS == 0) 608 aadj1 += 0.5; 609 #endif 610 } 611 y = word0 (rv) & Exp_mask; 612 613 /* Check for overflow */ 614 615 if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) 616 { 617 rv0.d = rv.d; 618 word0 (rv) -= P * Exp_msk1; 619 adj = aadj1 * ulp (rv.d); 620 rv.d += adj; 621 if ((word0 (rv) & Exp_mask) >= 622 Exp_msk1 * (DBL_MAX_EXP + Bias - P)) 623 { 624 if (word0 (rv0) == Big0 && word1 (rv0) == Big1) 625 goto ovfl; 626 #ifdef _DOUBLE_IS_32BITS 627 word0 (rv) = Big1; 628 #else 629 word0 (rv) = Big0; 630 word1 (rv) = Big1; 631 #endif 632 goto cont; 633 } 634 else 635 word0 (rv) += P * Exp_msk1; 636 } 637 else 638 { 639 #ifdef Sudden_Underflow 640 if ((word0 (rv) & Exp_mask) <= P * Exp_msk1) 641 { 642 rv0.d = rv.d; 643 word0 (rv) += P * Exp_msk1; 644 adj = aadj1 * ulp (rv.d); 645 rv.d += adj; 646 #ifdef IBM 647 if ((word0 (rv) & Exp_mask) < P * Exp_msk1) 648 #else 649 if ((word0 (rv) & Exp_mask) <= P * Exp_msk1) 650 #endif 651 { 652 if (word0 (rv0) == Tiny0 653 && word1 (rv0) == Tiny1) 654 goto undfl; 655 word0 (rv) = Tiny0; 656 word1 (rv) = Tiny1; 657 goto cont; 658 } 659 else 660 word0 (rv) -= P * Exp_msk1; 661 } 662 else 663 { 664 adj = aadj1 * ulp (rv.d); 665 rv.d += adj; 666 } 667 #else 668 /* Compute adj so that the IEEE rounding rules will 669 * correctly round rv.d + adj in some half-way cases. 670 * If rv.d * ulp(rv.d) is denormalized (i.e., 671 * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid 672 * trouble from bits lost to denormalization; 673 * example: 1.2e-307 . 674 */ 675 if (y <= (P - 1) * Exp_msk1 && aadj >= 1.) 676 { 677 aadj1 = (double) (int) (aadj + 0.5); 678 if (!dsign) 679 aadj1 = -aadj1; 680 } 681 adj = aadj1 * ulp (rv.d); 682 rv.d += adj; 683 #endif 684 } 685 z = word0 (rv) & Exp_mask; 686 if (y == z) 687 { 688 /* Can we stop now? */ 689 L = aadj; 690 aadj -= L; 691 /* The tolerances below are conservative. */ 692 if (dsign || word1 (rv) || word0 (rv) & Bndry_mask) 693 { 694 if (aadj < .4999999 || aadj > .5000001) 695 break; 696 } 697 else if (aadj < .4999999 / FLT_RADIX) 698 break; 699 } 700 cont: 701 Bfree (ptr, bb); 702 Bfree (ptr, bd); 703 Bfree (ptr, bs); 704 Bfree (ptr, delta); 705 } 706 retfree: 707 Bfree (ptr, bb); 708 Bfree (ptr, bd); 709 Bfree (ptr, bs); 710 Bfree (ptr, bd0); 711 Bfree (ptr, delta); 712 ret: 713 if (se) 714 *se = (char *) s; 715 if (digits == 0) 716 ptr->_errno = EINVAL; 717 return sign ? -rv.d : rv.d; 718 } 719 720