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

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


/* DIVIDE TEN SHORTWORDS BY FIVE SHORTWORDS

The arguments are a four-shortword dividend x, and a six-
shortword divisor, y. The dividend and divisor are assumed to be 
already normalized. This assumption makes it unnecessary for this 
function to worry about normalizing. The result is also in y and 
is a six-shortword quotient resulting from dividing x by y. 

The basis for the algorithm is Algorithm D in section 4.3.1 of 
D.E. Knuth's 'Seminumerical Algorithms'. The discussion of the 
algorithm which follows is largely paraphrased from Knuth's 
description of the algorithm. */

#define UL unsigned long
#define US unsigned short

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

void _dknu (x, y)

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

US x[4];           /* dividend */
US y[6];           /* divisor and quotient */

{

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

 UL c; /* difference of first two terms of current dividend and 
          first partial product */
 US j; /* induction variable */
 US qh; /* trial quotient */
 US t[6]; /* remembers current dividend */
 US m[6]; /* full product, sum of partial products */
 UL o;    /* carry/borrow indicator */ 
 US u[10]; /* working dividend */
 US z[5]; /* working divisor */
 short i; /* induction variable */
 union v
     {                /* To allow data item to be      */
      UL a; US b[2];  /*  treated as either one        */
     };               /*   longword or two shortwords. */
 union v n; /* short c catenated with third term of current
               dividend */
 union v d; /* partial product */
 union v p; /* partial product */
 union v q; /* partial product */
 union v r; /* partial product */
 union v s; /* partial product */
 union v w; /* first two terms of current dividend */
 
 /************************* ALGORITHM **************************/
 
 /*A new dividend u is obtained from x by extending it on the left 
 with a leading zero, and on the right with five trailing zeros. 
 The divisor y is copied into z so that y may be used to store the 
 result. Only the first five shortwords of the divisor are copied 
 into z, since we know that the last shortword is zero. Knuth's 
 algorithm D provides a normalizing step, but this function 
 assumes the arguments are already normalized. 

                   ________________y0_y1_y2_y3_y4_
   z0 z1 z2 z3 z4 ) u0 u1 u2 u3 u4 u5 u6 u7 u8 u9
                   <-- z * y0 --->
                  _______________
                      <--- r --->

 Each step of the division process can be viewed as the 
 determination of the largest yi which when used to multiply z
 gives a product which is less than or equal to the current 
 portion of interest of the dividend. The complete quotient is the 
 list of all the y's. The art lies in determining at each step an 
 initial estimate for yi, called qh below, and then correcting this 
 estimate if necessary. */

 u[0] = 0;                     /* u is 0, x, 0 0 0 0 0 */
 u[1] = x[0]; u[2] = x[1]; u[3] = x[2]; u[4] = x[3];
 u[5] = u[6] = u[7] = u[8] = u[9] = 0;
 z[0] = y[0]; z[1] = y[1]; z[2] = y[2]; z[3] = y[3]; z[4] = y[4];
#if DEBUG
printf ("u is %d %d %d %d\r\n", u[1], u[2], u[3], u[4]);
printf ("z is %d %d %d %d %d\r\n", z[0], z[1], z[2], z[3], z[4]);
#endif

 for (j = 0; j <= 4; j++)      /* induction on j */
     {
      /*The initial estimate of qh is guaranteed to be either 
      accurate or high by just one. It is formed by taking either 
      65535 (one less than the number system base, since the 
      elements are all short words) if the leading element of the 
      divisor and the leading element of the current dividend 
      agree, or the quotient of the first two elements of the 
      current dividend and the leading element of the divisor. 
      Before accepting this value of qh, we test its product with 
      z[1] against the difference of the first three digits of the 
      divisor, less qh times z[0] times the base 65536. If qh*z[1] 
      is greater, we decrease qh by 1 and repeat the test. This 
      determines most of the cases when qh is one too large, and 
      ALL cases where it is two too large. The justification for 
      this step is given in Knuth's Section 4.3.1; see Theorem A 
      and step D3 of Algorithm D. */

      w.b[0] = u[j];
      w.b[1] = u[j+1];

      if (u[j] == z[0])        /* Depending on test, */
          qh = 65535;          /* qh is maximum      */
      else                     /* or trial quotient. */
          qh = (((UL)u[j] << 16) + u[j+1]) / z[0];
#if DEBUG
printf ("j1 is %d, qh1 is %d\r\n", j, qh);
#endif

      for (;;)                 /* May have to abate qh. */
          {
           s.a = z[0] * qh;    /* s.a needed later. */
           c = w.a - s.a;
           if (c >= 65536) break;
           n.b[0] = (US) c;
           n.b[1] = u[j+2];
           if (z[1] * qh <= n.a) break;
           qh--;
#if DEBUG
printf ("qh1a is %d\r\n", qh);
#endif
          }

      /* We make a copy, t, of the current portion of the dividend 
      for later use, since we are about to subtract from this 
      portion of the dividend in order to form the next dividend. 
      */ 

      for (i = 0; i <= 5; i++) /* Copy current dividend. */
          t[i] = u[j + i];

      /*We multiply and subtract in order to form the new 
      dividend. The multiplication forms qh times y, giving the 
      product m, which is subtracted from the current dividend to 
      form the beginning portion of the new dividend. */

      d.a = (z[4]) ? z[4] * qh : 0; /* Develop                */
      p.a = (z[3]) ? z[3] * qh : 0; /*  remaining             */
      q.a = (z[2]) ? z[2] * qh : 0; /*   partial products     */
      r.a = z[1] * qh;              /*    (already have s.a). */

      m[5] =          d.b[1]                  ; /* Make     */
      m[4] = d.b[0] + p.b[1]                  ; /* sum of   */
      m[3] = p.b[0] + q.b[1] + (m[4] < d.b[0]); /* partial  */
      m[2] = q.b[0] + r.b[1] + (m[3] < p.b[0]); /* products */
      m[1] = r.b[0] + s.b[1] + (m[2] < q.b[0]); /* with     */
      m[0] = s.b[0]          + (m[1] < r.b[0]); /* carries. */
#if DEBUG
printf ("m is %d %d %d %d %d %d\r\n", m[0], m[1], m[2], m[3], m[4], m[5]);
#endif
     
      for (o = 0, i = 5; i >= 0; i--)
          {                /* Subtract product from dividend. */
           o = (UL)u[j + i] - ((UL)m[i] + o); /* Subtract with borrow. */
           u[j + i] = (US)o; /* Store as short. */
           o = o >> 31;    /* Borrow indicator. */
          }
#if DEBUG
printf ("uA is %d %d %d %d %d %d\r\n", 
         u[j], u[j+1], u[j+2], u[j+3], u[j+4], u[j+5]);
#endif
      /* The tentative value qh is stored as the j-th element of 
      the quotient list y. */ 

      y[j] = qh;               /* Save quotient term. */
#if DEBUG
printf ("j3 is %d, y[j] is %d\r\n", j, y[j]);
#endif

      /*Now we have to determine whether, after all, the value of 
      qh we just stored is correct. If qh is zero, we can't be too 
      large. Otherwise, we look from the left through t and the 
      new dividend until we find corresponding elements that are 
      not equal. When we find one (as we must, since qh is 
      nonzero) we compare them, and if the element from t is 
      smaller it means that the subtraction step went negative, 
      and so z[j] must be reduced by one, and z must be added to 
      the current dividend. It isn't necessary to bother adding 
      the carry (if any) to the leading term of the dividend, 
      because we make no further use of this term. */ 

      if (qh)                  /* Only if qh isn't zero: */
          {
           for (i = 0; t[i] == u[j + i]; i++);
           if (t[i] < u[j + i]) /* Test first unequal terms */
               {                 /* and if old term is smaller  */
                y[j]--;          /* abate quotient term         */
#if DEBUG
printf ("j4 is %d, y[j] is %d\r\n", j, y[j]);
#endif
                for (o = 0, i = 4; i >= 0; i--)
                    {            /* and add z back to dividend. */
                     o = u[j + i + 1] + o + z[i];
                     u[j + i + 1] = (US)o;
                     o = o >> 16;    /* Carry indicator. */
                    }
                u[j] += o;  /* last carry */
#if DEBUG
printf ("uB is %d %d %d %d %d %d\r\n",
         u[j+1], u[j+2], u[j+3], u[j+4], u[j+5], u[j+6]);
#endif
               }
          } 
     }
}
#undef UL
#undef US

