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

Subversion Repositories or1k

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

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 y, z;
  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)
        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)
        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 (s00, se);
}
 
#endif

Browse

Tools

Subversion Repositories or1k

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

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