// ************************************************************************** //
//                                                                            //
//    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                                             //
//                                                                            //
// ************************************************************************** //
// This module implements the restoring division algorithm as a further       //
// improvement of NatDivModSeq1. The idea is to subtract the given dividend   //
// as long as the resulting difference is non-negative. Once the result became//
// negative, a restore operation takes place, where the given dividend is     //
// added. The advantage over NatDivModSeq1 is that only a simple subtraction  //
// is required, i.e., one without a previous multiplication with a digit q[i].//
//   As all other digit selection algorithms for division, also this version  //
// is based on the computation of the following 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(i).     //
// 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.                                                        //
//   As in NatDivModSeq1, the numbers z(i) are stored in r[N-1..0]@x[i-1..0]. //
// The quotient digit q[i-1] is determined naively by simply subtracting      //
// y[N-1..0] * B^(i-1) from z(i) until the result becomes negative. The number//
// of subtractions is then the digit q[i-1]. After the result became negative,//
// an addition z(i) + y[N-1..0] * B^(i-1) is performed to restore the correct //
// partial remainder z(i-1).                                                  //
//   The sign of the result of the subtraction is computed as follows: Note   //
// that z(i), i.e., s[N..0] has one more digit than y[N-1..0]. Thus, the      //
// leftmost digit of the sum is s[N] = (-c[N] + s[N]) % B and the next carry  //
// would be c[N+1] = abs((-c[N] + s[N]) / B). As usual for subtraction, the   //
// result is negative iff the final carry is not zero, so our result is neg.  //
// iff (-c[N] + s[N]) / B !=0. Since c[N] is either 0 or 1, it follows that   //
// this is the case iff c[N]==1 & s[N]==0 holds. This allows us to determine  //
// the sign even without computing c[N+1].                                    //
// ************************************************************************** //

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 NatDivRestore([M]nat{B} ?x,[N]nat{B} ?y,
                        [M]nat{B}  q,[N]nat{B}  r,event !rdy) {
    [N+1]nat{B} s;    // sum digits of difference
    [N+1]nat{2} c;    // carry digits of subtraction
    event isNegative; // is set when subtraction yields negative result
    // -------------------------------------------------------------------------
    // prepare first partial remainder z(m) = r[N-1..0]@x[M-1]
    // -------------------------------------------------------------------------

   for(i=0..N-1)
        r[i] = 0;

    for(iup=0..M-1)
        let(i=M-1-iup) {// thus i=M-1..0
        // ---------------------------------------------------------------------
        // r[N-1..0]@x[i] is the initial partial remainder z(i+1)
        // ---------------------------------------------------------------------
        s[0] = x[i];
        for(j=1..N)
            s[j] = r[j-1];
        // ---------------------------------------------------------------------
        // now subtract y[N-1..0] from z(i+1) until the result becomes negative
        // ---------------------------------------------------------------------
        next(q[i]) = 1;
        do {
            // compute next(s[N..0]) = s[N..0] - y[N-1..0]
            c[0] = 0;
            for(j=0..N-1)
                let(sm = +s[j] - (y[j] + c[j])) {
                next(s[j]) = sm % B;
                c[j+1] = -(sm / B);
                }
            next(s[N]) = (+s[N] - c[N]) % B;
            next(isNegative) = (s[N]<c[N]);
            sub: pause;
            if(!isNegative & q[i]!=B-1) 
                next(q[i]) = q[i]+1;
        } while(!isNegative & q[i]!=B-1);
        // ---------------------------------------------------------------------
        // restore if s became negative; and in each case, put the 
        // result in r, so that the next partial remainder is obtained
        // ---------------------------------------------------------------------
        if(isNegative) { 
            // restore step next(r[N-1..0]) = y[N-1..0] + s[N-1..0]
            c[0] = 0;
            for(j=0..N-1)
                // note that 0<=sm<=2*B-1, thus 0<=sm/B<=1
                let(sm = c[j] + s[j] + y[j]) {
                next(r[j]) = sm % B;
                c[j+1] = sm / B;
                }
            next(q[i]) = q[i]-1;
            add: pause;
        } else {
            // shift the resulting sum into r
            for(j=0..N-1)
                next(r[j]) = s[j];
            sft: pause;
        }
    }
    emit(rdy);
    // -------------------------------------------------------------------------
    // check assertions
    // -------------------------------------------------------------------------
    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,0,1];
    y = [2,1,3];
    await(rdy);
}