// ************************************************************************** //
//                                                                            //
//    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                                             //
//                                                                            //
// ************************************************************************** //
// The following module implements a very simple algorithm for division of    //
// radix-B numbers: Given radix-B numbers x and y with M and N digits, the    //
// algorithm computes a sequence of numbers z(m),...,z(0):                    //
//  (1) z(m) = x                                                              //
//  (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 q[i−1]*y[N-1..0]*B^(i-1) ≤ z         //
// It can be easily proved that the following invariants are then true:       //
//  (1) x = z(i) + q[M−1,...,i]*y*B^i                                         //
//  (2) z(i) < y*B^i                                                          //
// Hence, for i=0, it follows that z(0) < y and x = q*y+z(0), so the the final//
// 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   //
// with these numbers 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.           //
//   For the definition of isNegative, note that the leading digit of s would //
// be c[N+1] = (-c[N] + z[i+N]) / B, which is either -1 or 0. Thus, we have   //
// isNegative <-> c[N+1]==-1 <-> -c[N] + z[i+N] < 0 <-> z[i+N] < c[N].        //
// ************************************************************************** //

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 natval(x,m) = sum(i=0..m-1) (x[i] * exp(B,i));


module NatDivSeq0([M]nat{B} ?x,[N]nat{B} ?y,
                     [M]nat{B}  q,[N]nat{B}  r,event !rdy) {
    [M+N]nat{B} z;    // preliminary dividend
    [N+1]nat{B} s,c;  // auxiliary variables
    event isNegative; // is set when subtraction yields negative result
    // -------------------------------------------------------------------------
    // copy x to right end of z
    // -------------------------------------------------------------------------
    for(j=0..M+N-1)
        z[j] = (j<M ? x[j] : 0);
    //-------------------------------------------------------------------------
    // computing the quotient digits
    //-------------------------------------------------------------------------
    for(iup=0..M-1)
        let(i=M-1-iup) {// thus i=M-1..0
        // we try digits q[i] = B-1,B-2,... until we reach one where 
        // z[N+i..i] - q[i] * y[N-1..0] remains non-negative
        next(q[i]) = B-1;
        do {
            w0: pause;
            //------------------------------------------------------------------
            // compute s[N..0] = z[N+i..i] - q[i] * y[N-1..0]
            //------------------------------------------------------------------
            c[0] = 0;
            for(j=0..N-1)
                // note that -B*(B-1)<=sm<=B-1, thus sm has type int{B}
                // and -cout has type nat{B} due to -(B-1)<=sm/B<=0
                let(sm = +z[i+j] - q[i] * y[j] + c[j]) {
                s[j]   = sm % B;
                c[j+1] = abs(sm / B);
                }
            s[N] = abs((-c[N] + z[i+N]) % B);
            isNegative = z[i+N]<c[N]; // determine sign instead of c[N+1]
            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].
        //----------------------------------------------------------------------
        w1: pause;
        for(j=0..N)
            z[i+j] = s[j];
    }
    //-------------------------------------------------------------------------
    // copy the remainder from z
    //-------------------------------------------------------------------------
    for(i=0..N-1)
        r[i] = z[i];
    emit(rdy);
    //-------------------------------------------------------------------------
    // check the assertions
    //-------------------------------------------------------------------------
    if(exists(i=0..N-1) (y[i]!=0)) {
        assert(natval(x,M) == natval(y,N) * natval(q,M) + natval(r,N));
        assert(natval(r,N) < natval(y,N));
    }
}
drivenby Test01 {
    x = [1,3,2,1];
    y = [2,1,3];
    await(rdy);
}