/*
 * 5799-WZQ (C) COPYRIGHT IBM CORPORATION 1987
 * LICENSED MATERIALS - PROPERTY OF IBM
 * REFER TO COPYRIGHT INSTRUCTIONS FORM NUMBER G120-2083
 */
/* $Header:atof.c 12.0$ */
/* $ACIS:atof.c 12.0$ */
/* $Source: /ibm/acis/usr/src/lib/libc/ca/gen/RCS/atof.c,v $ */

#ifndef lint
static char *rcsid = "$Header:atof.c 12.0$";
#endif

#ifdef DEBUG
#define pp(x)	printf("x = %u\n",x)
#define msgs(x)	printf("x\n")
#else
#define pp(x)
#define msgs(x)
#endif DEBUG

#include <ctype.h>
#include <ieee.h>
#include <machine/fp.h>
/*****************************************************************

                ATOF -- ASCII TO DOUBLE CONVERSION
 
      Converts the string 'str' to ieee double form and returns 
      it. The result may be plus or minus zero, denormal, normal 
      or infinity. Maximum error in normalized result: .511 units
      in the last place (UNVERIFIED). 
 
      The syntax represents a rather generous interpretation of 
      "number": 

       {<sp>} [+|-] {<digit>} [.{<digit>}] [(e|E) [+|-] {<digit>}] 

      where <sp> is a whitespace character, <digit> is zero 
      through nine, [] is zero or one, and {} is zero or more. The 
      constant is delimited by any character that can't be part of 
      the constant, for example '\0' or other nondigit (including  
      +|-|.|e|E  when they are not syntactically valid).  A string 
      with no digits is possible; the mantissa is taken to be 0. 

*****************************************************************/
 

extern int _c19d ();
    /* args: arg1 = pointer to 2-word destination area
             arg2 = pointer to digit string (must be 19 or fewer)
             arg3 = digit count
       result is normalization amount (0 to -63, or -64 if dest=0)
    */ 

double atof(p)
char    *p; 
{ 
 char c; 
 int expon; /* Decimal exponent */
 int szc = 0; /* Count of leading zeros skipped */
 int ifg = 1; /* 0 once decimal point scanned */
 int idc = 0; /* Count of digits before decimal point */
 unsigned long twol[2]; /* 64-bit representation of fraction */
 char d[19]; /* Buffer for digits to be converted */
 char *dp = d;  /* Pointer to next position of d to fill */
 int dct;                /* Count of digits in d */
 long y = 1023 + 64; /* binary exponent */
/* 1023 is the bias for IEEE doubles; the 64 comes from the fact
   that we want to add the 64s complement of the shift amount q
   at a certain point; instead of doing that, we set y to an
   initial value 64 greater than the bias, then subtract q. We
   need the 64-s complement because we are shifting a right
   justified number in a 64-bit field until it is left justified.
   The greater the amount of left shift required, the less we have
   to add to y (since the greater the shift, the smaller the
   value.) */
 long cc; /* used in developing binary exponent */
 unsigned long t[3];      /* power of 10; product/quotient, etc. */
 unsigned short q; /* flag */
 unsigned short ts; /* shift amount */
 unsigned long sign = 0; /* BIT31 if '-' scanned */
 unsigned long gb; /* 1 if guard bit is 1 */
 unsigned long ub; /* 1 if ulp bit is 1 */
 unsigned long sb = 0; /* 1 if any bits beyond guard bit are 1 */
 unsigned long ud; /* 1 if UP and sign 0 or DOWN and sign 1 */
 unsigned long nr; /* 1 if TONEAREST */
 unsigned long r; /* round amount */
 unsigned long mode; /* rounding mode desired */
 union {
        double dres;
        struct {
                unsigned long dresh, dresl;
               } u2l;
       } ures;
 
 while (isspace(c = *p++)); 	/* skip leading blanks */ 
 
 switch (c) /* scan optional sign */ 
     { 
      case '-': sign = BIT31; /* NOTE: fall through! */
      case '+': c = *p++; 
     } 

 switch (c)
     {
      case 'I': case 'i': /* INF */
          if ( ((*p == 'N') || (*p == 'n')) &&
	       ((*(p+1) == 'F') || (*(p+1) == 'f')) )
               {
                ures.u2l.dresh = DEXPBITS | sign;
                ures.u2l.dresl = 0x0;
                return(ures.dres);
               }
          break;
      case 'N': case 'n':  /* NaN() */
          if ( ((*p == 'A') || (*p == 'a')) &&
               ((*(p+1) == 'N') || (*(p+1) == 'n')) &&
               (*(p+2) == '(') && (*(p+3) == ')') )
               {
                ures.u2l.dresh = DEXPBITS | DoubleStoQ | sign;
                ures.u2l.dresl = 0x0;
                return(ures.dres);
               }
     }

 for (;;) /* start collecting the number */ 
     { 
      if (c == '.' && ifg) 
          ifg = 0; /* accept one decimal point */
      else if (c < '0' || c > '9')
                break; 
           else
                {
                 idc += ifg; /* bump count of integer digits */
                 if (c == '0' && dp == d)
                    szc++; /* ignore leading zeros */
                 else if (dp < d+19)
                         *dp++ = c;
                           /* store only the first 19 digits; */
                           /* that many fit as a 2-word integer */
                }
      c = *p++;
     } 

 dct = dp-d;
/* The function _c19d(), given a pointer to a character buffer d,
    and a count of digits dct, picks up ASCII digits and converts
    them to a right-justified integer, then shifts this integer
    left until it is left adjusted (normal); the result of the
    function is the negative of the amount of shift, 0 to -63, or
    -64 if the number is zero. */  
 y += _c19d(twol,d,dct);
 msgs(***** after conversion by c19d *****);
 pp(twol[0]);
 pp(twol[1]);
 if (y == 1023)
     {/* further scaling won't help much, */
      ures.u2l.dresh = sign;	/* but give zero the proper sign    */
      ures.u2l.dresl = 0;
      return(ures.dres);
     }

 expon = idc-szc-dct;
 if ((c == 'e') || (c == 'E'))
     { /* test for exponent */ 
      if (!isspace(*p)) /* disallow  E whitespace number     */
          expon += atoi(p); 
     }
 
 if (expon < -345)
     expon = -345;
 if (expon > 340)
     expon = 340;

      /* The normal number in twol represents a number having its least 
      significant bit scaled by 1e0, or 1. If expon isn't zero, it 
      is necessary to scale the number by the power of ten 
      appropriate to the argument expon. Scaling is by 
      multiplication by 10**expon if expon is positive, or by
      dividing by 10**-expon if expon is negative. Both the 
      multiplication and division processes used yield three-
      longword results. Division by a positive power of ten, 
      rather than multiplication by a negative power of ten, is 
      used in scaling numbers less than one in magnitude. In order 
      to satisfy the requirements of IEEE standard 754 it has been 
      found that multiplying by negative powers of 10, even when carried 
      to 96 bits of precision, is inadequate to correctly convert 
      all numbers of the form M times 1eN, where N is between -28 
      and -1. See section 5.6 of IEEE 754. As examples of numbers 
      which are particularly difficult to convert, consider: 

              1.1273711918250669e-11  
       or    1.8ca891d088cce 7fff ffff ffff fffa in hexadecimal

        and   4.8002944080187915e-11
        or    1.a63d0559dac31 8000 0000 0000 1b in hexadecimal. */

      if (expon)
          {
           _TenToThe (((expon > 0) ? expon : -expon), t); /* t is result. */
	   msgs(***** power of 10 *****);
	   pp(t[0]);
	   pp(t[1]);
	   pp(t[2]);
           if (expon > 0) /* avoid scaling by 10**0 */
               { /* Scale by multiplying. */
                _mulu32 (t, twol); /* t holds product. */
		msgs(***** after _mulu32 *****);
		pp(t[0]);
		pp(t[1]);
		pp(t[2]);

                 /* The result in t is three words, with
                    the first word either normal or with
                    at most one leading zero bit. We must
                    normalize t if necessary, and subtract
                    the normalizing amount from y. */

                q = _shlun (t);
		msgs(***** after _shlun *****);
		pp(t[0]);
		pp(t[1]);
		pp(t[2]);
                y -= q;
                sb |= t[2] << q; /* Remember ignored bits */
               } /* t is now normal. */
           else /* expon is negative */
               { /* Scale by dividing. */
/* Normalized powers of 10 from 1 through 6 can fit in a half
   word. Dividing by these powers of ten is made efficient by
   using the function _dknu16(). For powers of ten greater than
   6 it is necessary to use the more general division routine
   _dknu(). */
                if (expon > -7)
                    _dknu16 (twol, t); /* t is 1 halfword. */
                else
                    _dknu (twol, t); /* t is 2 to 5 halfwords */

                /*  The result in t is three words, with
                    the first word in one of the forms
                    0000gxxx or 0001xxxx, where g is one
                    of 89abcdef. We must normalize this
                    by shifting left either 15 or 16, and
                    subtracting 1 from y if the form is
                    0001xxxx. Only the first two words of
                    t need to be shifted, since only 54
                    bits are needed to give the properly
                    rounded fraction. */
		    msgs(***** after _dknu *****);
		    pp(t[0]);
		    pp(t[1]);
		    pp(t[2]);

                q = t[0] > 65535;
                ts = 16 - q; /* Shift amount, 15 or 16. */
                t[0] = t[0] << ts | t[1] >> 32 - ts;
                t[1] = t[1] << ts | t[2] >> 32 - ts;
		msgs(***** after 15/16 shift *****);
		pp(t[0]);
		pp(t[1]);
		pp(t[2]);
                y -= !q;
                sb |= t[2] << ts; /* remember lost bits */
               }
          }
      else  /* expon is 0, no scaling needed. */
          {
           y--;           /* But we have to abate y */
/* Why y-- ? Well, we fell through to this spot because the
   expon value is zero, and we wanted to avoid the extra
   (and useless) work of scaling by 10**0, which is 1. 
   However, if we HAD scaled the number, using either times
   or divide, we would have had to shift the argument to
   normalize it, whatever value the argument had, from the
   smallest number to the largest possible number. In order
   to account for this normalization which we avoided doing
   yet must account for, we subtract 1 from y. */
           t[0] = twol[0]; /* and move number */
           t[1] = twol[1]; /* to t. */
          }

      /* The binary exponent y has to be adjusted by the base-2 
      equivalent of the ten's exponent expon. This is obtained by 
      multiplying the tens' exponent expon by an approximation to the 
      base-2 log of 10, and taking the integer part of the 
      product. A sufficiently accurate approximation to the base-2 
      log of 10 is given by 3.3219, and multiplication by this 
      value may be achieved largely by shifting and adding rather 
      than by actual multiplication. The value 3.3219 is given by 
      108853/32768. The result, which may be positive, negative, 
      or zero, is added to y. */ 
           
      /* Consider multiplication by 108853: this is the same as
      multiplication by 53 plus multiplication by 256 plus
      multiplication by 108544; and 108854 is 53 * 2048. We
      trust the compiler does an efficient job in generating
      the multiplies by the smaller numbers, and then we combine
      them to get the desired end product. */

      cc = expon * 53; /*Good C compilers shift & add for this. */
      y += (cc + (expon * 256) + (cc * 2048)) >> 15;

      /* A normal number will have a binary exponent y between 1 
      and 2046, inclusive. If y is less than 1, the result will 
      probably be either subnormal or zero. */ 
 
      if (y < 1)
          { /* Probably subnormal or zero result. */
           if (y < -31)
               {t[2] = t[1]; /* Save sticky bits */
                t[1] = t[0]; /* Shift right 32. */
                t[0] = 0;
                if (y < -52) 
/* If y is less than -52 we have a number too small to represent,
   so we substitute zero for the unrepresentable too-small value. */

                    t[1] = y = 0; /* Force zero result. */
                else y += 32;
               }

                    /* We shift t right by 1-y rather
                       than by -y because a subnormal
                       number has a bias of only 1022,
                       not 1023. */

           t[2] = t[1] << (32 - (1 - y)) | t[2] >> (1 - y);
           t[1] = t[0] << (32 - (1 - y)) | t[1] >> (1 - y);
           t[0] >>= (1 - y);
           y = 0; /* Ensure zero exponent */
          }

      /* If y is greater than 2046 the result is infinite. */

      if (y >= MAX_DEXPON)
          {
           t[0] = t[1] = 0;
           y = MAX_DEXPON;
          }

      /* Two things need to be done to the fraction to put it in 
      final form. First, the leading 1-bit must be deleted, since 
      for all but subnormal numbers this is implied. (We have 
      already accounted for subnormal numbers properly by giving a 
      1-bit extra shift, to put the leading bit out of harm's 
      way). */ 

      t[0] &= LO31BITS;

      /* The fraction is now shifted right so that its first 52 
      bits are right-adjusted in the result. */ 

      twol[0] = t[0] >> 11;
      twol[1] = t[0] << 21 | t[1] >> 11;
      msgs(***** after right asjusting fraction *****);
      pp(twol[0]);
      pp(twol[1]);

      /* The rounding amount r is determined and added to twol. 
      This is 1 if and only if (a) the rounding mode is TONEAREST 
      and there is a guard digit and either the ulp bit or the 
      sticky bit is 1; or (b) the rounding mode is UPWARDS with 
      sign 0 or DOWNWARDS with sign 1, and either the ulp or the 
      sticky bit is 1. */ 
      
      /* For those unfamiliar with the jargon in the above 
      comment: the ulp bit is the last bit in the result fraction; 
      the guard bit is the first truncated bit, and will be 1 if 
      the truncated part is a half or more. The sticky bit is a 
      fictitious summary bit and will be 1 if any of the bits 
      following the guard bit are 1. There are four rounding modes 
      provided: TONEAREST rounds to the nearest value: up if the 
      truncated portion is greater than a half, down if it's less 
      than a half; if it's exactly a half the result fraction is 
      rounded up if it's odd (thereby making it even) and down if 
      it's already even. This is called unbiased rounding, since 
      in principle half of such numbers will be rounded up and half 
      will be rounded down. Unbiased rounding prevents the drift 
      phenomenon which occurs with biased rounding, where 1 is 
      added to the fraction whenever the truncated portion is 
      exactly a half. See D.E. Knuth, Seminumerical Algorithms, 
      section 4.2.2.A, before and after Theorem D. UPWARD rounding 
      gives the value closest to and no less than the untruncated 
      value; that is, it rounds up whenever there are 1's in the 
      truncated portion and the sign is positive. DOWNWARD 
      rounding gives the value closest to and no more than the 
      untruncated value; that is, it rounds down if the sign is 
      positive, effectively just truncating the fraction. When the 
      sign is negative, the rules for UPWARD and DOWNWARD are 
      reversed. In TOZERO rounding the fraction is always 
      truncated. */
      
      swapround (mode = swapround (0)); /* get round mode */
      gb = t[1] & BIT10; /* Isolate the guard bit */
      ub = t[1] & BIT11; /* Isolate ulp bit */
      sb |= (t[1] & LO10BITS); /* Isolate sticky bits */
      ud = ((mode == UPWARD) && !sign) || ((mode == DOWNWARD) && 
           sign); /* 1 if UP and positive or DOWN and negative */
      nr = mode == TONEAREST; /* 1 if TONEAREST */
      r = (nr && gb && (ub || sb)) || (ud && (gb || sb));

      /* Round with carry if needed; this may propagate a 1 into 
      the exponent field. This is fine since when we add the 
      exponent field in it adjusts the exponent properly, perhaps 
      even converting a subnormal number into a normal number, or 
      a normal number into an infinity. */ 

      twol[1] += r;
      twol[0] += twol[1] < r; 
      msgs(***** after rounding *****);
      pp(twol[0]);
      pp(twol[1]);
      
      /* The binary exponent is shifted left 20 to occupy bit 
      positions 1-11 of a longword, and is added to the result 
      first longword. It is added rather than or-ed because we 
      have to add 1 to y if the rounding operation produced an 
      all-zero fraction (now we see the utility of possibly having 
      propagated a 1 into the low-order bit of the exponent field 
      in the rounding step). */ 

      twol[0] += y << 20;
      msgs(***** after exponent field inserted *****);
      pp(twol[0]);

 /* The last step ors the sign bit into the first bit position of 
 the first word of the result. We know this bit is currently a 
 zero because either the result is all zero, or we shifted the 
 words holding the unsigned long fraction values right just a 
 moment ago. */ 

 twol[0] |= sign;

 /* The result is complete. We return it as a floating-point 
 double. */

 ures.u2l.dresh = twol[0];
 ures.u2l.dresl = twol[1];
 return (ures.dres);
}
