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

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

/* SPECIAL CONVERSION TO BASE 1E9

The argument is three unsigned longwords, considered as the three 
"digits" of a number in radix two to the 32. The result is three 
unsigned longwords, representing the same number in radix 1e9.  

The algorithm is highly specialized toward this particular 
problem, and makes extensive use of the properties of various 
constants, in order to obtain efficiency.

For information concerning the range of values which the argument 
might take on, see the comments to function dfb9.

A straightforward implementation of this function would use an 
algorithm in which three successive divisions of two longwords by 
one longword occurred. Because divisions are relatively 
expensive, I attempt to gain efficiency by replacing them 
wherever possible by multiplications, additions, and shifts. The 
basis for the technique uses the following observation: consider 
the table below, in which the argument and result are shown 
side-by-side.  
                                     
  argument         result            
    1 0 0    18 446744073 709551616     (hi)
    0 1 0     0         4 294967296     (mid)
    0 0 1     0         0         1     (lo)

Because of the way the calling program dfb9 scales its argument, 
I know that its leading element will be between 0 and 5. Thus, 
to obtain the result summand corresponding to the first digit of 
the argument, it suffices to select the appropriate entry from a 
six-row by three-column table that contains the 0 through 5 
multiples of (hi). Note that (hi) is the base 1e9 representation 
of the 64th power of two.

Because (lo) is equal to 1, the result summand corresponding to 
the last digit of the argument can be formed by, for each whole 
multiple of 1e9 in the last digit, up to a maximum of four times, 
simultaneously adding 1 to the central result digit, and 
subtracting 1e9 from the last digit of the argument, and  then 
adding the (possibly reduced) last digit of the argument to the 
last digit of the result. 

The result summand corresponding to the middle term of the 
argument can be developed by adding the appropriate multiple, in 
radix 1e9, of (mid) to the result for each 1 bit in the binary 
representation of the middle term of the argument. It so happens 
that the first 28 multiples of (mid) have 0 for the first 
element, so only two terms have to be added for these multiples; 
the last 4 multiples have nonzero first elements, so for these 
all three terms must be added to the result.

The first term can never be very big, so carries need not be 
worried about for it. The second term can get to exceed 1e9, but 
will always be less than 2 power 32 (since the sum of all the 
middle multiples is less than 2 power 32), even with maximum 
carry into the middle result cell from the last result cell, so 
carry beyond 1e9 in the middle term, which can easily occur, need 
only be checked for at the very end. The last result term, 
however, can equal or exceed 2 power 32, so each time the last 
result term is added to a test must be made for carry beyond 1e9. 
*/

static unsigned long hitable[6][3] = 
     {         0,         0,         0,
              18, 446744073, 709551616,
              36, 893488147, 419103232,
              55, 340232221, 128654848,
              73, 786976294, 838206464,
              92, 233720368, 547758080
     };

static unsigned long midtbl[28][2] =
     {         4, 294967296,
               8, 589934592,
              17, 179869184,
              34, 359738368,
              68, 719476736,
             137, 438953472,
             274, 877906944,
             549, 755813888,
            1099, 511627776,
            2199,  23255552,
            4398,  46511104,
            8796,  93022208,
           17592, 186044416,
           35184, 372088832,
           70368, 744177664,
          140737, 488355328,
          281474, 976710656,
          562949, 953421312,
         1125899, 906842624,
         2251799, 813685248,
         4503599, 627370496,
         9007199, 254740992,
        18014398, 509481984,
        36028797,  18963968,
        72057594,  37927936,
       144115188,  75855872,
       288230376, 151711744,
       576460752, 303423488
     };

static unsigned long midtbl2[4][3] =
     {         1, 152921504, 606846976,
               2, 305843009, 213693952,
               4, 611686018, 427387904,
               9, 223372036, 854775808
     };

#define ONEBILL 1000000000

void _spdiv8 (tni0, tni1, tni2)
unsigned long *tni0, *tni1, *tni2;
{unsigned long z[3];
 unsigned short i;

 z[0] = hitable[*tni0][0];
 z[1] = hitable[*tni0][1];
 z[2] = hitable[*tni0][2];

 for (i = 0; i < 28; i++)
     {if (1 << i & *tni1)
          {z[2] += midtbl[i][1];
           if (z[2] >= ONEBILL)
                 {z[2] -= ONEBILL;
                  z[1]++;
                 }
           z[1] += midtbl[i][0];
          }
     }

 for (i = 28; i < 32; i++)
     {if (1 << i & *tni1)
          {z[2] += midtbl2[i-28][2];
           if (z[2] >= ONEBILL)
                 {z[2] -= ONEBILL;
                  z[1]++;
                 }
           z[1] += midtbl2[i-28][1];
           z[0] += midtbl2[i-28][0];
          }
     }

 while (*tni2 >= ONEBILL)
     {z[1]++;
     *tni2 -= ONEBILL;
     }
 z[2] += *tni2;
 if (z[2] >= ONEBILL)
     {z[1]++;
      z[2] -= ONEBILL;
     }

 while (z[1] >= ONEBILL)
     {z[0]++;
      z[1] -= ONEBILL;
     }

 *tni0 = z[0];
 *tni1 = z[1];
 *tni2 = z[2];
}

