3 <<strtod>>, <<strtodf>>---string to double or float
14 double strtod(const char *<[str]>, char **<[tail]>);
15 float strtodf(const char *<[str]>, char **<[tail]>);
17 double _strtod_r(void *<[reent]>,
18 const char *<[str]>, char **<[tail]>);
22 double strtod(<[str]>,<[tail]>)
26 float strtodf(<[str]>,<[tail]>)
30 double _strtod_r(<[reent]>,<[str]>,<[tail]>)
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.
53 This implementation returns the nearest machine number to the
54 input decimal string. Ties are broken by using the IEEE
57 The alternate function <<_strtod_r>> is a reentrant version.
58 The extra argument <[reent]> is a pointer to a reentrancy structure.
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.
68 Supporting OS subroutines required: <<close>>, <<fstat>>, <<isatty>>,
69 <<lseek>>, <<read>>, <<sbrk>>, <<write>>.
72 /****************************************************************
74 * The author of this software is David M. Gay.
76 * Copyright (c) 1991 by AT&T.
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.
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.
89 ***************************************************************/
91 /* Please send bug reports to
93 AT&T Bell Laboratories, Room 2C-463
95 Murray Hill, NJ 07974-2070
97 dmg@research.att.com or research!dmg
106 _DEFUN (_strtod_r, (ptr, s00, se),
107 struct _Jv_reent *ptr _AND
108 _CONST char *s00 _AND
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. */
115 _CONST char *s, *s0, *s1;
116 double aadj, aadj1, adj;
119 union double_union rv, rv0;
121 _Jv_Bigint *bb, *bb1, *bd, *bd0, *bs, *delta;
159 for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++)
163 y = 10 * y + c - '0';
165 z = 10 * z + c - '0';
173 for (; c == '0'; c = *++s)
178 if (c > '0' && c <= '9')
188 for (; c >= '0' && c <= '9'; c = *++s)
196 for (i = 1; i < nz; i++)
199 else if (nd <= DBL_DIG + 1)
203 else if (nd <= DBL_DIG + 1)
211 if (c == 'e' || c == 'E')
213 if (!nd && !nz && !nz0)
227 if (c >= '0' && c <= '9')
231 if (c > '0' && c <= '9')
235 while ((c = *++s) >= '0' && c <= '9')
236 e = 10 * e + c - '0';
238 /* Avoid confusion from exponents
239 * so large that e might overflow.
248 /* No exponent after an 'E' : that's an error. */
249 ptr->_errno = EINVAL;
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
270 k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1;
273 rv.d = tens[k - 9] * rv.d + z;
290 /* rv.d = */ rounded_product (rv.d, tens[e]);
295 if (e <= Ten_pmax + i)
297 /* A fancier test would sometimes let us do
298 * this for larger i values.
303 /* VAX exponent range is so narrow we must
304 * worry about overflow here...
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))
312 word0 (rv) += P * Exp_msk1;
314 /* rv.d = */ rounded_product (rv.d, tens[e]);
319 #ifndef Inaccurate_Divide
320 else if (e >= -Ten_pmax)
322 /* rv.d = */ rounded_quotient (rv.d, tens[-e]);
329 /* Get starting approximation = rv.d * 10**e1 */
338 if (e1 > DBL_MAX_10_EXP)
341 ptr->_errno = ERANGE;
343 /* Force result to IEEE infinity. */
344 word0 (rv) = Exp_mask;
353 for (j = 0; e1 > 1; j++, e1 >>= 1)
356 /* The last multiplication could overflow. */
357 word0 (rv) -= P * Exp_msk1;
359 if ((z = word0 (rv) & Exp_mask)
360 > Exp_msk1 * (DBL_MAX_EXP + Bias - P))
362 if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P))
364 /* set to largest number */
365 /* (Can't trust DBL_MAX) */
367 #ifndef _DOUBLE_IS_32BITS
372 word0 (rv) += P * Exp_msk1;
385 if (e1 >= 1 << n_bigtens)
387 for (j = 0; e1 > 1; j++, e1 >>= 1)
390 /* The last multiplication could underflow. */
401 ptr->_errno = ERANGE;
406 #ifndef _DOUBLE_IS_32BITS
412 /* The refinement below will clean
413 * this approximation up.
419 /* Now the hard part -- adjusting rv to the correct value.*/
421 /* Put digits into bd: true value = bd * 10^e */
423 bd0 = s2b (ptr, s0, nd0, nd, y);
427 bd = Balloc (ptr, bd0->_k);
429 bb = d2b (ptr, rv.d, &bbe, &bbbits); /* rv.d = bb * 2^bbe */
447 #ifdef Sudden_Underflow
449 j = 1 + 4 * P - 3 - bbbits + ((bbe + bbbits - 1) & 3);
454 i = bbe + bbbits - 1; /* logb(rv.d) */
455 if (i < Emin) /* denormal */
456 j = bbe + (P - Emin);
462 i = bb2 < bd2 ? bb2 : bd2;
473 bs = pow5mult (ptr, bs, bb5);
474 bb1 = mult (ptr, bs, bb);
479 bb = lshift (ptr, bb, bb2);
481 bd = pow5mult (ptr, bd, bd5);
483 bd = lshift (ptr, bd, bd2);
485 bs = lshift (ptr, bs, bs2);
486 delta = diff (ptr, bb, bd);
487 dsign = delta->_sign;
492 /* Error is less than half an ulp -- check for
493 * special case of mantissa a power of two.
495 if (dsign || word1 (rv) || word0 (rv) & Bndry_mask)
497 delta = lshift (ptr, delta, Log2P);
498 if (cmp (delta, bs) > 0)
504 /* exactly half-way between */
507 if ((word0 (rv) & Bndry_mask1) == Bndry_mask1
508 && word1 (rv) == 0xffffffff)
510 /*boundary case -- increment exponent*/
511 word0 (rv) = (word0 (rv) & Exp_mask)
517 #ifndef _DOUBLE_IS_32BITS
523 else if (!(word0 (rv) & Bndry_mask) && !word1 (rv))
526 /* boundary case -- decrement exponent */
527 #ifdef Sudden_Underflow
528 L = word0 (rv) & Exp_mask;
537 L = (word0 (rv) & Exp_mask) - Exp_msk1;
539 word0 (rv) = L | Bndry_mask1;
540 #ifndef _DOUBLE_IS_32BITS
541 word1 (rv) = 0xffffffff;
550 if (!(word1 (rv) & LSB))
559 #ifndef Sudden_Underflow
567 if ((aadj = ratio (delta, bs)) <= 2.)
571 else if (word1 (rv) || word0 (rv) & Bndry_mask)
573 #ifndef Sudden_Underflow
574 if (word1 (rv) == Tiny1 && !word0 (rv))
582 /* special case -- power of FLT_RADIX to be */
583 /* rounded down... */
585 if (aadj < 2. / FLT_RADIX)
586 aadj = 1. / FLT_RADIX;
595 aadj1 = dsign ? aadj : -aadj;
596 #ifdef Check_FLT_ROUNDS
599 case 2: /* towards +infinity */
602 case 0: /* towards 0 */
603 case 3: /* towards -infinity */
611 y = word0 (rv) & Exp_mask;
613 /* Check for overflow */
615 if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1))
618 word0 (rv) -= P * Exp_msk1;
619 adj = aadj1 * ulp (rv.d);
621 if ((word0 (rv) & Exp_mask) >=
622 Exp_msk1 * (DBL_MAX_EXP + Bias - P))
624 if (word0 (rv0) == Big0 && word1 (rv0) == Big1)
626 #ifdef _DOUBLE_IS_32BITS
635 word0 (rv) += P * Exp_msk1;
639 #ifdef Sudden_Underflow
640 if ((word0 (rv) & Exp_mask) <= P * Exp_msk1)
643 word0 (rv) += P * Exp_msk1;
644 adj = aadj1 * ulp (rv.d);
647 if ((word0 (rv) & Exp_mask) < P * Exp_msk1)
649 if ((word0 (rv) & Exp_mask) <= P * Exp_msk1)
652 if (word0 (rv0) == Tiny0
653 && word1 (rv0) == Tiny1)
660 word0 (rv) -= P * Exp_msk1;
664 adj = aadj1 * ulp (rv.d);
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 .
675 if (y <= (P - 1) * Exp_msk1 && aadj >= 1.)
677 aadj1 = (double) (int) (aadj + 0.5);
681 adj = aadj1 * ulp (rv.d);
685 z = word0 (rv) & Exp_mask;
688 /* Can we stop now? */
691 /* The tolerances below are conservative. */
692 if (dsign || word1 (rv) || word0 (rv) & Bndry_mask)
694 if (aadj < .4999999 || aadj > .5000001)
697 else if (aadj < .4999999 / FLT_RADIX)
716 ptr->_errno = EINVAL;
717 return sign ? -rv.d : rv.d;