/* numLongDiv: Divide two NumberT's using long division

30sep96 MJG -- Not finished, not tested

**
** return quo = num / den  and  rem = num MOD den
**
** num (in): numerator or dividend, must be twice the size of den
** den (in): denominator or divisor
** quo (out): quotient
** rem (out): remainder
**
** The quotient and remainder have n - m units
** where n is the number of units of the dividend
** where m is the number of units of the divisor.
**
** It is assumed that the length of num is NUMBER_LENGTH*2
** and that the lengths of den, quo, and rem are NUMBER_LENGTH.
**
** See D.E.Knuth, The Art of Computer Programming, 2nd ed, 1973, section 4.3.1

NOTE that Knuth's indices for positional notation of big numbers is
backwards from mine.  For Knuth, the most significant nonzero unit
has index 0 or 1, whereas in my notation, index 0 refers to the least
significant unit.  In my case, the most significant nonzero unit must
be searched for explicitly.  (Actually, in any implementation of
Knuth's routine, the most significant nonzero unit would have to be
searched for also.)  Also note that Knuth's positional indices run in
the opposite direction from my positional indices.

*/
void
numLongDiv(const NumberT num, const NumberT den, NumberT quo)
{
  int lnz_num;       /* index of leading non-zero unit of numerator */
  int lnz_den;       /* index of leading non-zero unit of denominator */
  unsigned long dl;  /* normalization value */
  NumberT un;        /* temp space numerator */
  NumberT vn;        /* temp space denominator */
  NumberT big_tmp;   /* temp space */
  int ji;            /* loop variable */
  unsigned long qh;  /* trial quotient */
  int borrow;        /* borrow flag */

#ifdef VERBOSE
  printf("num = ");  numPrint(num, NUMBER_LENGTH*2);
  printf("den = ");  numPrint(den, NUMBER_LENGTH);
#endif

  /* Allocate space for temp space variables */
  numAlloc(&vn, NUMBER_LENGTH);
  numAlloc(&un, NUMBER_LENGTH * 2 + 1);
  numAlloc(&big_tmp, NUMBER_LENGTH * 2 + 1);

  memset(quo, 0, sizeof(NumberPrimT) * NUMBER_LENGTH);

  /* Find leading nonzero unit for denominator */
  for(lnz_num=NUMBER_LENGTH*2; (lnz_num >= 0) && (!num[lnz_num]); lnz_num--)
    ;
  for(lnz_den=NUMBER_LENGTH-1; (lnz_den >= 0) && (!den[lnz_den]); lnz_den--)
    ;

  /* D1: Normalize */
  dl = NUMBER_BASE / (den[lnz_den] + 1);
  dl = 1;

  numMultPrimUnit(dl, num, un, NUMBER_LENGTH*2, 0);
  numMultPrimUnit(dl, den, vn, NUMBER_LENGTH, 0);

  /* Find leading nonzero units */
  lnz_num ++;
  for(lnz_den=NUMBER_LENGTH-1; (lnz_den >= 0) && (!vn[lnz_den]); lnz_den--)
    ;

  /* Steps D2 - D7 */
  for(/*D2*/ ji = lnz_num; /*D7*/ ji > 0; ji --) {

    /* -- D3: Calculate qh -- */

    if(un[ji] == vn[lnz_den]) {
      qh = NUMBER_BASE - 1;
    } else {
      qh = (un[ji] * NUMBER_BASE + un[ji-1]) / vn[lnz_den];
    }

    while(vn[lnz_den-1]*qh >
            (un[ji]*NUMBER_BASE + un[ji-1] - qh*vn[lnz_den])*NUMBER_BASE
             + un[ji-2])
    {
      qh --;
    }


    /* -- D4: Multiply and subtract -- */

multiply:
    /* Multiply to form subtrahend */
    memset(big_tmp, 0, sizeof(NumberPrimT) * (NUMBER_LENGTH * 2 + 1));
    numMultPrimUnit(qh, vn, big_tmp, NUMBER_LENGTH, 0);

    /* Subtract */
    if(numCompare(un, big_tmp, NUMBER_LENGTH*2+1) < 0) {
      /* Subtract would give a negative: */
      /* Find complementary subtraction
      ** i.e. add NUMBER_BASE^(NUMBER_LENGTH+1) to the minuend
      ** (where NUMBER_LENGTH is the length of the denominator)
      */
      un[NUMBER_LENGTH] ++;

      /* Remember "borrow" */
      borrow = 1;
    }
    numSubPrim(un, big_tmp, un, NUMBER_LENGTH*2+1);


    /* -- D5: -- */
    quo[ji-1] = qh;


    /* -- D6: Test remainder  and  Add back -- */

    if(borrow) {
      quo[ji-1] --;
      borrow = 0;
      goto multiply;
    }
  }

  numFree(&vn);
  numFree(&un);
  numFree(&big_tmp);
}