URL https://opencores.org/ocsvn/or1k/or1k/trunk

Subversion Repositories or1k

[/] [or1k/] [trunk/] [newlib-1.10.0/] [newlib/] [libc/] [stdlib/] [strtod.c] - Blame information for rev 1773

Go to most recent revision | Details | Compare with Previous | View Log


/*
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 <_ansi.h>
#include <reent.h>
#include <string.h>
#include "mprec.h"
 
double
_DEFUN (_strtod_r, (ptr, s00, se),
        struct _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;
  long e;
  _CONST char *s, *s0, *s1;
  double aadj, aadj1, adj;
  long L;
  unsigned long z;
  __ULong y;
  union double_union rv, rv0;
 
  _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')
    {
      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 < 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)
            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')
            {
              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
            e = 0;
        }
      else
        s = s00;
    }
  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) != 0)
        rv.d *= tens[i];
      if (e1 &= ~15)
        {
          if (e1 > DBL_MAX_10_EXP)
            {
            ovfl:
              ptr->_errno = ERANGE;
#ifdef _HAVE_STDC
              rv.d = HUGE_VAL;
#else
              /* Can't trust HUGE_VAL */
#ifdef IEEE_Arith
              word0 (rv) = Exp_mask;
#ifndef _DOUBLE_IS_32BITS
              word1 (rv) = 0;
#endif
#else
              word0 (rv) = Big0;
#ifndef _DOUBLE_IS_32BITS
              word1 (rv) = Big1;
#endif
#endif
#endif
              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) != 0)
        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;
  return sign ? -rv.d : rv.d;
}
 
#ifndef NO_REENT
 
double
_DEFUN (strtod, (s00, se),
        _CONST char *s00 _AND char **se)
{
  return _strtod_r (_REENT, s00, se);
}
 
float
_DEFUN (strtodf, (s00, se),
        _CONST char *s00 _AND
        char **se)
{
  return _strtod_r (_REENT, s00, se);
}
 
#endif

Browse

Tools

Subversion Repositories or1k

[/] [or1k/] [trunk/] [newlib-1.10.0/] [newlib/] [libc/] [stdlib/] [strtod.c] - Blame information for rev 1773

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