// ************************************************************************** //
//                                                                            //
//    eses                   eses                                             //
//   eses                     eses                                            //
//  eses    eseses  esesese    eses   Embedded Systems Group                  //
//  ese    ese  ese ese         ese                                           //
//  ese    eseseses eseseses    ese   Department of Computer Science          //
//  eses   eses          ese   eses                                           //
//   eses   eseses  eseseses  eses    University of Kaiserslautern            //
//    eses                   eses                                             //
//                                                                            //
// ************************************************************************** //
//          ** Note: We are only interested in even values of B **            // 
// -------------------------------------------------------------------------- //
// The following module implements the basic digit recurrence algorithm for   //
// division of integers: Given B-complement numbers x and y, the algorithm    //
// computes the quotient digits of the quotient x/y to radix B as well as the //
// remainder r. To this end, the algorithm computes a sequence of numbers     //
// z(m),...,z(0):                                                             //
//  (1) z(m) = x - alpha(q[M]) * |y[N-1..0]| * B^M                            //
//  (2) z(i-1) = z(i) - q[i-1] * |y[N-1..0]| * B^(i-1) where q[i-1] is the    //
//      maximal digit ∈ {0,...,B−1} with 0 ≤ z(i-1)                           //
// It can be easily proved that the following invariants are then true:       //
//  (1) x = z(i) + intval(q[M−1,…,i])*|y|*B^i                                 //
//  (2) 0 ≤ z(i) < |y|*B^i                                                    //
// Hence, for i=0, it follows that 0<= z(0) < |y| and x = intval(q)*|y|+z(0), //
// so that z(0) is the remainder and the q[i] are the digits of the quotient  //
// of the division of x by y.                                                 //
//   The numbers z(i) are stored in the same variable z that is overwritten   //
// step by step. The quotient digit q[i-1] is determined naively by simply    //
// calculating z(i) - q[i-1] * |y[N-1..0]| * B^(i-1) for each digit q[i-1]    //
// starting from B-1 until the result is non-negative (note to this end, that //
// the numbers z(i) are always nonnegative).                                  //
// ************************************************************************** //

macro B = 4;    // base of the radix numbers
macro M = 4;    // number of digits used for x
macro N = 3;    // number of digits used for y
macro alpha(x) = (x<B/2 ? +x : +x-B);
macro gamma(y) = (y<0 ? y+B : y);
macro dval(x,i,k) = (i==k-1 ? alpha(x[i]) : +x[i]);
macro natval(x,m) = sum(i=0..m-1) (x[i] * exp(B,i));
macro intval(x,k) = sum(i=0..k-1) (dval(x,i,k) * exp(B,i));


module IntDivSeq0([M]nat{B} ?x,[N]nat{B} ?y,[M+1]nat{B} q,[N]nat{B} r,event !rdy) {
    [M+N]nat{B} z;    // partial remainders
    [N]nat{B} s;      // auxiliary variable for sums
    [N+1]int{B} c;    // auxiliary variable for carry digits
    int{2} sgn_y;     // sign of y
    event isNegative; // is set when subtraction yields negative result

    //-------------------------------------------------------------------------
    // copy x to z using digit extension
    //-------------------------------------------------------------------------
    for(j=0..M+N-1)
        z[j] = (j<M ? x[j] : (x[M-1]<B/2 ? 0 : B-1));

    //-------------------------------------------------------------------------
    // determine leftmost digit of the quotient and the first partial remainder
    //-------------------------------------------------------------------------
    sgn_y = (y[N-1]<B/2 ? +1 :  -1);
    q[M]  = (x[M-1]<B/2 ?  0 : B-1);
    // compute z[N+M-1..M] = z[N+M-1..M] - alpha(q[M]) * sgn_y * y
    c[0] = 0;
    for(j=0..N-1)
        // note that -(B-1)<=sm<=2*B-1, thus -1<=sm/B<=1
        let(yJ = (j==N-1 ? alpha(y[j]) : y[j]))
        let(sm = z[M+j] + c[j] - alpha(q[M]) * sgn_y * yJ) {
        next(z[M+j]) = sm % B;
        c[j+1] = sm / B;
        }
    // assignment to z[M+N] is not required since z[M+N]==0 would hold by def.
    // thus, z rather represents a radix-B number (we omit the leftmost digit)
    w0: pause;

    //-------------------------------------------------------------------------
    // computing the remaining quotient digits
    //-------------------------------------------------------------------------
    for(iup=0..M-1)
        let(i=M-1-iup) {// thus i=M-1..0
        // try digits q[i] = B-1,B-2,... until we reach one where 
        // z[N+i..i] - q[i] * sgn_y * y[N-1..0] remains non-negative
        next(q[i]) = B-1;
        do {
            w1: pause;
            // compute s[N..0] = z[N+i..i] - q[i] * sgn_y * y[N-1..0]
            c[0] = 0;
            for(j=0..N-1)
                // note that -B*(B-1)<=sm<=B-1, thus -(B-1)<=sm/B<=0
                let(yJ = (j==N-1 ? alpha(y[j]) : y[j]))
                let(sm = z[i+j] + c[j] - q[i] * sgn_y * yJ) {
                s[j]   = sm % B;
                c[j+1] = sm / B;
                }
            isNegative = (z[N+i] + c[N] < 0);
            if(isNegative) 
                next(q[i]) = q[i] - 1;
        } while(isNegative);
        // Here, q[i] is the largest digit with non-negative s[N..0].
        // We therefore commit the subtraction in that we copy the preliminary
        // result from s to z before computing the next digit q[i].
        w2: pause;
        z[N+i] = 0;
        for(j=0..N-1)
            z[j+i] = s[j];
    }
    //-------------------------------------------------------------------------
    // copy the remainder from z
    //-------------------------------------------------------------------------
    for(i=0..N-1)
        r[i] = z[i];
    //-------------------------------------------------------------------------
    // here we have x = q * |y| + r; so if y<0 holds, we have to negate q 
    //-------------------------------------------------------------------------
    if(y[N-1]>=B/2) {
        [M+1]nat{2} cq;
        cq[0]=0;
        for(i=0..M-1) 
            let(sm = -(cq[i] + q[i])) {
            cq[i+1] = -(sm / B);
            next(q[i]) = sm % B;
        }
        let(sm = -(cq[M] + alpha(q[M])))
        next(q[M]) = gamma(sm % B);
        w3:pause;
    }
    emit(rdy);
    //-------------------------------------------------------------------------
    // final assertion
    //-------------------------------------------------------------------------
    if(intval(y,N)!=0) {
        assert(intval(r,N) == intval(z,M+N));
        assert(intval(x,M) == intval(q,M+1) * intval(y,N) + intval(r,N));
        assert(intval(r,N) < abs(intval(y,N)));
    }
}