Macintosh-Tools/Retro68
1
0
Fork 0
mirror of https://github.com/autc04/Retro68.git synced 2024-12-13 03:29:50 +00:00
Code Issues Projects Releases Wiki Activity
6e0a59fa92
Retro68/gcc/libjava/classpath/native/fdlibm/strtod.c
Wolfgang Thaller aaf905ce07 add gcc 4.70
2012-03-28 01:13:14 +02:00

720 lines
15 KiB
C
Raw Blame History

/*
FUNCTION
<<strtod>>, <<strtodf>>---string to double or float
INDEX
strtod
INDEX
_strtod_r
INDEX
strtodf
ANSI_SYNOPSIS
#include <stdlib.h>
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 <stdlib.h>
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 <<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 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]>.
<<strtodf>> 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.
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) 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 <string.h>
#include <float.h>
#include <errno.h>
#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 = NULL, *bb1, *bd = NULL, *bd0, *bs = NULL, *delta = NULL;
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;
}
Reference in New Issue View Git Blame Copy Permalink