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

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

/* TENTOTHE  --  COMPUTE A POWER OF TEN

Returns ten to the power "t", in the first 72 bits of three words. 
The values of t for which this program has been tested range from 
-352 to 352, inclusive.

The key to correct conversions between external decimal and 
internal floating-point binary forms is a correct scale factor, 
10*scale. Jerome Coonen's Ph. D. dissertation gives an acccurate 
analysis of the requirements for such conversions, but his 
Algorithm Q, for obtaining a nonnegative power of 10 in double 
format, fails to give exact results for twenty arguments out of 
340. For ten of these, his results are too low by one ulp, and 
for the other ten his results are too high by one ulp.  

The algorithm given here gives results, accurate from 69 to 73 
bits, for all arguments in the range -352 to 352, inclusive. It 
obtains negative powers of 10 as well as positive powers, and thus 
allows scaling by negative powers of ten by using an extended 
precision multiplication algorithm, rather than the more 
complicated extended precision divide algorithm that Coonen's 
method would require. The two negative powers of ten at the 
beginning of the table and the three positive powers of ten at the 
end are "accurate table" values. An accurate table value is one 
which is particularly close to an exact value. For example, the 
value for 1e55, although given to only 64 bits, is within 1 part 
in 2*76 of being exact. Using accurate table values to extend a 
table makes it possible to get better than 64-bit accuracy with 64 
bit values. The values for powers 1 through 27 are exact.

The powers of ten are kept as left justified 64-bit values. 
"Left-justified" means the same as "normalized", and it means 
that the leading bit of the value is 1. In the tabulation below 
an indication is given which shows if the value is exact, or if 
not, what the values of the next several elided bits are. */

#define COLS 2                  /* Number of columns in table*/
#define HIGHBITONE 0x80000000   /* Smallest no. with bit 0=1 */
#define ROWS 33                 /* Number of powers in table */
#define DENSEMIN 0              /* Minimum dense argument    */
#define DENSEMAX 27             /* Maximum dense argument    */
#define DENSESTART 2            /* Index of DENSEMIN in ptx  */
#define DENSEEND DENSESTART+DENSEMAX-DENSEMIN
                                /* Index of DENSEMAX in ptx  */

static long pt[ROWS][COLS] =    /* power of 10   elided bits */
     {0xcd42a113, 0x46f34f7d,   /*    -352          0092     */
      0x9fd561f1, 0xfd0f9bd4,   /*    -145       (3)feb6     */
      0x80000000, 0x00000000,   /*       0         exact     */
      0xa0000000, 0x00000000,   /*       1         exact     */
      0xc8000000, 0x00000000,   /*       2         exact     */
      0xfa000000, 0x00000000,   /*       3         exact     */
      0x9c400000, 0x00000000,   /*       4         exact     */
      0xc3500000, 0x00000000,   /*       5         exact     */
      0xf4240000, 0x00000000,   /*       6         exact     */
      0x98968000, 0x00000000,   /*       7         exact     */
      0xbebc2000, 0x00000000,   /*       8         exact     */
      0xee6b2800, 0x00000000,   /*       9         exact     */
      0x9502f900, 0x00000000,   /*      10         exact     */
      0xba43b740, 0x00000000,   /*      11         exact     */
      0xe8d4a510, 0x00000000,   /*      12         exact     */
      0x9184e72a, 0x00000000,   /*      13         exact     */
      0xb5e620f4, 0x80000000,   /*      14         exact     */
      0xe35fa931, 0xa0000000,   /*      15         exact     */
      0x8e1bc9bf, 0x04000000,   /*      16         exact     */
      0xb1a2bc2e, 0xc5000000,   /*      17         exact     */
      0xde0b6b3a, 0x76400000,   /*      18         exact     */
      0x8ac72304, 0x89e80000,   /*      19         exact     */
      0xad78ebc5, 0xac620000,   /*      20         exact     */
      0xd8d726b7, 0x177a8000,   /*      21         exact     */
      0x87867832, 0x6eac9000,   /*      22         exact     */
      0xa968163f, 0x0a57b400,   /*      23         exact     */
      0xd3c21bce, 0xcceda100,   /*      24         exact     */
      0x84595161, 0x401484a0,   /*      25         exact     */
      0xa56fa5b9, 0x9019a5c8,   /*      26         exact     */
      0xcecb8f27, 0xf4200f3a,   /*      27         exact     */
      0xd0cf4b50, 0xcfe20766,   /*      55       (5)fff4     */    
      0xc1a12d2f, 0xc3978937,   /*     123          0052     */
      0xc1d4ce1f, 0x63f57d73};  /*     269       (2)fd1c     */

/* With this table of values it is possible to compute any power 
of 10 between -352 and 352 with at most 4 extended-precision 
multiplies, and only one power requires that many (-146). All of 
the commonly occurring powers can be computed simply by selection 
from the table. */ 

/* The algorithm requires the presence of a vector having one 
element for each row in the table, giving the power of 10 
associated with that row. */

static short ptx[ROWS] = {-352, -145, /* start dense entries */
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 
19, 20, 21, 22, 23, 24, 25, 26, 27,
/* end dense entries */ 55, 123, 269};
/* With this table of values it is possible to compute any power 
of 10 between -352 and 352 with at most 8 extended-precision 
multiplies. */ 

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

void _TenToThe (t, w)    /* Compute the t-th power of 10.   */

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

long t;                 /* Argument, a signed integer.     */
unsigned long w[3];     /* Result.                         */

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

{register int i;        /* Index of 1st element of ptx  
                           t is greater than or equal to.  */
 unsigned long z[3];    /* Developing product.             */

 /************************* ALGORITHM **************************/
       
/* Because fairly small powers of 10 are expected to be common, 
the first thing done is to check to see if t is one of those in 
the dense set. If it is, row t of table pt is returned directly 
as the result. */

if ((t >= DENSEMIN) && (t <= DENSEMAX))
     {i = DENSESTART + t - DENSEMIN;
      w[0] = pt[i][0];
      w[1] = pt[i][1];
      w[2] = 0;
      return;
     }

/* The row index of the initial value of the result is now 
determined. This initial value is the index of the largest 
element of ptx which t is greater than or equal to.  The search 
begins with one of two possible values: if t is negative, the 
search begins with the value in ptx just below the start of the 
dense set of entries. Since there are only two entries below this 
start, the choice is easy to make. Otherwise, the search begins 
with the largest value in ptx. Since there are four possibilities 
here, the search is a simple scan of the four possibilities. The 
search begins with the largest value. */ 

if (t<0) i = (t >= ptx[1]); /* Zeroth or first. */
else
     {i = ROWS-1; /* last value is 1st choice */
      while (t < ptx[i]) --i; /* reduce i until ptx[i] greater 
               or equal. At most four comparisons will be made. */
     } 

/* The initial value of the result can now be obtained. */

z[0] = w[0] = pt[i][0];
z[1] = w[1] = pt[i][1];
z[2] = w[2] = 0;

/* Reduce t by the size of the exponent just used. If this causes 
t to equal 0, return w. */

t -= ptx[i]; 
if (!t) return;

/* The preceding was the easy part, and handled those cases where 
no multiplication was required. Now starts the part where a power 
of ten is needed which is not in the table pt, but must be 
calculated by taking the product of two or more powers of ten. 
The examples below show how the desired exponent is produced as 
the sum of other exponents: 

  desired exponent   summand exponents
         23                 23        
         28                27 1       
         99              55 27 17     
        122              55 55 12     
        268             123 123 22    
         -1            -145 123 22        
       -146          -352 123 55 27 1 

The idea is to start with the smallest power of ten which t is 
greater than or equal to. We've already done this, and this is 
the current value of z, and we've reduced t by the amount of the 
exponent used. We now find the next smallest power of ten which 
the new t is greater than or equal to. We select this entry from 
table pt, and multiply z by it, forming a new z. We reduce t by 
the amount of the exponent just used, and if this makes t zero,
we're through.  Otherwise we repeat the process until t becomes 
zero.  
 
Before starting the multiplication process, however, there is one 
possible adjustment that must be made. If the original value of t 
was negative, i now points to one of the rows of pt containing a 
negative power of 10. Since the strategy with negative values is 
to obtain it using only one negative power of ten, and from then 
on using only positive powers of ten, we have to change i so that 
it points to the last row of ptx in order that that our searches 
lead us from the largest positive power of ten down through the 
smaller positive powers of ten. */

if (i < DENSESTART) i = ROWS - 1;

/* At last we're ready for the multiplication process. */

while (t)
     {if (t >= ptx[i]) 
          {_mulu32 (z,pt[i]); 
           t -= ptx[i];
          }
      else
          {--i;
           if (t <= DENSEMAX) i = t + DENSESTART - DENSEMIN;
          }
     }
/* In words, what the above sequence does is roughly as follows: 
as long as t isn't zero we'll do either of two things: if t is 
greater than or equal to the exponent of ten pointed at by i, 
we'll use that power of ten to form a new z from the old z, and 
we'll reduce t by that exponent value; if t is less than that 
exponent, we'll either reduce i by one so that it points at a 
smaller exponent, or, if we're using the long table, and t is one 
of the dense set of exponent values, we'll make i point right away 
at the entry for t, rather than looping down through these dense 
values one at a time. */ 

/* At this point, the three-longword z contains an accurate value 
of ten to the original t, accurate to between 69 and 73 bits. 
Now we have to insure that z is normalized, or left-adjusted, and 
then we'll truncate the last element of z by setting to zero all 
but its first 8 bits. */

/* Normalizing requires a triple-precision shift. */

while (z[0] < HIGHBITONE)
 {z[0] = z[0]<<1 | z[1]>>31;  
  z[1] = z[1]<<1 | z[2]>>31;
  z[2] <<= 1;
 }

/* Truncation simply ands the first 8 bits of z[2] in. */

z[2] &= 0xff000000;

/* The result is w. */

w[0] = z[0];     
w[1] = z[1];
w[2] = z[2];
}

/********************** END OF ALGORITHM ************************/
