/* $Header$ */

/*
 * Copyright © 1988-2004 Keith Packard and Bart Massey.
 * All Rights Reserved.  See the file COPYING in this directory
 * for licensing information.
 */

/*
 * divide.c
 *
 * perform natural number division
 */

#include	"nickle.h"

#undef CHECK
#undef DEBUG

#if defined(DEBUG) || defined(CHECK)
void
pr (Natural *n)
{
	int	i;

	for (i = 0; i < NaturalLength(n); i++)
		(void) printf ("%u ", (unsigned int) NaturalDigits(n)[i]);
	(void) putchar ('\n');
}

void
prs (char *s, Natural *n)
{
	(void) printf ("%s ", s);
	pr (n);
}
#endif

/*
 * Return shift amount needed to normalize d (MSB (d << shift) == 1)
 */
static int normalize (digit d)
{
    int	    n;
    digit   msb;

    n = 0;
    msb = 1 << (LBASE2 - 1);
    while (!(d & msb))
    {
	n++;
	d <<= 1;
    }
    return n;
}

Natural *
NaturalDivide (Natural *a, Natural *b, Natural **remp)
{
    ENTER ();
    int		    quolen, remlen;
    Natural	    *quo, *rem;
    int		    index;
    double_digit    temp;
    digit	    temp2;
    Natural	    *partial;
    int		    offset;
    digit	    dividend1;	/* first digit of dividend */
    digit	    dividend2;	/* second digit of dividend */
    digit	    dividend3;	/* third digit of dividend */
    digit	    dividend4;	/* fourth digit of dividend */
    digit	    divisor1;	/* first digit of divisor */
    digit	    divisor2;	/* second digit of divisor */
    digit	    divisor3;	/* third digit of divisor */
    digit	    divisorc;	/* combination of first two digits */
    digit	    divisorc2;	/* combination of digits 2 and 3 */
    digit	    d;
    digit	    carry;
    int		    normal;
    
    if (NaturalLess (a, b)) {
	quo = zero_natural;
	rem = a;
    } else if (oneNp (b)) {
	quo = a;
	rem = zero_natural;
    } else if ((offset = NaturalPowerOfTwo (b)) >= 0) {
	quo = NaturalRsl (a, offset);
	rem = NaturalMask (a, offset);
    } else {
	quolen = a->length - b->length + 1;
	quo = AllocNatural (quolen);
	NaturalDigits(quo)[a->length - b->length] = 0;

	partial = AllocNatural (b->length + 1);

	remlen = a->length + 1;
	rem = AllocNatural (remlen);
	NaturalCopy (a, rem);
	NaturalDigits(rem)[remlen - 1] = 0;

	divisor1 = NaturalDigits(b)[b->length - 1];
	divisor2 = 0;
	divisor3 = 0;
	if (b->length > 1)
	{
	    divisor2 = NaturalDigits(b)[b->length - 2];
	    if (b->length > 2)
		divisor3 = NaturalDigits(b)[b->length - 3];
	}
	/*
	 * Compute a scale factor to use as many bits of
	 * the divisor as possible, then scale the
	 * first and second digit of the divisor by that
	 * scale factor
	 */
	normal = normalize (divisor1);
	divisorc = divisor1 << normal;
	if (normal)
	    divisorc |= divisor2 >> (LBASE2 - normal);
	divisorc2 = divisor2 << normal;
	if (normal)
	    divisorc2 |= divisor3 >> (LBASE2 - normal);
#ifdef DEBUG
	printf ("divisor 1 %u divisor2 %u divisorc %u\n",
		(unsigned int) divisor1, (unsigned int) divisor2,
		(unsigned int) divisorc);
#endif
	/* 
	 * division just like humans, estimate each digit,
	 * correct by checking the partial product,
	 * then subtract the resultant product from
	 * the dividend (which has been copied into rem)
	 */
	for (index = remlen-1;
	     rem->length > 0 && index >= b->length;
	     index--)
	{
#ifdef DEBUG
	    prs ("rem:    ", rem);
	    prs ("b:      ", b);
#endif
	    /*
	     * estimate this digit
	     */
#ifdef DEBUG
	    printf ("digit1: %u digit2: %u divisorc: %u divisorc2: %u\n",
		    (unsigned int) NaturalDigits(rem)[index],
		    (unsigned int) NaturalDigits(rem)[index-1],
		    (unsigned int) divisorc,
		    (unsigned int) divisorc2);
#endif
	    /* 
	     * Using the scale factor and scaled divisor computed
	     * above, compute an estimate for this digit
	     *
	     * This is computing
	     *
	     *	    d1 = top dividend digit
	     *	    d2 = next dividend digit
	     *	    v1 = top divisor digit
	     *	    v2 = next divisor digit
	     *
	     *	      (d1 * base + d2)	<- original temp
	     * temp = ----------------
	     *        (v1 + v2 / base)	<- divisorc is scaled from this
	     *
	     * To make this computation work in integers,
	     * the top and bottom are scaled by the largest amount
	     * possible (base / (v1 + 1)); for powers of two, shift
	     * the divisor until the MSB is one.
	     */
	    dividend1 = NaturalDigits(rem)[index];
	    dividend2 = dividend3 = dividend4 = 0;
	    if (index)
	    {
		dividend2 = NaturalDigits(rem)[index - 1];
		if (index > 1)
		{
		    dividend3 = NaturalDigits(rem)[index - 2];
		    if (index > 2)
			dividend4 = NaturalDigits(rem)[index - 3];
		}
	    }
	    temp = (double_digit) dividend1 << LBASE2;
	    temp |= (double_digit) dividend2;
	    if (normal)
	    {
		temp = (temp << normal) | (dividend3 >> (LBASE2 - normal));
	    }
	    if ((digit) (temp >> LBASE2) == divisorc)
		d = (digit) 0xffffffffL;
	    else
		d = temp / divisorc;
#ifdef DEBUG	    
	    printf ("temp 0x%08x%08x divisorc 0x%08x d 0x%08x\n", 
		    (unsigned int) (temp >> 32), (unsigned int) temp,
		    divisorc, d);
#endif
	    offset = index - b->length;
	    if (d)
	    {
		temp2 = dividend3;
		if (normal)
		    temp2 = (temp2 << normal) | (dividend4 >> (LBASE2 - normal));
		temp = temp - (double_digit) d * (double_digit) divisorc;
		while (temp < BASE &&
		       (double_digit) divisorc2 * d > temp * BASE + temp2)
		{
		    --d;
		    temp += (double_digit) divisorc;
		}
		if (d == 1)
		{
		    carry = NaturalSubtractOffset (rem, b, offset);
		}
		else
		{
		    NaturalDigitMultiply (b, d, partial);
		    carry = NaturalSubtractOffset (rem, partial, offset);
		}
		if (carry)
		{
#ifdef DEBUG
		    printf ("add back 0x%x\n",
			    (unsigned int) d);
		    prs ("a", a);
		    prs ("b", b);
		    prs ("rem", rem);
		    prs ("partial", partial);
#endif
		    NaturalAddOffset (rem, b, offset);
#ifdef DEBUG
		    prs ("rem", rem);
#endif
		    d--;
		}
	    }
#ifdef DEBUG
	    prs ("partial", partial);
	    printf ("result[%d] = %u\n", (int) offset, 
		    (unsigned int) d);
#endif
	    NaturalDigits(quo)[offset] = d;
	}
	/*
	 * clean up the rest of the digits
	 */
	for (offset = index - b->length; offset >= 0; offset --) {
#ifdef DEBUG
	    printf ("result[%d] zeroed\n", offset);
#endif
	    NaturalDigits(quo)[offset] = 0;
	}
	for (index = quo->length - 1; index >= 0; index--)
	    if (NaturalDigits(quo)[index] != 0) {
		quo->length = index + 1;
		break;
	    }
#ifdef DEBUG
	prs ("quo: ", quo);
	prs ("rem: ", rem);
#endif
    }
#ifdef CHECK
    {
	Natural	*check;

	check = NaturalPlus (NaturalTimes (quo, b), rem);
	if (!NaturalEqual (check, a) ||
	    !NaturalLess (rem, b))
	{
	    prs ("a:     ", a);
	    prs ("b:     ", b);
	    prs ("quo:   ", quo);
	    prs ("rem:   ", rem);
	    prs ("check: ", check);
	    printf ("divide failed\n");
	}
    }
#endif
    EXIT ();
    /*
     * Stack guaratees that this will work -- it saves a place
     * for pushing at least as far as the previous StackReset, or
     * a stack chunk, whichever is smaller.
     */
    REFERENCE (quo);
    REFERENCE (rem);
    *remp = rem;
    return quo;
}