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

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

#define DEBUG 0
/****************************************************************


   DECIMAL FROM IEEE DOUBLE FORMAT BINARY FLOATING-POINT

The function ftoa converts an IEEE 754 double format floating-
point number to decimal. The conversion produces a maximum of 
seventeen significant digits.  When rounding to nearest, the error 
in the converted result is no more than 0.002 units in the 
destination's least significant digit.

This function will round either to a specified number of decimal 
digits (FORTRAN E format), or to a specified number of decimal 
places (FORTRAN F format). The calling parameter eflag 
distinguishes between these two.

Two functions, ecvt and fcvt, are provided, which serve to call
ftoa with E format and F format specified, respectively.

All of the rounding modes specified by IEEE 754 are provided: 
TONEAREST, UPWARD, DOWNWARD, and TOWARDZERO. In addition, when the 
precision warrants it and CVT_ROUND is 0, a technique 
competitive with TONEAREST is provided, called TOROUNDEST. The 
difference between these two modes can be shown by example (assume 
that maximum precision of seventeen digits is asked for): 

If the number 0.07 is input, the closest-fitting floating-point 
number to it is represented by the hexadecimal form 3fb1 eb85 
1eb8 51ec. This number, being finite, has an exact decimal 
representation, which is:

0.070000 000000 000006 661338 147750 939242 541790 008544 921875

We can compare how TONEAREST and TOROUNDEST will treat the three 
consecutive floating-point numbers of which this is the center 
one. We'll call them OH7m, OH7, and OH7p.  

At some point in the conversion process a three-digit number in 
radix 1e9 is developed. Usually, this will represent a twenty-
digit decimal number; sometimes it will represent a twenty-one 
digit decimal number. From this radix 1e9 number it is easy to 
convert to a decimal number. In the table below a bar (|) is 
inserted in the third radix 1e9 digit, just after the seventeenth 
decimal digit. In the column headed "radix 1e9 form" are shown the 
values as computed by this function. In the column headed 
"accurate" are shown the actual next six digits following the bar.  

   number         radix 1e9 form             accurate     
    OH7m        69 999999999 999992|782       |783550...
    OH7         70 000000000 000006|660       |661338...
    OH7p        70 000000000 000020|538       |539125...

In TONEAREST mode, rounding these numbers requires adding 1 to the 
last digit kept if the dropped digits, considered as a fraction, 
are more than a half, or are equal to a half and the last digit 
kept is odd. Following this rule, we get the values: 

   number      rounded to nearest 17   number of digits
    OH7m        69 999999999 999993          17
    OH7         70 000000000 000007          17
    OH7p        70 000000000 000021          17

Notice that the above table implies that not only the number .07, 
but also all of the seventeen-digit numbers beginning with .07 
and ending with 1 (.070000000000000001), 2, 3, 4, 5, 6, and 7, 
when entered will result in the same internal number being 
generated; furthermore, an equal range of numbers on the other 
side will also be represented by the same internal number (those 
ending with 8, 9, and so forth). We'll call this internal number 
the "internal representative" of the range of external numbers.  
Conversely, we'll call the range of external numbers which map 
into the same internal representative "the set of external 
representatives".  

In the TOROUNDEST mode, on the other hand, rounding requires 
consideration of two things: what the interval is between two 
consecutive twenty-digit numbers, and what seventeen-digit numbers 
are in the set of external representatives that lies half this 
amount on either side of the number. For example, the interval 
between the twenty-decimal versions of OH7m and OH7, which is the 
same as the interval between OH7 and OH7p, is 13|878 units. Half 
of this is 6|939 units. In order to keep two adjacent numbers from 
claiming the same round number halfway between them, we abate this 
half-interval by a small amount, actually extending an interval of 
only 6|932 units on either side of the number. We then look for 
the roundest number of seventeen or fewer significant digits in 
the set of external representatives which lies in this interval. 
The roundest number is the one with the largest number of trailing 
zeros in it. The range of possibilities for the three numbers 
above is as follows: 

       number    low end    high end   roundest value
                 of range   of range      in range
        OH7m      85|850     99|714        90|000
        OH7       99|728     13|592        00|000
        OH7p      13|606     27|470        20|000

Applying this procedure, we get the following results:   

   number     rounded to roundest 17   number of digits
    OH7m        69 999999999 99999           16
    OH7         7                             1
    OH7p        70 000000000 00002           16

The IEEE 754 standard requires that conversion from internal form 
to external form and back again shall be the identity as long as 
the decimal string is carried to the maximum precision of 
seventeen digits. The TOROUNDEST mode chooses to represent the 
internal form of a floating point number by that one of the set of 
external representatives which is easiest to look at, and it will 
be true that if that roundest external representative is again 
entered, it will result in the same internal representative being 
developed. It is of course also true that the external 
representative found by the TONEAREST mode, if entered, will 
result in the development of the same internal representative. The 
difference, as stated above, is that the external representative 
found by the round to roundest mode may be easier to look at. */ 

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

A list of the powers of ten from 1e1 through 1e8: */

static unsigned long p[] = {
                            10, 100, 1000, 10000, 100000, 
                            1000000, 10000000, 100000000
                           };

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

The definitions of the IEEE rounding modes TONEAREST, etc., and
the swapround() function are given in ieee.h. Definitions for
various machine word bit configurations are given in
machine/fp.h. */ 

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

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

The static variable _ieeerounding controls whether IEEE-style
rounding will be made (if its value is CVT_IEEE) or whether
rounding to roundest will be attempted (if its value is
CVT_ROUND). Its default is CVT_ROUND, and it can be changed via
cvtrounding(), which changes _ieeerounding and returns its
previous value. */

static int _ieeerounding = CVT_ROUND;

int cvtrounding(class)
int class;
{int tmp;


	tmp = _ieeerounding;
	_ieeerounding = class;
	return (tmp);
}
 
/*****************************************************************

A definition for a minimum macroinstruction will be useful. */
#ifdef __HIGHC__
#define min(a,b) _min(a,b)
#else
#define min(a,b) ((a)<(b)?(a):(b))
#endif

/*****************************************************************
 
     The functions ecvt and fcvt are provided simply to make it 
     convenient to use ftoa. Their sole purpose is to provide for 
     setting the eflag parameter in the call to ftoa. */ 

char* ftoa();

char* ecvt(arg, ndigits, decpt, sign)
double arg;
int ndigits, *decpt, *sign;
{
 return ftoa(arg, ndigits, decpt, sign, 1);
}

char* fcvt(arg, ndigits, decpt, sign)
double arg;
int ndigits, *decpt, *sign;
{
 return ftoa(arg, ndigits, decpt, sign, 0);
}

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

The function ftoa gives as result a pointer to a character vector 
of ASCII digits representing the properly rounded decimal value 
corresponding to the argument, depending on the setting of the 
rounding parameter, _ieeerounding, the number of digits requested, 
and the eflag setting.
 
If the result is not zero, the leading  digit will be nonzero.

If the result is zero, a single character "0" is pointed to. 

When rounding TONEAREST or TOROUNDEST the function is accurate to 
about 0.002 units in the seventeenth digit. As side effects the 
function gives pointers to the location of the result's decimal 
place and  its sign. */ 

/************ FUNCTION SYNTAX AND ARGUMENTS *******************/

char *ftoa (arg, ndigits, decpt, sign, eflag)

/**************** ARGUMENT DECLARATIONS ***********************/

double arg;        /* Argument: A floating point  number */
int   *decpt;      /* Side effect: Location of decimal point in 
                      result. Its value gives the digit in the 
                      result to the right of which the decimal 
                      point occurs. A zero value means the decimal 
                      point occurs immediately to the left of 
                      the most significant digit; a negative value
                      means that the decimal point occurs the 
                      number of places to the left of the most 
                      significant digit given by the magnitude of 
                      decpt. */
int    eflag;      /* Parameter: 0 = F format Fortran style 
                      conversion; 1 = E format. */
int    ndigits;    /* Parameter: Required precision (if eflag is 
                      0) or number of fraction digits (if eflag is 
                      1). Meaningful values are 0 or greater. */
int   *sign;       /* Side effect: Result's sign: 1 if negative,
                      0 otherwise. */

/************ LOCAL VARIABLE DECLARATIONS *********************/

/* This union permits converting the sense of a value from double 
floating point to a pair of unsigned long integers. */

{
 union
     {
      double           d;
      unsigned long    i[2];
     }                       a;

 /* aldi is used to record the actual leading digit index in tnf; 
 for subnormal numbers this may be different from the value in 
 ldi. */

 int                         aldi;

 /* borrow is a logical value recording the borrow status of rem 
 with respect to tnh, and is 1 if rem is less than tnh. */ 

 int                         borrow;

 /* carry is a logical value recording the carry status of rem + 
 tnh with respect to tnt, and is 1 if this sum exceeds tnt. */ 

 int                         carry;

 /* dtri gives the index in tnf of the digit which is to be 
 rounded in the case of ndigits rounding. */

 int                         dtri;

 /* dtrir gives the index in tnf of the digit to be rounded
 in the case of rounding to roundest. */

 int                         dtrir;

 /* exp is the unbiassed exponent of the argument. */

 int                         exp;

 /* fraction is the fractional part of the argument. */

 unsigned long               fraction[2];
 
 /* guard is four times the value of the digit following the least 
 significant digit kept. */ 

 unsigned long               guard;

 /* i is an induction variable used to control the padding out of 
 the result with trailing character "0"s and in the computation 
 of sigdigs. */

 int                         i;

 /* j is an induction variable used to locate the leading non '0' 
 in the result tnf. */

 int                         j;

 /* ldi is the index in tnf of the leading nonzero character 
 digit. */

 int                         ldi;

 /* lzcsn is the number of leading zeros in the fraction of the 
 argument. This is 11 for normal numbers, but varies from 12 to 63 
 for subnormal numbers. */ 

 int                         lzcsn;
  
 /* maketnh is a copy of the first element of the scaling power of 
 ten ts; it is used to make the half-ulp value used in "round to 
 roundest" mode. */ 

 unsigned long               maketnh;

 /* n is the unsigned long two-element form of the argument. */

 unsigned long               n[2];

 /* normal is the exponent part of the argument, and is nonzero if 
 the argument is normal, and zero otherwise. */ 

 unsigned long               normal;

 /* np is the factor, derived from the argument exponent, used in 
 determining the shift amount needed to form a 20- or 21-decimal 
 integer. */ 

 long                        np;

 /* nr is 1 if the rounding mode is TONEAREST, and is 0 otherwise. 
 */

 unsigned long               nr;

 /* ns is a two-longword form of the argument significand. */ 

 unsigned long               ns[2];

 /* oops21 is a logical value having the value one if we have 
 produced a 21-digit result rather than a 20-digit result. */ 

 int                         oops21;

 /* parity is 0 if the least significant digit kept is even before 
 rounding, and is 1 if it is odd. */

 unsigned long               parity;

 /* poweroftwo is a logical value, and is 1 if the argument is an 
 exact power of two. */ 

 int                         poweroftwo;

 /* ptn determines the power of ten used in scaling the significand 
 of the argument so that it can be formed into a 20-decimal (or 
 occasionally a 21-decimal) integer. */ 

 int                         ptn;

 /* r is initially the rounding amount, and subsequently the carry 
 amount. */

 unsigned long               r;

 /* rem is the remainder upon division of the last element of tni 
 with tnt, and is used in determining how to round to roundest. */ 

 unsigned long               rem;

 /* rm is the rounding mode indicator, and should be one of 
 TONEAREST, UPWARD, DOWNWARD and TOWARDZERO. */

 unsigned long               rm;

 /* sticky is a summary value of the digits following the guard 
 digit, and is 2 if any of these digits are nonzero, and is 0 
 otherwise. */

 unsigned long               sticky;

 /* sum records the result of evaluating the polynomial:

      4*guard + 2*sticky + parity

 and is used in determining the amount used in rounding. */
 
 unsigned long               sum;  
 
 /* tnf is the ASCII character form of the result significand. */

 static char                 tnf[27];

 /* tnh is a little less than half the difference between any two 
 consecutive numbers having the same exponent as the argument, and 
 is central to the "round to roundest" process. */ 

 unsigned long               tnh;

 /* tni is the three-element integer, representing a 20- or 21- 
 decimal number, formed from ts by performing a triple-precision 
 shift by the amount tnishift. */ 

 unsigned long               tni[3];

 /* tnishift is the amount by which the scaled number must be 
 shifted right in order to shift out the entire fractional part, 
 thus leaving only the integer part. */ 

 int                         tnishift;

 /* tnp is the sum of the shift amounts np and tp. */

 long                        tnp;

 /* tnt is the next power of 10 greater than tnh, and is used in 
 producing quo and rem from the last element of tni. */ 

 int                         tnt;

 /* tp is the factor, derived from ptn, used in determining the 
 amount of shift needed to form the 20 (or 21) decimal integer 
 part of a number. */ 

 int                         tp;

 /* ts is initially the normalized three-longword representation 
 of the power of ten used in scaling the significand (exact for 
 powers 1 through 27, and accurate to amounts ranging from 69 to 
 72 bits for other powers); subsequently it holds the scaled 
 product. */ 
 
 unsigned long               ts[3];
 
 /* tx records the ten's exponent of the argument; in rare cases 
 it may be high by one, because of the approximation used in 
 forming it. */ 

 int                         tx;

 /* ud is 1 if the rounding mode is either UPWARD with sign 0, or 
 DOWNWARD with sign 1. */

 unsigned long               ud;

 /* whichp is the index in p of the smallest integral power of ten 
 equal to or greater than tnh. */

 int                         whichp;

 /* x is a temporary used in computing tp. */

 int                         x;
 
/*********************** ALGORITHM ****************************
     
     Convert the argument from floating point type to unsigned 
     long integer type. */

 a.d = arg;
 n[0] = a.i[0]; 
 n[1] = a.i[1];

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

     Decompose the argument into its constituent parts: unbiassed 
     exponent, sign, and fraction. */

 exp = ((n[0] & DEXPBITS) >> 20) -1023;
 *sign = n[0] >> 31;          /*** Side effect of function ***/
 fraction[0] = n[0] & LO20BITS;
 fraction[1] = n[1];

/***************************************************************
     
     Note whether the fraction is all zero (indicating that the 
     argument is an exact power of two). */

 poweroftwo = !(fraction[0] || fraction[1]);

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

     Dispatch the cases of zero, infinity, and "not a number". */

 *decpt = 0;               /* Default for 0, INF, and NAN */
 
 if ((exp == -1023) && poweroftwo)    
     return("0");                  /* Zero */

 if (exp == 1024)
     {if (poweroftwo)
          return("INF");           /* Infinity */
      else
          return("NAN()");         /* Not a Number */
     }

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

     Set dtrir and r to useful initial values (dtrir to 25 so
     that dtri will never exceed it, and r to 0 so that if
     rounding amount never is 1 we don't attempt rounding. */

 dtrir = 25;
 r = 0;

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

     Find tx, the integer exponent of the power of ten closest to 
     the argument. The approximation used will occasionally err by 
     giving a value low by one. */

 tx = _l10(n);

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

     Find ptn, the exponent of the power of ten to be used in 
     scaling the argument. This exponent is computed by 19 - tx, 
     and is chosen so that the integer part of the product of the 
     argument and ten to the ptn power is a twenty-decimal number.  
     For example, if tx is 100, ptn will be -81; if tx is -100, ptn 
     will be 119. Because of the error possible in tx, the product 
     will occasionally have twenty-one digits. This produces no 
     great complications. */

 ptn = 19 - tx;
   
/****************************************************************

     Find ts, three unsigned long integers giving an approximation 
     to the significand of ten to the ptn power. We may need the 
     first element of ts in order to compute the half-ulp value 
     later, so save it in maketnh. */

 _TenToThe ( (long) ptn, ts);
 maketnh = ts[0];
 
/****************************************************************

     Find ns, the left-adjusted significand of the argument. For 
     normal numbers, the significand is the fraction with a 
     leading 1-bit added (since in the IEEE 754 representation the 
     leading 1 bit of normal numbers is implied); for subnormal 
     numbers -- numbers with biassed exponent 0, and unbiassed 
     exponent -1023 -- the significand is the same as the 
     fraction. The _shlun function left-adjusts its argument and 
     returns the amount of shift, and this is saved in lzcsn. */ 

 normal = exp != -1023;
 ns[0] = fraction[0] | normal << 20;
 ns[1] = fraction[1];
 lzcsn = _shlun(ns);

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

     Find ts, the first three longwords of the product of ts and 
     ns. */

 _mulu32 (ts, ns);
     
/****************************************************************

     Form tni, the integer part of this product, by shifting ts 
     right until the fractional part is lost. The amount of shift 
     required is a function of three quantities: 
     
        96 - (tp + np)
   
     The 96 arises from the fact that ts is three longwords, or 96 
     bits long. tp is the integral number of bits in the power of 
     ten ptn, and is given by 1 plus the floor of 3.3219 times 
     ptn. The value 3.3219 approximates the ratio (ln 10)/(ln 2); 
     adding 1 to to the floor of this gives the integral number of 
     bits in ptn. 
     
     Multiplication by 3.3219 is equivalent to multiplication by 
     108,853 followed by a right shift of 15, since 3.3219 is 
     close to 108853 / 32768. Furthermore, multiplication by 
     108,853 can be achieved by shifts and adds in a relatively 
     straightforward way, since the binary representation of 
     108,853 is 
     
              1 1 0 1 0 1 0 0 1 0 0 1 1 0 1 0 1
              \_________/           \_________/
     
     and the indicated groups match; this allows us to form ptn 
     times 53 by three shifts and adds (53 is 1 1 0 1 0 1 in 
     binary) and then ptn times 108544 by a single shift (108544 in 
     binary is 1 1 0 1 0 1 0 0 0 0 0 0 0 0 0 0 0), and then ptn 
     times 256 by a single shift (256 is 1 0 0 0 0 0 0 0 0 in 
     binary); adding these three gives ptn times 108853, since 53 + 
     256 + 108544 gives 108853.

     For normal arguments, np is the unbiassed exponent of the 
     argument; for subnormal arguments it is this exponent minus 
     one less than the number of leading zeros in the significand 
     (one less because the bias of subnormal numbers is only 1022, 
     not 1023). For normal numbers the high order digit of tni 
     will vary between 0 and 5. The expression 12 + !normal - 
     lzcsn evaluates to 0 for normal numbers. */ 
   
 x = ptn * 53;
 tp = 1 + ((x + (x << 11) + (ptn << 8)) >> 15); /* 1 + floor of 
                                             ptn times 3.3219 */
 np = exp + 12 + !normal - lzcsn;
 tnp = np + tp;
 tnishift = 96 - tnp;
 tni[2] = ts[2] >> tnishift | (ts[1] << (32 - tnishift));
 tni[1] = ts[1] >> tnishift | (ts[0] << (32 - tnishift));
 tni[0] = ts[0] >> tnishift;
     
/****************************************************************

     Record in sticky whether any of the bits shifted out of ts 
     are nonzero. */

 sticky = 0 != ts[2] << (32 - tnishift);

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

     Convert tni to radix 1e9. The high order digit of tni will 
     vary between 10 and 101. About one out of 256 conversions 
     will produce a tni with twenty-one decimals, that is, a tni 
     with a high order digit greater than 99. The three radix-1e9 
     digits of tni represent in an easily converted form the 20 
     (or occasionally 21) digits of the decimal form of the 
     result, accurate to within two ulps. After converting, set 
     oops21 if there are 21 decimals in tni, to aid in adjusting 
     tx. We augment tx by oops21 to allow for the fact that the 
     approximation computed by _l10 was off by 1. We also compute 
     the index of the leading nonzero digit in the buffer tnf, by 
     subtracting oops21 from 7 (the usual index). */ 
   
 _spdiv8(tni, tni + 1, tni + 2);
 oops21 = tni[0] > 99;
 tx += oops21;

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

     If the number is subnormal and _ieeerounding is CVT_ROUND,
     we have to go over much of what 
     we have already done, but with a revised argument. The reason 
     for this is that everything done so far has assumed that the 
     maximum significance for a number is 17. For subnormal 
     numbers, when rounding to roundest, the maximum significance
     varies with the number of 
     bits in the significand (which ranges between 1 and 52). The 
     smallest subnormal number 00000000 00000001, has one bit in 
     it, and, to 17 digits, has the value 4.940656458412465e-324. 
     The next larger subnormal number 00000000 00000002, has the 
     value 9.8813129168249308e-324; in other words, is twice as 
     large. To avoid giving the appearance of spurious accuracy, 
     we gradually reduce the maximum significance of subnormal 
     numbers. With full accuracy asked for, we will give these two 
     numbers as 5e-324 and 1e-323, respectively, since these
     short numbers are the roundest that will yield the desired
     internal form. What we have done 
     so far has, for subnormal numbers, almost the sole purpose of 
     allowing us to determine whether the exponent of 10, tx, is 
     accurate or off by one. Now we go through essentially the 
     same sequence as before, but with a significand which is not 
     left-adjusted, only shifted left 12 (the length of the 
     combined sign and exponent fields), possibly leaving some 
     leading zeros in the significand. We also save some computing 
     because the scaling power of 10 is constant (327), as well as 
     are the values tp (1087), np (-1022), tnp (1087 + -1022 = 
     65), and tnishift (96 - 65 = 31). This strategy allows us 
     gradually to diminish the significance of subnormal numbers. 
     */ 

 if (normal < (_ieeerounding == CVT_ROUND))
     {ts[0] = 0x9a65406d;   /* Approximation */
      ts[1] = 0x44a5c903;   /*      to       */
      ts[2] = 0x75000000;   /*     1e327     */
      maketnh = ts[0];  
      ns[0] = fraction[0] << 12 | fraction[1] >> 20;
      ns[1] = fraction[1] << 12;  
      _mulu32(ts, ns);
      tnp = 65;    /* Needed later in forming tnh from maketnh. */
      tni[2] = ts[2] >> 31 | (ts[1] << 1);
      tni[1] = ts[1] >> 31 | (ts[0] << 1);
      tni[0] = ts[0] >> 31;
      sticky |= 0 != ts[2] << 1; 
      _spdiv8(tni, tni + 1, tni + 2);
     }

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

     Convert the radix 1e9 integers in tni into characters and 
     place these characters in the result buffer. Note that the 
     _rad2b function, even though it will stop early when only 
     zeros are left to store, nonetheless insures that the result 
     field contains trailing "0"s. */

 _rad2b (tnf + 18, tni[2]);
 _rad2b (tnf + 9, tni[1]);
 _rad2b (tnf, tni[0]);
 #if DEBUG
 printf ("tnfA is %s\n", tnf);
 #endif

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

     Find the index of the leading nonzero digit in tnf. */

 for (ldi = 7 - oops21; tnf[ldi] == '0'; ldi++);

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

     Determine dtri, the index in tnf of the digit position which 
     will be rounded. The sole use made of eflag occurs here, and 
     it has the effect (if eflag is 0, that is, F  format is 
     desired) of adding tx + 1 to ldi, thereby locating the digit 
     position to the right of the decimal point. Since tx + 1 may 
     be negative, the result of this addition could be negative. 
     Adding ndigits - 1 to this sum locates the index in tnf of 
     the digit to be rounded. We constrain this value not to 
     exceed 26, because we never give more than 17 digits 
     of precision (further digits if asked for will be made 
     zeros).  */

 dtrir = dtri = min(26, ldi + (ndigits - 1) + ((tx + 1) & (eflag-1)));

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

     We use TOROUNDEST rounding if _ieeerounding is CVT_ROUND. */

 if (CVT_ROUND == _ieeerounding)

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

     TOROUNDEST rounding may be understood by studying the 
     following cases:

     tni[2] = 123456789 123456489 123456500 123457500
     tnh =         2345      2345      2345      2345
     carry =          0         0         0         0
     borrow =         0         0         0         0
     result = 123457000 123456000 123456000 123458000

     tni[2] = 123451234 123456789 123456789 123446789
     tnh =         5678      4567      9876      9876
     carry =          0         1         1         1
     borrow =         1         0         1         1
     result = 123450000 123460000 123460000 123440000  */

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

     Find tnh, half the interval between this tni and the one that 
     would have resulted from the next higher argument.  
     Intuitively, tnh may be obtained by multiplying ts by an ns 
     which had a single 1 bit in position 65. In practice, we 
     obtain tnh by shifting maketnh, the first element of the 
     scaling power of ten, to the right by the amount 86 - tnp, 
     if the argument is normal, or 86 - (tnp + 1) if the argument
     is subnormal (because subnormal numbers are effectively one
     place to the right of normal numbers, lacking the leading
     implied one bit). */
 
     {tnh = maketnh >> (86 - (tnp + !normal));
     
/****************************************************************  

     Abate tnh by one part in 1024, in order that two adjacent 
     numbers don't both claim the same round number between them. 
     The value of tnh will range roughly between 1100 and 11300. 
     The amount of abatement has been established empirically to 
     satisfy the conflicting requirements of maximizing accuracy 
     while avoiding overlap. */
   
      tnh -= tnh >> 10;
       
/****************************************************************

     Find whichp, the index in p of the smallest power of ten 
     greater than tnh. This we know empirically to be one of 
     1000, 10000, or 100000. */
   
      whichp = 2 + (tnh > 1000) + (tnh > 10000);

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

     We make use of p[whichp] several times later, so record its 
     value in tnt. We create dtrir, as well. The digit to round
     in for rounding to roundest is determined by p[whichp]; if
     this is 1000, we round in tnf[23]; if it is 10000, we round
     in tnf[22]; if it is 100000, we round in tnf[21]. */

      tnt = p[whichp];
      dtrir = 25 - whichp;
     
/****************************************************************

     Find the remainder rem upon division of the last part of tni, 
     that is tni[2], by tnt. */ 
   
      rem = tni[2] % tnt;
     
/****************************************************************

     Find the "carry/borrow status" of the remainder. Carry is set 
     if tnt is less than rem + tnh. If the argument is not an 
     integral power of two, borrow is set if rem is less than tnh; 
     if it is an integral power of two, borrow is set if rem is 
     less than tnh/2. This is because the difference between a tni 
     that is an exact power of two and the next smaller number is 
     only half that between other numbers having the same 
     exponent. */
 
      carry = tnt < (rem + tnh);
      borrow = rem < (tnh >> poweroftwo);

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

     We round if either or both of carry or borrow is set. */

      if (carry || borrow)
          

/*****************************************************************
      
     Determine the rounding amount r. This is 1 if there is carry, 
     but no borrow, or if there is both carry and borrow and the 
     last digit of quo is odd. */ 

       r = (carry > borrow) || carry && borrow && (1 & tnf[dtrir]);
     
/*****************************************************************

     If there is neither carry nor borrow add 1 to dtrir, 
     indicating that we are now concerned to round in the next 
     digit position. For an explanation of the formation and 
     meaning of parity, guard, sticky, sum, and r, see the 
     comments below in connection with rounding TONEAREST. */ 

      else
          {dtrir++;
           parity = 1 & tnf[dtrir]; 
           guard = tnf[dtrir + 1] - '0';
           for (i = dtrir + 2; i < 27; i++)
               if (tnf[i] != '0')
                    {sticky = 1;
                     break;
                    }
           sum = (guard << 2) + (sticky << 1) + parity;
           r = sum > 20;
          }

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

     Set all the characters in tnf after tnf[dtrir] to '0'. Then,
     if r (the rounding amount) is 1, do the required rounding
     to roundest. */

      for (aldi = 6; tnf[aldi] == '0'; aldi++);
      for (i = dtrir + 1; i < 27; i++)
          tnf[i] = '0';
      for (i = dtrir; r; i--)
          {tnf[i]++;
           if (tnf[i] > '9')
                tnf[i] = '0';
           else break;
          }
      tx += i < aldi;
      aldi -= i < aldi;
     }
#if DEBUG
printf ("tnfB is %s\n", tnf);
#endif

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

     Round tnf[dtri] if necessary. In order for this to be
     possible dtri must be less than dtrir (to avoid double
     rounding). */

 if ((dtri < dtrir) || (CVT_ROUND != _ieeerounding))
     {

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

     However, we want to avoid rounding if dtri is to the left of
     ldi (it could even be negative, and might cause addressing
     problems). */

      for (ldi = 6; tnf[ldi] == '0'; ldi++);
      if (dtri >= ldi)
           {
/****************************************************************

     IEEE rounding can be described as follows: 

     In TONEAREST rounding, one is added to tnf[dtri] if guard, 
     tnf[dtri + 1], is 6 or greater, or if it is 5 and sticky is 
     1, or if it is 5 and sticky is zero and dtri is odd. The 
     value sticky is 1 if any of the digits following guard are 
     nonzero.

     In UPWARD rounding, one is added to tnf[dtri] if the sign is 
     zero and either guard or sticky is nonzero. 

     In DOWNWARD rounding, one is added to tnf[dtri] if the sign 
     is one and either guard or sticky is nonzero. 

     In TOWARDZERO rounding, one is never added.

     If we multiply the guard digit by four, and the sticky value 
     by two, and let the parity of tnf[dtri] be 0 or 1 according 
     to whether it is even or odd, we can make the following 
     tables: 

     guard 4*guard     0       2   = 2*sticky    TO      UPWARD&+
     digit  digit    0   1   0   1 = parity    NEAREST  DOWNWARD&-

       0      0      0   1   2   3             0 0 0 0    0 0 1 1
       1      4      4   5   6   7             0 0 0 0    1 1 1 1
       2      8      8   9  10  11             0 0 0 0    1 1 1 1
       3     12     12  13  14  15             0 0 0 0    1 1 1 1
       4     16     16  17  18  19             0 0 0 0    1 1 1 1
       5     20     20  21  22  23             0 1 1 1    1 1 1 1
       6     24     24  25  26  27             1 1 1 1    1 1 1 1
       7     28     28  29  30  31             1 1 1 1    1 1 1 1
       8     32     32  33  34  35             1 1 1 1    1 1 1 1
       9     36     36  37  38  39             1 1 1 1    1 1 1 1

     The entries in the tables TO/NEAREST or UPWARD&+/DOWNWARD&- 
     give the amount added to tnf[dtri]. This suggests that the 
     rounding amount can be determined by summing four times the 
     value of the guard digit, twice the value of sticky, and the 
     parity of tnf[dtri]. If rounding is TONEAREST and the sum 
     exceeds 20, 1 will be added; if rounding is UPWARD&+ or 
     DOWNWARD&-, and the sum exceeds 1, 1 will be added. */ 

	   rm = fptestround();
           parity = 1 & tnf[dtri]; /* 0 or 1 */
           guard = tnf[dtri + 1] - '0';
           for (i =  dtri + 2; i < 27; i++)
               if (tnf[i] != '0')
                    {sticky = 1;
                     break;
                    } 
           sum = (guard << 2) + (sticky << 1) + parity;
           ud = ((rm == UPWARD)&& !*sign)
                || ((rm == DOWNWARD) && *sign);
           nr = (rm == TONEAREST) || (_ieeerounding == CVT_ROUND);
#if DEBUG
printf("ud is %d, nr is %d\n", ud,nr);
printf("sum is %d\n", sum);
printf("parity is %d, guard is %d, sticky is %d\n", parity,guard,sticky);
printf("rm is %d\n", rm);
#endif
           r = (nr && sum > 20) || (ud && sum > 1);
	   for (i = dtri + 1; i < 27; i++)
		tnf[i] = '0';
	   for (aldi = 6; tnf[aldi] == '0'; aldi++);

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

     We do the actual rounding, at long last. In order to detect 
     the case when carry causes a new decimal digit position to be 
     formed (that is, when 09...9 becomes 10...0) we use aldi, the
     actual leading digit position. For normal numbers, this is
     the same as ldi; for subnormal numbers it may not be. */

           for (i = dtri; r; i--)
                {tnf[i]++;
                 if (tnf[i] > '9')           /* with carry. */
                     tnf[i] = '0';
                 else break;
                } 

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

     Augment tx if carry has caused a new decimal place to be 
     added. */

            tx += i < aldi;
	    aldi -= i < aldi;
#if DEBUG
printf ("tnfC is %s\n", tnf);
#endif
           }

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

     The digit to round is to the left of the leading digit. This
     means that we have no significance to represent. If E format
     is specified we return an empty vector. If F format is
     specified we return a one element vector containing the
     character digit zero. */

     else
	   {*decpt = tx + 1;
	    if (eflag)
		 return ("");
            else 
		 return ("0");
           }
      }
     
/*****************************************************************

     We locate the first nonzero digit at or to the left of
     tnf[dtri], and store a binary zero in the position following
     this, as a string terminator. We also locate the most
     significant digit in order to return its address as the
     formal result. However, two circumstances must be attended
     to which would interfere with this normal situation. The
     first abnormal case is where E format is asked for, with
     no significant digits wanted. In this case, we return a
     null string, by storing the binary zero (arbitrarily) at
     tnf[0], and by setting j to 0, thereby giving the address of
     the string terminator as the result, which effectively
     designates the null string. The second abnormal case is where
     F format is wanted, and all the precision is to the right of
     the last position retained. In this case, we store the binary
     zero in tnf[1], and set j to 0, thus returning a character
     string consisting of one character '0'. */

      if (eflag && !ndigits)
          {tnf [0] = 0;
           j = 0;
          }
      else
          if (!eflag && (dtri < aldi) && (tnf[dtri] == '0'))
               {tnf[1] = 0;
                j = 0;
               }
          else {
                tnf[dtri + 1] = 0;  
		for (j = 6; tnf[j] == '0'; j++);
               } 
#if DEBUG
printf ("tnfD is %s\n", tnf);
printf ("j is %d\n", j);
#endif

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

     Store the location of the decimal point as a side effect. */

 *decpt = tx + 1;         /* decpt is one more than 10's exponent */

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

     Return the address of the most significant digit in the 
     character vector result. */

 return (tnf + j);
}


