//  ============================================================
//
//  divide.cpp
//
//  Copyright 2002, Dennis Meilicke and Rene Peralta
//
//  ============================================================
//
//  Description:
//
//  Implementation functions for division.
//
//  ============================================================
#include <lnv3/ln3impl.h>


//  ============================================================


void
divide( const Base *U, digit_t lv, Base *Q, digit_t *r )
//
//  Scalar division.
//
//  Calculate Q = floor( U / lv ), remainder r
//      where lv and r are scalars.
//
{
    if(lv == 0) { int x=1; x /= (x - 1); }

    if( U->size == 1 )
    {
        //
        //  There could be a sign problem, here
        //
        Q->digit[ 0 ] = U->digit[ 0 ] / lv;
        Q->size = ( Q->digit[0] == 0 ) ? 0 : 1;
        *r = U->digit[ 0 ] % lv;
        return;
    }

    ddSplit_t llr; llr.all = 0;

    short j;
    for( j=U->size - 1 ; j>=0 ; j-- )
    {
        llr.dd.high = llr.dd.low;
        llr.dd.low = U->digit[ j ];
        Q->digit[ j ] = (digit_t)( llr.all / lv );
        llr.all = llr.all % lv;
    }

    Q->size = U->size;
    while( (Q->size>=0) && (Q->digit[Q->size - 1]==0) )
        Q->size--;

    if( ( Q->size == 1 ) && ( Q->digit[0] == 0 ) )
        Q->size = 0;

    *r = llr.dd.low;

    return;
}


//  ============================================================


static
digit_t
getQHat( Base *u, Base *v, size_t j )
//
//  Helper function for "large" divide.
//
//  Estimate a digit of the quotient.
//
//  There is too much subscripting in this function.  All
//  of the subscripting should be changed to pointer math.
//  Also, some of the values are used often enough that 
//  they should be cached.
//
{
    ddSplit_t firstTwo;
    firstTwo.dd.high = u->digit[ j ];
    firstTwo.dd.low  = u->digit[ j - 1 ];

    digit_t qHat;
    if( u->digit[ j ] == v->digit[ v->size - 1 ] )
        qHat = 0xffffffff;
    else
    {
        // ddigit_t denom = 0;
        // lowDigit( denom ) = v->digit[ v->size - 1 ] ;
        // ddigit_t quo;
        // quo = firstTwo / denom;
        // qHat = lowDigit( quo );
        qHat = (digit_t)( firstTwo.all / (ddigit_t)v->digit[ v->size - 1 ] );
    }

    ddigit_t  v1qHat = (ddigit_t)qHat * (ddigit_t)v->digit[ v->size - 1 ];
    ddigit_t  v2qHat = (ddigit_t)qHat * (ddigit_t)v->digit[ v->size - 2 ];

    ddSplit_t rHat;
    rHat.all = firstTwo.all - v1qHat;

    ddSplit_t testVal;
    while( qHat )
    {
        if( rHat.dd.high )
            break;

        testVal.dd.high = rHat.dd.low;
        testVal.dd.low  = u->digit[ j - 2 ];
        if( v2qHat <= testVal.all )
            break;

        qHat--;
        v2qHat   -= (ddigit_t)v->digit[ v->size - 2 ];
        rHat.all += (ddigit_t)v->digit[ v->size - 1 ];
    }

    return qHat;
}


//  ============================================================


static
bool
multAndSub( Base *u, Base *v, size_t j, digit_t qHat )
//
//  Compute:
//      x = qHat * v
//
//  then:
//      u_j u_j-1 u_j-2  ...  u_j-n
//    - x_n x_n-1 x_n-2  ...  x_0
//      ---------------------------
//      u_j u_j-1 u_j-2  ...  u_j-n
//
//  Note that the two calculations are interleaved.
//
//  Return:  true  if qHat was too big.  That is, if the
//                 would have resulted in a negative number.
//           false otherwise.
//
{
    if( qHat == 0 ) return false;

    digit_t *a = &u->digit[ j - v->size ];
    digit_t *b =  v->digit;

    ddSplit_t temp;  temp.all  = 0;
    ddSplit_t x;     x.all     = 0;
    ddSplit_t prevX; prevX.all = 0;
    unsigned short borrow = 0, prevBorrow = 0;

    for( size_t i=0 ; i<v->size ; i++, a++, b++ )
    {
        borrow = 0;

        x.all = (ddigit_t)*b * (ddigit_t)qHat;
        
        temp.dd.high = 0;
        temp.dd.low  = *a;
        if( temp.dd.low < x.dd.low )
        {
            temp.dd.high = 1;
            borrow++;
        }
        temp.all -= (ddigit_t)(x.dd.low);

        if( temp.dd.low < prevX.dd.high )
        {
            temp.dd.high = 1;
            borrow++;
        }
        temp.all -= (ddigit_t)(prevX.dd.high);

        if( temp.dd.low < prevBorrow )
        {
            temp.dd.high = 1;
            borrow++;
        }
        temp.all -= (ddigit_t)prevBorrow;

        *a = temp.dd.low;

        prevX = x;
        prevBorrow = borrow;
    }

    borrow = 0;
    temp.dd.high = 0;
    temp.dd.low = *a;
    if( temp.dd.low < prevX.dd.high )
    {
        temp.dd.high = 1;
        borrow++;
    }
    temp.all -= (ddigit_t)(prevX.dd.high);

    if( temp.dd.low < prevBorrow )
    {
        temp.dd.high = 1;
        borrow++;
    }
    temp.all -= (ddigit_t)prevBorrow;

    *a = temp.dd.low;

    return ( borrow != 0 );
}


//  ============================================================


static
void
addBack( Base *u, Base *v, size_t j )
//
//  Add v back in, ignoring the final carry:
//
//     u_j u_j-1 u_j-2  ...  u_j-n
//  +  0   v_n-1 v_n-2  ...  v_0
//
{
    digit_t *a = &(u->digit[ j - v->size ]);
    digit_t *b = v->digit;

    ddSplit_t accum; accum.all = 0;
    for( size_t i=0 ; i<v->size ; i++, a++, b++ )
    {
        //  highDigit( accum ) is the carry
        accum.all = (ddigit_t)*a + (ddigit_t)*b + accum.dd.high;
        *a = accum.dd.low;
    }
    *a += accum.dd.high;
}


//  ============================================================


void
divide( const Base *N, const Base *D, Base *Q, Base *R )
//
//  Compute N / D = ( Q, R ), where
//
//      N = D * Q + R
//
//  This function assumes that the input is proper.  That is:
//
//  1.  D != 0
//  2.  N > D
//  3.  D->size > 1
//
{
    Base *u = R; *u = *N;
    u->count = 1; // fix leak?
    Base *v = new Base( *D );

    size_t sizeU = u->size;

    size_t d = 0;
    d = digitSizeInBits - NumBits( v->digit[ v->size - 1] );

    if( d )
    {
        ShiftLeftEqual( u, d );
        ShiftLeftEqual( v, d );
    }

    if( u->size == sizeU )
    {
        u->digit[ u->size ] = 0;
        u->size++;
    }

    Q->size = u->size - v->size;
    digit_t *q = Q->digit + Q->size - 1;

    for( size_t j=(u->size - 1); j>=v->size ; j--, q-- )
    {
        digit_t qHat;
        qHat = getQHat( u, v, j );

        bool qHatTooBig;
        qHatTooBig = multAndSub( u, v, j, qHat );

        if( qHatTooBig )
        {
            addBack( u, v, j );
            qHat--;
        }

        *q = qHat;
    }

    if( d )
        ShiftRightEqual( u, d );

    //  Check for leading zeros in R = u
    while( u->size && ( u->digit[ u->size - 1 ] == 0 ) )
        u->size--;

    //  Check for leading zeros in Q
    while( Q->size && ( Q->digit[ Q->size - 1 ] == 0 ) )
        Q->size--;

    if( N->sign == D->sign )
        Q->sign = positive;
    else
        Q->sign = negative;

    delete v;
}


//  ============================================================


void
Divide( const Base *N, const Base *D, Base *Q, Base *R )
{
    if( D->size == 0 ) { int x=1; x/=(x-1); } // raise divide-by-zero exception
    if( N->size == 0 )
    {
        Q->size = 0;
        R->size = 0;
        return;
    }
    if( D->size == 1 )
    { 
        digit_t r;
        divide( N, D->digit[0], Q, &r );
        if( r == 0 )
        {
            R->sign = positive;
            R->size = 0;
        }
        else
        {
            R->sign = N->sign;
            R->size = 1;
            R->digit[ 0 ] = r;
        }
        ( N->sign == D->sign ) ? Q->sign = positive : Q->sign = negative;
        return;
    }

    compare_t comp = AbsCompare( N, D );
    if( comp == lessThan ) { Q->size = 0; *R = *N; R->count = 1; return ;}
    if( comp == equalTo  ) { Q->size = 1; Q->digit[0] = 1; R->size = 0; return; }

    divide( N, D, Q, R );
}


//  ============================================================


void
Divide( const Base *N, long d, Base *Q, long *r )
{
    if( d == 0 ) { int x=1; x/=(x-1); }  //  raise divide-by-zero exception
    if( N->size == 0 )
    {
        Q->size = 0;
        *r = 0;
        return;
    }

    divide( N, d, Q, (digit_t *)r );

    if( ( r < 0 ) && ( N->sign == positive ) )
        *r = -(*r);
    else if( ( r > 0 ) && ( N->sign == negative ) )
        *r = -(*r);

    ( N->sign == sign( d ) ) ? Q->sign = positive : Q->sign = negative;
    return;
}


//  ============================================================


void
Divide( const Base *N, unsigned long d, Base *Q, unsigned long *r )
//
//  Note that this function doesn't really make sense in the 
//  case where N->sign == negative.  This would imply that the
//  sign of r be negative, but since r is unsigned, this is 
//  impossible.
//
//  Take great care when calling this function.
//
{
    if( d == 0 ) { int x=1; x/=(x-1); }  //  raise divide-by-zero exception

    divide( N, d, Q, r );

    Q->sign = N->sign;
    return;
}