/* FUNCTION <>, <>---string to double or float INDEX strtod INDEX _strtod_r INDEX strtodf ANSI_SYNOPSIS #include double strtod(const char *<[str]>, char **<[tail]>); float strtodf(const char *<[str]>, char **<[tail]>); double _strtod_r(void *<[reent]>, const char *<[str]>, char **<[tail]>); TRAD_SYNOPSIS #include double strtod(<[str]>,<[tail]>) char *<[str]>; char **<[tail]>; float strtodf(<[str]>,<[tail]>) char *<[str]>; char **<[tail]>; double _strtod_r(<[reent]>,<[str]>,<[tail]>) char *<[reent]>; char *<[str]>; char **<[tail]>; DESCRIPTION The function <> parses the character string <[str]>, producing a substring which can be converted to a double value. The substring converted is the longest initial subsequence of <[str]>, beginning with the first non-whitespace character, that has the format: .[+|-]<[digits]>[.][<[digits]>][(e|E)[+|-]<[digits]>] The substring contains no characters if <[str]> is empty, consists entirely of whitespace, or if the first non-whitespace character is something other than <<+>>, <<->>, <<.>>, or a digit. If the substring is empty, no conversion is done, and the value of <[str]> is stored in <<*<[tail]>>>. Otherwise, the substring is converted, and a pointer to the final string (which will contain at least the terminating null character of <[str]>) is stored in <<*<[tail]>>>. If you want no assignment to <<*<[tail]>>>, pass a null pointer as <[tail]>. <> is identical to <> except for its return type. This implementation returns the nearest machine number to the input decimal string. Ties are broken by using the IEEE round-even rule. The alternate function <<_strtod_r>> is a reentrant version. The extra argument <[reent]> is a pointer to a reentrancy structure. RETURNS <> returns the converted substring value, if any. If no conversion could be performed, 0 is returned. If the correct value is out of the range of representable values, plus or minus <> is returned, and <> is stored in errno. If the correct value would cause underflow, 0 is returned and <> is stored in errno. Supporting OS subroutines required: <>, <>, <>, <>, <>, <>, <>. */ /**************************************************************** * * The author of this software is David M. Gay. * * Copyright (c) 1991 by AT&T. * * Permission to use, copy, modify, and distribute this software for any * purpose without fee is hereby granted, provided that this entire notice * is included in all copies of any software which is or includes a copy * or modification of this software and in all copies of the supporting * documentation for such software. * * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED * WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR AT&T MAKES ANY * REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE MERCHANTABILITY * OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR PURPOSE. * ***************************************************************/ /* Please send bug reports to David M. Gay AT&T Bell Laboratories, Room 2C-463 600 Mountain Avenue Murray Hill, NJ 07974-2070 U.S.A. dmg@research.att.com or research!dmg */ #include #include #include #include "mprec.h" double _DEFUN (_strtod_r, (ptr, s00, se), struct _Jv_reent *ptr _AND _CONST char *s00 _AND char **se) { int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, dsign, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; int digits = 0; /* Number of digits found in fraction part. */ long e; _CONST char *s, *s0, *s1; double aadj, aadj1, adj; long L; unsigned long y, z; union double_union rv, rv0; _Jv_Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; sign = nz0 = nz = 0; rv.d = 0.; for (s = s00;; s++) switch (*s) { case '-': sign = 1; /* no break */ case '+': if (*++s) goto break2; /* no break */ case 0: s = s00; goto ret; case '\t': case '\n': case '\v': case '\f': case '\r': case ' ': continue; default: goto break2; } break2: if (*s == '0') { digits++; nz0 = 1; while (*++s == '0') digits++; if (!*s) goto ret; } s0 = s; y = z = 0; for (nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) { digits++; if (nd < 9) y = 10 * y + c - '0'; else if (nd < 16) z = 10 * z + c - '0'; } nd0 = nd; if (c == '.') { c = *++s; if (!nd) { for (; c == '0'; c = *++s) { digits++; nz++; } if (c > '0' && c <= '9') { digits++; s0 = s; nf += nz; nz = 0; goto have_dig; } goto dig_done; } for (; c >= '0' && c <= '9'; c = *++s) { digits++; have_dig: nz++; if (c -= '0') { nf += nz; for (i = 1; i < nz; i++) if (nd++ < 9) y *= 10; else if (nd <= DBL_DIG + 1) z *= 10; if (nd++ < 9) y = 10 * y + c; else if (nd <= DBL_DIG + 1) z = 10 * z + c; nz = 0; } } } dig_done: e = 0; if (c == 'e' || c == 'E') { if (!nd && !nz && !nz0) { s = s00; goto ret; } s00 = s; esign = 0; switch (c = *++s) { case '-': esign = 1; case '+': c = *++s; } if (c >= '0' && c <= '9') { while (c == '0') c = *++s; if (c > '0' && c <= '9') { e = c - '0'; s1 = s; while ((c = *++s) >= '0' && c <= '9') e = 10 * e + c - '0'; if (s - s1 > 8) /* Avoid confusion from exponents * so large that e might overflow. */ e = 9999999L; if (esign) e = -e; } } else { /* No exponent after an 'E' : that's an error. */ ptr->_errno = EINVAL; e = 0; s = s00; goto ret; } } if (!nd) { if (!nz && !nz0) s = s00; goto ret; } e1 = e -= nf; /* Now we have nd0 digits, starting at s0, followed by a * decimal point, followed by nd-nd0 digits. The number we're * after is the integer represented by those digits times * 10**e */ if (!nd0) nd0 = nd; k = nd < DBL_DIG + 1 ? nd : DBL_DIG + 1; rv.d = y; if (k > 9) rv.d = tens[k - 9] * rv.d + z; bd0 = 0; if (nd <= DBL_DIG #ifndef RND_PRODQUOT && FLT_ROUNDS == 1 #endif ) { if (!e) goto ret; if (e > 0) { if (e <= Ten_pmax) { #ifdef VAX goto vax_ovfl_check; #else /* rv.d = */ rounded_product (rv.d, tens[e]); goto ret; #endif } i = DBL_DIG - nd; if (e <= Ten_pmax + i) { /* A fancier test would sometimes let us do * this for larger i values. */ e -= i; rv.d *= tens[i]; #ifdef VAX /* VAX exponent range is so narrow we must * worry about overflow here... */ vax_ovfl_check: word0 (rv) -= P * Exp_msk1; /* rv.d = */ rounded_product (rv.d, tens[e]); if ((word0 (rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) goto ovfl; word0 (rv) += P * Exp_msk1; #else /* rv.d = */ rounded_product (rv.d, tens[e]); #endif goto ret; } } #ifndef Inaccurate_Divide else if (e >= -Ten_pmax) { /* rv.d = */ rounded_quotient (rv.d, tens[-e]); goto ret; } #endif } e1 += nd - k; /* Get starting approximation = rv.d * 10**e1 */ if (e1 > 0) { if ((i = e1 & 15)) rv.d *= tens[i]; if (e1 &= ~15) { if (e1 > DBL_MAX_10_EXP) { ovfl: ptr->_errno = ERANGE; /* Force result to IEEE infinity. */ word0 (rv) = Exp_mask; word1 (rv) = 0; if (bd0) goto retfree; goto ret; } if (e1 >>= 4) { for (j = 0; e1 > 1; j++, e1 >>= 1) if (e1 & 1) rv.d *= bigtens[j]; /* The last multiplication could overflow. */ word0 (rv) -= P * Exp_msk1; rv.d *= bigtens[j]; if ((z = word0 (rv) & Exp_mask) > Exp_msk1 * (DBL_MAX_EXP + Bias - P)) goto ovfl; if (z > Exp_msk1 * (DBL_MAX_EXP + Bias - 1 - P)) { /* set to largest number */ /* (Can't trust DBL_MAX) */ word0 (rv) = Big0; #ifndef _DOUBLE_IS_32BITS word1 (rv) = Big1; #endif } else word0 (rv) += P * Exp_msk1; } } } else if (e1 < 0) { e1 = -e1; if ((i = e1 & 15)) rv.d /= tens[i]; if (e1 &= ~15) { e1 >>= 4; if (e1 >= 1 << n_bigtens) goto undfl; for (j = 0; e1 > 1; j++, e1 >>= 1) if (e1 & 1) rv.d *= tinytens[j]; /* The last multiplication could underflow. */ rv0.d = rv.d; rv.d *= tinytens[j]; if (!rv.d) { rv.d = 2. * rv0.d; rv.d *= tinytens[j]; if (!rv.d) { undfl: rv.d = 0.; ptr->_errno = ERANGE; if (bd0) goto retfree; goto ret; } #ifndef _DOUBLE_IS_32BITS word0 (rv) = Tiny0; word1 (rv) = Tiny1; #else word0 (rv) = Tiny1; #endif /* The refinement below will clean * this approximation up. */ } } } /* Now the hard part -- adjusting rv to the correct value.*/ /* Put digits into bd: true value = bd * 10^e */ bd0 = s2b (ptr, s0, nd0, nd, y); for (;;) { bd = Balloc (ptr, bd0->_k); Bcopy (bd, bd0); bb = d2b (ptr, rv.d, &bbe, &bbbits); /* rv.d = bb * 2^bbe */ bs = i2b (ptr, 1); if (e >= 0) { bb2 = bb5 = 0; bd2 = bd5 = e; } else { bb2 = bb5 = -e; bd2 = bd5 = 0; } if (bbe >= 0) bb2 += bbe; else bd2 -= bbe; bs2 = bb2; #ifdef Sudden_Underflow #ifdef IBM j = 1 + 4 * P - 3 - bbbits + ((bbe + bbbits - 1) & 3); #else j = P + 1 - bbbits; #endif #else i = bbe + bbbits - 1; /* logb(rv.d) */ if (i < Emin) /* denormal */ j = bbe + (P - Emin); else j = P + 1 - bbbits; #endif bb2 += j; bd2 += j; i = bb2 < bd2 ? bb2 : bd2; if (i > bs2) i = bs2; if (i > 0) { bb2 -= i; bd2 -= i; bs2 -= i; } if (bb5 > 0) { bs = pow5mult (ptr, bs, bb5); bb1 = mult (ptr, bs, bb); Bfree (ptr, bb); bb = bb1; } if (bb2 > 0) bb = lshift (ptr, bb, bb2); if (bd5 > 0) bd = pow5mult (ptr, bd, bd5); if (bd2 > 0) bd = lshift (ptr, bd, bd2); if (bs2 > 0) bs = lshift (ptr, bs, bs2); delta = diff (ptr, bb, bd); dsign = delta->_sign; delta->_sign = 0; i = cmp (delta, bs); if (i < 0) { /* Error is less than half an ulp -- check for * special case of mantissa a power of two. */ if (dsign || word1 (rv) || word0 (rv) & Bndry_mask) break; delta = lshift (ptr, delta, Log2P); if (cmp (delta, bs) > 0) goto drop_down; break; } if (i == 0) { /* exactly half-way between */ if (dsign) { if ((word0 (rv) & Bndry_mask1) == Bndry_mask1 && word1 (rv) == 0xffffffff) { /*boundary case -- increment exponent*/ word0 (rv) = (word0 (rv) & Exp_mask) + Exp_msk1 #ifdef IBM | Exp_msk1 >> 4 #endif ; #ifndef _DOUBLE_IS_32BITS word1 (rv) = 0; #endif break; } } else if (!(word0 (rv) & Bndry_mask) && !word1 (rv)) { drop_down: /* boundary case -- decrement exponent */ #ifdef Sudden_Underflow L = word0 (rv) & Exp_mask; #ifdef IBM if (L < Exp_msk1) #else if (L <= Exp_msk1) #endif goto undfl; L -= Exp_msk1; #else L = (word0 (rv) & Exp_mask) - Exp_msk1; #endif word0 (rv) = L | Bndry_mask1; #ifndef _DOUBLE_IS_32BITS word1 (rv) = 0xffffffff; #endif #ifdef IBM goto cont; #else break; #endif } #ifndef ROUND_BIASED if (!(word1 (rv) & LSB)) break; #endif if (dsign) rv.d += ulp (rv.d); #ifndef ROUND_BIASED else { rv.d -= ulp (rv.d); #ifndef Sudden_Underflow if (!rv.d) goto undfl; #endif } #endif break; } if ((aadj = ratio (delta, bs)) <= 2.) { if (dsign) aadj = aadj1 = 1.; else if (word1 (rv) || word0 (rv) & Bndry_mask) { #ifndef Sudden_Underflow if (word1 (rv) == Tiny1 && !word0 (rv)) goto undfl; #endif aadj = 1.; aadj1 = -1.; } else { /* special case -- power of FLT_RADIX to be */ /* rounded down... */ if (aadj < 2. / FLT_RADIX) aadj = 1. / FLT_RADIX; else aadj *= 0.5; aadj1 = -aadj; } } else { aadj *= 0.5; aadj1 = dsign ? aadj : -aadj; #ifdef Check_FLT_ROUNDS switch (FLT_ROUNDS) { case 2: /* towards +infinity */ aadj1 -= 0.5; break; case 0: /* towards 0 */ case 3: /* towards -infinity */ aadj1 += 0.5; } #else if (FLT_ROUNDS == 0) aadj1 += 0.5; #endif } y = word0 (rv) & Exp_mask; /* Check for overflow */ if (y == Exp_msk1 * (DBL_MAX_EXP + Bias - 1)) { rv0.d = rv.d; word0 (rv) -= P * Exp_msk1; adj = aadj1 * ulp (rv.d); rv.d += adj; if ((word0 (rv) & Exp_mask) >= Exp_msk1 * (DBL_MAX_EXP + Bias - P)) { if (word0 (rv0) == Big0 && word1 (rv0) == Big1) goto ovfl; #ifdef _DOUBLE_IS_32BITS word0 (rv) = Big1; #else word0 (rv) = Big0; word1 (rv) = Big1; #endif goto cont; } else word0 (rv) += P * Exp_msk1; } else { #ifdef Sudden_Underflow if ((word0 (rv) & Exp_mask) <= P * Exp_msk1) { rv0.d = rv.d; word0 (rv) += P * Exp_msk1; adj = aadj1 * ulp (rv.d); rv.d += adj; #ifdef IBM if ((word0 (rv) & Exp_mask) < P * Exp_msk1) #else if ((word0 (rv) & Exp_mask) <= P * Exp_msk1) #endif { if (word0 (rv0) == Tiny0 && word1 (rv0) == Tiny1) goto undfl; word0 (rv) = Tiny0; word1 (rv) = Tiny1; goto cont; } else word0 (rv) -= P * Exp_msk1; } else { adj = aadj1 * ulp (rv.d); rv.d += adj; } #else /* Compute adj so that the IEEE rounding rules will * correctly round rv.d + adj in some half-way cases. * If rv.d * ulp(rv.d) is denormalized (i.e., * y <= (P-1)*Exp_msk1), we must adjust aadj to avoid * trouble from bits lost to denormalization; * example: 1.2e-307 . */ if (y <= (P - 1) * Exp_msk1 && aadj >= 1.) { aadj1 = (double) (int) (aadj + 0.5); if (!dsign) aadj1 = -aadj1; } adj = aadj1 * ulp (rv.d); rv.d += adj; #endif } z = word0 (rv) & Exp_mask; if (y == z) { /* Can we stop now? */ L = aadj; aadj -= L; /* The tolerances below are conservative. */ if (dsign || word1 (rv) || word0 (rv) & Bndry_mask) { if (aadj < .4999999 || aadj > .5000001) break; } else if (aadj < .4999999 / FLT_RADIX) break; } cont: Bfree (ptr, bb); Bfree (ptr, bd); Bfree (ptr, bs); Bfree (ptr, delta); } retfree: Bfree (ptr, bb); Bfree (ptr, bd); Bfree (ptr, bs); Bfree (ptr, bd0); Bfree (ptr, delta); ret: if (se) *se = (char *) s; if (digits == 0) ptr->_errno = EINVAL; return sign ? -rv.d : rv.d; }