0,0 → 1,1197 |
/* |
FUNCTION |
<<strtod>>, <<strtof>>---string to double or float |
|
INDEX |
strtod |
INDEX |
_strtod_r |
INDEX |
strtof |
|
ANSI_SYNOPSIS |
#include <stdlib.h> |
double strtod(const char *<[str]>, char **<[tail]>); |
float strtof(const char *<[str]>, char **<[tail]>); |
|
double _strtod_r(void *<[reent]>, |
const char *<[str]>, char **<[tail]>); |
|
TRAD_SYNOPSIS |
#include <stdlib.h> |
double strtod(<[str]>,<[tail]>) |
char *<[str]>; |
char **<[tail]>; |
|
float strtof(<[str]>,<[tail]>) |
char *<[str]>; |
char **<[tail]>; |
|
double _strtod_r(<[reent]>,<[str]>,<[tail]>) |
char *<[reent]>; |
char *<[str]>; |
char **<[tail]>; |
|
DESCRIPTION |
The function <<strtod>> 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 one of these formats: |
.[+|-]<[digits]>[.[<[digits]>]][(e|E)[+|-]<[digits]>] |
.[+|-].<[digits]>[(e|E)[+|-]<[digits]>] |
.[+|-](i|I)(n|N)(f|F)[(i|I)(n|N)(i|I)(t|T)(y|Y)] |
.[+|-](n|N)(a|A)(n|N)[<(>[<[hexdigits]>]<)>] |
.[+|-]0(x|X)<[hexdigits]>[.[<[hexdigits]>]][(p|P)[+|-]<[digits]>] |
.[+|-]0(x|X).<[hexdigits]>[(p|P)[+|-]<[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, and cannot be parsed as infinity or NaN. If the platform |
does not support NaN, then NaN is treated as an empty substring. |
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]>. |
<<strtof>> is identical to <<strtod>> 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. However, <<strtof>> is currently subject to |
double rounding errors. |
|
The alternate function <<_strtod_r>> is a reentrant version. |
The extra argument <[reent]> is a pointer to a reentrancy structure. |
|
RETURNS |
<<strtod>> 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 <<HUGE_VAL>> is returned, and <<ERANGE>> is |
stored in errno. If the correct value would cause underflow, 0 |
is returned and <<ERANGE>> is stored in errno. |
|
Supporting OS subroutines required: <<close>>, <<fstat>>, <<isatty>>, |
<<lseek>>, <<read>>, <<sbrk>>, <<write>>. |
*/ |
|
/**************************************************************** |
|
The author of this software is David M. Gay. |
|
Copyright (C) 1998-2001 by Lucent Technologies |
All Rights Reserved |
|
Permission to use, copy, modify, and distribute this software and |
its documentation for any purpose and without fee is hereby |
granted, provided that the above copyright notice appear in all |
copies and that both that the copyright notice and this |
permission notice and warranty disclaimer appear in supporting |
documentation, and that the name of Lucent or any of its entities |
not be used in advertising or publicity pertaining to |
distribution of the software without specific, written prior |
permission. |
|
LUCENT DISCLAIMS ALL WARRANTIES WITH REGARD TO THIS SOFTWARE, |
INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS. |
IN NO EVENT SHALL LUCENT OR ANY OF ITS ENTITIES BE LIABLE FOR ANY |
SPECIAL, INDIRECT OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES |
WHATSOEVER RESULTING FROM LOSS OF USE, DATA OR PROFITS, WHETHER |
IN AN ACTION OF CONTRACT, NEGLIGENCE OR OTHER TORTIOUS ACTION, |
ARISING OUT OF OR IN CONNECTION WITH THE USE OR PERFORMANCE OF |
THIS SOFTWARE. |
|
****************************************************************/ |
|
/* Please send bug reports to David M. Gay (dmg at acm dot org, |
* with " at " changed at "@" and " dot " changed to "."). */ |
|
/* Original file gdtoa-strtod.c Modified 06-21-2006 by Jeff Johnston to work within newlib. */ |
|
#include <_ansi.h> |
#include <errno.h> |
#include <stdlib.h> |
#include <string.h> |
#include "mprec.h" |
#include "gdtoa.h" |
#include "gd_qnan.h" |
|
/* #ifndef NO_FENV_H */ |
/* #include <fenv.h> */ |
/* #endif */ |
|
#include "locale.h" |
|
#ifdef IEEE_Arith |
#ifndef NO_IEEE_Scale |
#define Avoid_Underflow |
#undef tinytens |
/* The factor of 2^53 in tinytens[4] helps us avoid setting the underflow */ |
/* flag unnecessarily. It leads to a song and dance at the end of strtod. */ |
static _CONST double tinytens[] = { 1e-16, 1e-32, 1e-64, 1e-128, |
9007199254740992.e-256 |
}; |
#endif |
#endif |
|
#ifdef Honor_FLT_ROUNDS |
#define Rounding rounding |
#undef Check_FLT_ROUNDS |
#define Check_FLT_ROUNDS |
#else |
#define Rounding Flt_Rounds |
#endif |
|
#ifndef NO_HEX_FP |
|
static void |
_DEFUN (ULtod, (L, bits, exp, k), |
__ULong *L _AND |
__ULong *bits _AND |
Long exp _AND |
int k) |
{ |
switch(k & STRTOG_Retmask) { |
case STRTOG_NoNumber: |
case STRTOG_Zero: |
L[0] = L[1] = 0; |
break; |
|
case STRTOG_Denormal: |
L[_1] = bits[0]; |
L[_0] = bits[1]; |
break; |
|
case STRTOG_Normal: |
case STRTOG_NaNbits: |
L[_1] = bits[0]; |
L[_0] = (bits[1] & ~0x100000) | ((exp + 0x3ff + 52) << 20); |
break; |
|
case STRTOG_Infinite: |
L[_0] = 0x7ff00000; |
L[_1] = 0; |
break; |
|
case STRTOG_NaN: |
L[_0] = 0x7fffffff; |
L[_1] = (__ULong)-1; |
} |
if (k & STRTOG_Neg) |
L[_0] |= 0x80000000L; |
} |
#endif /* !NO_HEX_FP */ |
|
#ifdef INFNAN_CHECK |
static int |
_DEFUN (match, (sp, t), |
_CONST char **sp _AND |
char *t) |
{ |
int c, d; |
_CONST char *s = *sp; |
|
while( (d = *t++) !=0) { |
if ((c = *++s) >= 'A' && c <= 'Z') |
c += 'a' - 'A'; |
if (c != d) |
return 0; |
} |
*sp = s + 1; |
return 1; |
} |
#endif /* INFNAN_CHECK */ |
|
|
double |
_DEFUN (_strtod_r, (ptr, s00, se), |
struct _reent *ptr _AND |
_CONST char *s00 _AND |
char **se) |
{ |
#ifdef Avoid_Underflow |
int scale; |
#endif |
int bb2, bb5, bbe, bd2, bd5, bbbits, bs2, c, decpt, dsign, |
e, e1, esign, i, j, k, nd, nd0, nf, nz, nz0, sign; |
_CONST char *s, *s0, *s1; |
double aadj, adj; |
U aadj1, rv, rv0; |
Long L; |
__ULong y, z; |
_Bigint *bb, *bb1, *bd, *bd0, *bs, *delta; |
#ifdef SET_INEXACT |
int inexact, oldinexact; |
#endif |
#ifdef Honor_FLT_ROUNDS |
int rounding; |
#endif |
|
delta = bs = bd = NULL; |
sign = nz0 = nz = decpt = 0; |
dval(rv) = 0.; |
for(s = s00;;s++) switch(*s) { |
case '-': |
sign = 1; |
/* no break */ |
case '+': |
if (*++s) |
goto break2; |
/* no break */ |
case 0: |
goto ret0; |
case '\t': |
case '\n': |
case '\v': |
case '\f': |
case '\r': |
case ' ': |
continue; |
default: |
goto break2; |
} |
break2: |
if (*s == '0') { |
#ifndef NO_HEX_FP |
{ |
static FPI fpi = { 53, 1-1023-53+1, 2046-1023-53+1, 1, SI }; |
Long exp; |
__ULong bits[2]; |
switch(s[1]) { |
case 'x': |
case 'X': |
/* If the number is not hex, then the parse of |
0 is still valid. */ |
s00 = s + 1; |
{ |
#if defined(FE_DOWNWARD) && defined(FE_TONEAREST) && defined(FE_TOWARDZERO) && defined(FE_UPWARD) |
FPI fpi1 = fpi; |
switch(fegetround()) { |
case FE_TOWARDZERO: fpi1.rounding = 0; break; |
case FE_UPWARD: fpi1.rounding = 2; break; |
case FE_DOWNWARD: fpi1.rounding = 3; |
} |
#else |
#define fpi1 fpi |
#endif |
switch((i = gethex(ptr, &s, &fpi1, &exp, &bb, sign)) & STRTOG_Retmask) { |
case STRTOG_NoNumber: |
s = s00; |
case STRTOG_Zero: |
break; |
default: |
if (bb) { |
copybits(bits, fpi.nbits, bb); |
Bfree(ptr,bb); |
} |
ULtod(rv.i, bits, exp, i); |
}} |
goto ret; |
} |
} |
#endif |
nz0 = 1; |
while(*++s == '0') ; |
if (!*s) |
goto ret; |
} |
s0 = s; |
y = z = 0; |
for(nd = nf = 0; (c = *s) >= '0' && c <= '9'; nd++, s++) { |
if (nd < DBL_DIG + 1) { |
if (nd < 9) |
y = 10*y + c - '0'; |
else |
z = 10*z + c - '0'; |
} |
} |
nd0 = nd; |
if (strncmp (s, _localeconv_r (ptr)->decimal_point, |
strlen (_localeconv_r (ptr)->decimal_point)) == 0) { |
decpt = 1; |
c = *(s += strlen (_localeconv_r (ptr)->decimal_point)); |
if (!nd) { |
for(; c == '0'; c = *++s) |
nz++; |
if (c > '0' && c <= '9') { |
s0 = s; |
nf += nz; |
nz = 0; |
goto have_dig; |
} |
goto dig_done; |
} |
for(; c >= '0' && c <= '9'; c = *++s) { |
have_dig: |
nz++; |
if (c -= '0') { |
for(i = 1; i < nz; i++) { |
if (nd <= DBL_DIG + 1) { |
if (nd + i < 10) |
y *= 10; |
else |
z *= 10; |
} |
} |
if (nd <= DBL_DIG + 1) { |
if (nd + i < 10) |
y = 10*y + c; |
else |
z = 10*z + c; |
} |
if (nd <= DBL_DIG + 1) { |
nf += nz; |
nd += nz; |
} |
nz = 0; |
} |
} |
} |
dig_done: |
e = 0; |
if (c == 'e' || c == 'E') { |
if (!nd && !nz && !nz0) { |
goto ret0; |
} |
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') { |
L = c - '0'; |
s1 = s; |
while((c = *++s) >= '0' && c <= '9') |
L = 10*L + c - '0'; |
if (s - s1 > 8 || L > 19999) |
/* Avoid confusion from exponents |
* so large that e might overflow. |
*/ |
e = 19999; /* safe for 16 bit ints */ |
else |
e = (int)L; |
if (esign) |
e = -e; |
} |
else |
e = 0; |
} |
else |
s = s00; |
} |
if (!nd) { |
if (!nz && !nz0) { |
#ifdef INFNAN_CHECK |
/* Check for Nan and Infinity */ |
__ULong bits[2]; |
static FPI fpinan = /* only 52 explicit bits */ |
{ 52, 1-1023-53+1, 2046-1023-53+1, 1, SI }; |
if (!decpt) |
switch(c) { |
case 'i': |
case 'I': |
if (match(&s,"nf")) { |
--s; |
if (!match(&s,"inity")) |
++s; |
dword0(rv) = 0x7ff00000; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = 0; |
#endif /*!_DOUBLE_IS_32BITS*/ |
goto ret; |
} |
break; |
case 'n': |
case 'N': |
if (match(&s, "an")) { |
#ifndef No_Hex_NaN |
if (*s == '(' /*)*/ |
&& hexnan(&s, &fpinan, bits) |
== STRTOG_NaNbits) { |
dword0(rv) = 0x7ff00000 | bits[1]; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = bits[0]; |
#endif /*!_DOUBLE_IS_32BITS*/ |
} |
else { |
#endif |
dword0(rv) = NAN_WORD0; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = NAN_WORD1; |
#endif /*!_DOUBLE_IS_32BITS*/ |
#ifndef No_Hex_NaN |
} |
#endif |
goto ret; |
} |
} |
#endif /* INFNAN_CHECK */ |
ret0: |
s = s00; |
sign = 0; |
} |
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; |
dval(rv) = y; |
if (k > 9) { |
#ifdef SET_INEXACT |
if (k > DBL_DIG) |
oldinexact = get_inexact(); |
#endif |
dval(rv) = tens[k - 9] * dval(rv) + z; |
} |
bd0 = 0; |
if (nd <= DBL_DIG |
#ifndef RND_PRODQUOT |
#ifndef Honor_FLT_ROUNDS |
&& Flt_Rounds == 1 |
#endif |
#endif |
) { |
if (!e) |
goto ret; |
if (e > 0) { |
if (e <= Ten_pmax) { |
#ifdef VAX |
goto vax_ovfl_check; |
#else |
#ifdef Honor_FLT_ROUNDS |
/* round correctly FLT_ROUNDS = 2 or 3 */ |
if (sign) { |
dval(rv) = -dval(rv); |
sign = 0; |
} |
#endif |
/* rv = */ rounded_product(dval(rv), 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. |
*/ |
#ifdef Honor_FLT_ROUNDS |
/* round correctly FLT_ROUNDS = 2 or 3 */ |
if (sign) { |
dval(rv) = -dval(rv); |
sign = 0; |
} |
#endif |
e -= i; |
dval(rv) *= tens[i]; |
#ifdef VAX |
/* VAX exponent range is so narrow we must |
* worry about overflow here... |
*/ |
vax_ovfl_check: |
dword0(rv) -= P*Exp_msk1; |
/* rv = */ rounded_product(dval(rv), tens[e]); |
if ((dword0(rv) & Exp_mask) |
> Exp_msk1*(DBL_MAX_EXP+Bias-1-P)) |
goto ovfl; |
dword0(rv) += P*Exp_msk1; |
#else |
/* rv = */ rounded_product(dval(rv), tens[e]); |
#endif |
goto ret; |
} |
} |
#ifndef Inaccurate_Divide |
else if (e >= -Ten_pmax) { |
#ifdef Honor_FLT_ROUNDS |
/* round correctly FLT_ROUNDS = 2 or 3 */ |
if (sign) { |
dval(rv) = -dval(rv); |
sign = 0; |
} |
#endif |
/* rv = */ rounded_quotient(dval(rv), tens[-e]); |
goto ret; |
} |
#endif |
} |
e1 += nd - k; |
|
#ifdef IEEE_Arith |
#ifdef SET_INEXACT |
inexact = 1; |
if (k <= DBL_DIG) |
oldinexact = get_inexact(); |
#endif |
#ifdef Avoid_Underflow |
scale = 0; |
#endif |
#ifdef Honor_FLT_ROUNDS |
if ((rounding = Flt_Rounds) >= 2) { |
if (sign) |
rounding = rounding == 2 ? 0 : 2; |
else |
if (rounding != 2) |
rounding = 0; |
} |
#endif |
#endif /*IEEE_Arith*/ |
|
/* Get starting approximation = rv * 10**e1 */ |
|
if (e1 > 0) { |
if ( (i = e1 & 15) !=0) |
dval(rv) *= tens[i]; |
if (e1 &= ~15) { |
if (e1 > DBL_MAX_10_EXP) { |
ovfl: |
#ifndef NO_ERRNO |
ptr->_errno = ERANGE; |
#endif |
/* Can't trust HUGE_VAL */ |
#ifdef IEEE_Arith |
#ifdef Honor_FLT_ROUNDS |
switch(rounding) { |
case 0: /* toward 0 */ |
case 3: /* toward -infinity */ |
dword0(rv) = Big0; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = Big1; |
#endif /*!_DOUBLE_IS_32BITS*/ |
break; |
default: |
dword0(rv) = Exp_mask; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = 0; |
#endif /*!_DOUBLE_IS_32BITS*/ |
} |
#else /*Honor_FLT_ROUNDS*/ |
dword0(rv) = Exp_mask; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = 0; |
#endif /*!_DOUBLE_IS_32BITS*/ |
#endif /*Honor_FLT_ROUNDS*/ |
#ifdef SET_INEXACT |
/* set overflow bit */ |
dval(rv0) = 1e300; |
dval(rv0) *= dval(rv0); |
#endif |
#else /*IEEE_Arith*/ |
dword0(rv) = Big0; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = Big1; |
#endif /*!_DOUBLE_IS_32BITS*/ |
#endif /*IEEE_Arith*/ |
if (bd0) |
goto retfree; |
goto ret; |
} |
e1 >>= 4; |
for(j = 0; e1 > 1; j++, e1 >>= 1) |
if (e1 & 1) |
dval(rv) *= bigtens[j]; |
/* The last multiplication could overflow. */ |
dword0(rv) -= P*Exp_msk1; |
dval(rv) *= bigtens[j]; |
if ((z = dword0(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) */ |
dword0(rv) = Big0; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = Big1; |
#endif /*!_DOUBLE_IS_32BITS*/ |
} |
else |
dword0(rv) += P*Exp_msk1; |
} |
} |
else if (e1 < 0) { |
e1 = -e1; |
if ( (i = e1 & 15) !=0) |
dval(rv) /= tens[i]; |
if (e1 >>= 4) { |
if (e1 >= 1 << n_bigtens) |
goto undfl; |
#ifdef Avoid_Underflow |
if (e1 & Scale_Bit) |
scale = 2*P; |
for(j = 0; e1 > 0; j++, e1 >>= 1) |
if (e1 & 1) |
dval(rv) *= tinytens[j]; |
if (scale && (j = 2*P + 1 - ((dword0(rv) & Exp_mask) |
>> Exp_shift)) > 0) { |
/* scaled rv is denormal; zap j low bits */ |
if (j >= 32) { |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = 0; |
#endif /*!_DOUBLE_IS_32BITS*/ |
if (j >= 53) |
dword0(rv) = (P+2)*Exp_msk1; |
else |
dword0(rv) &= 0xffffffff << (j-32); |
} |
#ifndef _DOUBLE_IS_32BITS |
else |
dword1(rv) &= 0xffffffff << j; |
#endif /*!_DOUBLE_IS_32BITS*/ |
} |
#else |
for(j = 0; e1 > 1; j++, e1 >>= 1) |
if (e1 & 1) |
dval(rv) *= tinytens[j]; |
/* The last multiplication could underflow. */ |
dval(rv0) = dval(rv); |
dval(rv) *= tinytens[j]; |
if (!dval(rv)) { |
dval(rv) = 2.*dval(rv0); |
dval(rv) *= tinytens[j]; |
#endif |
if (!dval(rv)) { |
undfl: |
dval(rv) = 0.; |
#ifndef NO_ERRNO |
ptr->_errno = ERANGE; |
#endif |
if (bd0) |
goto retfree; |
goto ret; |
} |
#ifndef Avoid_Underflow |
#ifndef _DOUBLE_IS_32BITS |
dword0(rv) = Tiny0; |
dword1(rv) = Tiny1; |
#else |
dword0(rv) = Tiny1; |
#endif /*_DOUBLE_IS_32BITS*/ |
/* The refinement below will clean |
* this approximation up. |
*/ |
} |
#endif |
} |
} |
|
/* 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,dval(rv), &bbe, &bbbits); /* rv = 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 Honor_FLT_ROUNDS |
if (rounding != 1) |
bs2++; |
#endif |
#ifdef Avoid_Underflow |
j = bbe - scale; |
i = j + bbbits - 1; /* logb(rv) */ |
if (i < Emin) /* denormal */ |
j += P - Emin; |
else |
j = P + 1 - bbbits; |
#else /*Avoid_Underflow*/ |
#ifdef Sudden_Underflow |
#ifdef IBM |
j = 1 + 4*P - 3 - bbbits + ((bbe + bbbits - 1) & 3); |
#else |
j = P + 1 - bbbits; |
#endif |
#else /*Sudden_Underflow*/ |
j = bbe; |
i = j + bbbits - 1; /* logb(rv) */ |
if (i < Emin) /* denormal */ |
j += P - Emin; |
else |
j = P + 1 - bbbits; |
#endif /*Sudden_Underflow*/ |
#endif /*Avoid_Underflow*/ |
bb2 += j; |
bd2 += j; |
#ifdef Avoid_Underflow |
bd2 += scale; |
#endif |
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); |
#ifdef Honor_FLT_ROUNDS |
if (rounding != 1) { |
if (i < 0) { |
/* Error is less than an ulp */ |
if (!delta->_x[0] && delta->_wds <= 1) { |
/* exact */ |
#ifdef SET_INEXACT |
inexact = 0; |
#endif |
break; |
} |
if (rounding) { |
if (dsign) { |
adj = 1.; |
goto apply_adj; |
} |
} |
else if (!dsign) { |
adj = -1.; |
if (!dword1(rv) |
&& !(dword0(rv) & Frac_mask)) { |
y = dword0(rv) & Exp_mask; |
#ifdef Avoid_Underflow |
if (!scale || y > 2*P*Exp_msk1) |
#else |
if (y) |
#endif |
{ |
delta = lshift(ptr, delta,Log2P); |
if (cmp(delta, bs) <= 0) |
adj = -0.5; |
} |
} |
apply_adj: |
#ifdef Avoid_Underflow |
if (scale && (y = dword0(rv) & Exp_mask) |
<= 2*P*Exp_msk1) |
dword0(adj) += (2*P+1)*Exp_msk1 - y; |
#else |
#ifdef Sudden_Underflow |
if ((dword0(rv) & Exp_mask) <= |
P*Exp_msk1) { |
dword0(rv) += P*Exp_msk1; |
dval(rv) += adj*ulp(dval(rv)); |
dword0(rv) -= P*Exp_msk1; |
} |
else |
#endif /*Sudden_Underflow*/ |
#endif /*Avoid_Underflow*/ |
dval(rv) += adj*ulp(dval(rv)); |
} |
break; |
} |
adj = ratio(delta, bs); |
if (adj < 1.) |
adj = 1.; |
if (adj <= 0x7ffffffe) { |
/* adj = rounding ? ceil(adj) : floor(adj); */ |
y = adj; |
if (y != adj) { |
if (!((rounding>>1) ^ dsign)) |
y++; |
adj = y; |
} |
} |
#ifdef Avoid_Underflow |
if (scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1) |
dword0(adj) += (2*P+1)*Exp_msk1 - y; |
#else |
#ifdef Sudden_Underflow |
if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) { |
dword0(rv) += P*Exp_msk1; |
adj *= ulp(dval(rv)); |
if (dsign) |
dval(rv) += adj; |
else |
dval(rv) -= adj; |
dword0(rv) -= P*Exp_msk1; |
goto cont; |
} |
#endif /*Sudden_Underflow*/ |
#endif /*Avoid_Underflow*/ |
adj *= ulp(dval(rv)); |
if (dsign) |
dval(rv) += adj; |
else |
dval(rv) -= adj; |
goto cont; |
} |
#endif /*Honor_FLT_ROUNDS*/ |
|
if (i < 0) { |
/* Error is less than half an ulp -- check for |
* special case of mantissa a power of two. |
*/ |
if (dsign || dword1(rv) || dword0(rv) & Bndry_mask |
#ifdef IEEE_Arith |
#ifdef Avoid_Underflow |
|| (dword0(rv) & Exp_mask) <= (2*P+1)*Exp_msk1 |
#else |
|| (dword0(rv) & Exp_mask) <= Exp_msk1 |
#endif |
#endif |
) { |
#ifdef SET_INEXACT |
if (!delta->x[0] && delta->wds <= 1) |
inexact = 0; |
#endif |
break; |
} |
if (!delta->_x[0] && delta->_wds <= 1) { |
/* exact result */ |
#ifdef SET_INEXACT |
inexact = 0; |
#endif |
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 ((dword0(rv) & Bndry_mask1) == Bndry_mask1 |
&& dword1(rv) == ( |
#ifdef Avoid_Underflow |
(scale && (y = dword0(rv) & Exp_mask) <= 2*P*Exp_msk1) |
? (0xffffffff & (0xffffffff << (2*P+1-(y>>Exp_shift)))) : |
#endif |
0xffffffff)) { |
/*boundary case -- increment exponent*/ |
dword0(rv) = (dword0(rv) & Exp_mask) |
+ Exp_msk1 |
#ifdef IBM |
| Exp_msk1 >> 4 |
#endif |
; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = 0; |
#endif /*!_DOUBLE_IS_32BITS*/ |
#ifdef Avoid_Underflow |
dsign = 0; |
#endif |
break; |
} |
} |
else if (!(dword0(rv) & Bndry_mask) && !dword1(rv)) { |
drop_down: |
/* boundary case -- decrement exponent */ |
#ifdef Sudden_Underflow /*{{*/ |
L = dword0(rv) & Exp_mask; |
#ifdef IBM |
if (L < Exp_msk1) |
#else |
#ifdef Avoid_Underflow |
if (L <= (scale ? (2*P+1)*Exp_msk1 : Exp_msk1)) |
#else |
if (L <= Exp_msk1) |
#endif /*Avoid_Underflow*/ |
#endif /*IBM*/ |
goto undfl; |
L -= Exp_msk1; |
#else /*Sudden_Underflow}{*/ |
#ifdef Avoid_Underflow |
if (scale) { |
L = dword0(rv) & Exp_mask; |
if (L <= (2*P+1)*Exp_msk1) { |
if (L > (P+2)*Exp_msk1) |
/* round even ==> */ |
/* accept rv */ |
break; |
/* rv = smallest denormal */ |
goto undfl; |
} |
} |
#endif /*Avoid_Underflow*/ |
L = (dword0(rv) & Exp_mask) - Exp_msk1; |
#endif /*Sudden_Underflow}*/ |
dword0(rv) = L | Bndry_mask1; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = 0xffffffff; |
#endif /*!_DOUBLE_IS_32BITS*/ |
#ifdef IBM |
goto cont; |
#else |
break; |
#endif |
} |
#ifndef ROUND_BIASED |
if (!(dword1(rv) & LSB)) |
break; |
#endif |
if (dsign) |
dval(rv) += ulp(dval(rv)); |
#ifndef ROUND_BIASED |
else { |
dval(rv) -= ulp(dval(rv)); |
#ifndef Sudden_Underflow |
if (!dval(rv)) |
goto undfl; |
#endif |
} |
#ifdef Avoid_Underflow |
dsign = 1 - dsign; |
#endif |
#endif |
break; |
} |
if ((aadj = ratio(delta, bs)) <= 2.) { |
if (dsign) |
aadj = dval(aadj1) = 1.; |
else if (dword1(rv) || dword0(rv) & Bndry_mask) { |
#ifndef Sudden_Underflow |
if (dword1(rv) == Tiny1 && !dword0(rv)) |
goto undfl; |
#endif |
aadj = 1.; |
dval(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; |
dval(aadj1) = -aadj; |
} |
} |
else { |
aadj *= 0.5; |
dval(aadj1) = dsign ? aadj : -aadj; |
#ifdef Check_FLT_ROUNDS |
switch(Rounding) { |
case 2: /* towards +infinity */ |
dval(aadj1) -= 0.5; |
break; |
case 0: /* towards 0 */ |
case 3: /* towards -infinity */ |
dval(aadj1) += 0.5; |
} |
#else |
if (Flt_Rounds == 0) |
dval(aadj1) += 0.5; |
#endif /*Check_FLT_ROUNDS*/ |
} |
y = dword0(rv) & Exp_mask; |
|
/* Check for overflow */ |
|
if (y == Exp_msk1*(DBL_MAX_EXP+Bias-1)) { |
dval(rv0) = dval(rv); |
dword0(rv) -= P*Exp_msk1; |
adj = dval(aadj1) * ulp(dval(rv)); |
dval(rv) += adj; |
if ((dword0(rv) & Exp_mask) >= |
Exp_msk1*(DBL_MAX_EXP+Bias-P)) { |
if (dword0(rv0) == Big0 && dword1(rv0) == Big1) |
goto ovfl; |
dword0(rv) = Big0; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv) = Big1; |
#endif /*!_DOUBLE_IS_32BITS*/ |
goto cont; |
} |
else |
dword0(rv) += P*Exp_msk1; |
} |
else { |
#ifdef Avoid_Underflow |
if (scale && y <= 2*P*Exp_msk1) { |
if (aadj <= 0x7fffffff) { |
if ((z = aadj) <= 0) |
z = 1; |
aadj = z; |
dval(aadj1) = dsign ? aadj : -aadj; |
} |
dword0(aadj1) += (2*P+1)*Exp_msk1 - y; |
} |
adj = dval(aadj1) * ulp(dval(rv)); |
dval(rv) += adj; |
#else |
#ifdef Sudden_Underflow |
if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) { |
dval(rv0) = dval(rv); |
dword0(rv) += P*Exp_msk1; |
adj = dval(aadj1) * ulp(dval(rv)); |
dval(rv) += adj; |
#ifdef IBM |
if ((dword0(rv) & Exp_mask) < P*Exp_msk1) |
#else |
if ((dword0(rv) & Exp_mask) <= P*Exp_msk1) |
#endif |
{ |
if (dword0(rv0) == Tiny0 |
&& dword1(rv0) == Tiny1) |
goto undfl; |
#ifndef _DOUBLE_IS_32BITS |
dword0(rv) = Tiny0; |
dword1(rv) = Tiny1; |
#else |
dword0(rv) = Tiny1; |
#endif /*_DOUBLE_IS_32BITS*/ |
goto cont; |
} |
else |
dword0(rv) -= P*Exp_msk1; |
} |
else { |
adj = dval(aadj1) * ulp(dval(rv)); |
dval(rv) += adj; |
} |
#else /*Sudden_Underflow*/ |
/* Compute adj so that the IEEE rounding rules will |
* correctly round rv + adj in some half-way cases. |
* If rv * ulp(rv) 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.) { |
dval(aadj1) = (double)(int)(aadj + 0.5); |
if (!dsign) |
dval(aadj1) = -dval(aadj1); |
} |
adj = dval(aadj1) * ulp(dval(rv)); |
dval(rv) += adj; |
#endif /*Sudden_Underflow*/ |
#endif /*Avoid_Underflow*/ |
} |
z = dword0(rv) & Exp_mask; |
#ifndef SET_INEXACT |
#ifdef Avoid_Underflow |
if (!scale) |
#endif |
if (y == z) { |
/* Can we stop now? */ |
L = (Long)aadj; |
aadj -= L; |
/* The tolerances below are conservative. */ |
if (dsign || dword1(rv) || dword0(rv) & Bndry_mask) { |
if (aadj < .4999999 || aadj > .5000001) |
break; |
} |
else if (aadj < .4999999/FLT_RADIX) |
break; |
} |
#endif |
cont: |
Bfree(ptr,bb); |
Bfree(ptr,bd); |
Bfree(ptr,bs); |
Bfree(ptr,delta); |
} |
#ifdef SET_INEXACT |
if (inexact) { |
if (!oldinexact) { |
dword0(rv0) = Exp_1 + (70 << Exp_shift); |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv0) = 0; |
#endif /*!_DOUBLE_IS_32BITS*/ |
dval(rv0) += 1.; |
} |
} |
else if (!oldinexact) |
clear_inexact(); |
#endif |
#ifdef Avoid_Underflow |
if (scale) { |
dword0(rv0) = Exp_1 - 2*P*Exp_msk1; |
#ifndef _DOUBLE_IS_32BITS |
dword1(rv0) = 0; |
#endif /*!_DOUBLE_IS_32BITS*/ |
dval(rv) *= dval(rv0); |
#ifndef NO_ERRNO |
/* try to avoid the bug of testing an 8087 register value */ |
if (dword0(rv) == 0 && dword1(rv) == 0) |
ptr->_errno = ERANGE; |
#endif |
} |
#endif /* Avoid_Underflow */ |
#ifdef SET_INEXACT |
if (inexact && !(dword0(rv) & Exp_mask)) { |
/* set underflow bit */ |
dval(rv0) = 1e-300; |
dval(rv0) *= dval(rv0); |
} |
#endif |
retfree: |
Bfree(ptr,bb); |
Bfree(ptr,bd); |
Bfree(ptr,bs); |
Bfree(ptr,bd0); |
Bfree(ptr,delta); |
ret: |
if (se) |
*se = (char *)s; |
return sign ? -dval(rv) : dval(rv); |
} |
|
#ifndef _REENT_ONLY |
|
double |
_DEFUN (strtod, (s00, se), |
_CONST char *s00 _AND char **se) |
{ |
return _strtod_r (_REENT, s00, se); |
} |
|
float |
_DEFUN (strtof, (s00, se), |
_CONST char *s00 _AND |
char **se) |
{ |
double retval = _strtod_r (_REENT, s00, se); |
if (isnan (retval)) |
return nanf (NULL); |
return (float)retval; |
} |
|
#endif |