// ************************************************************************** //
//                                                                            //
//    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 **            // 
// -------------------------------------------------------------------------- //
// This module implements restoring B-complement division in analogy to       //
// the restoring radix-B division given in NatDivRestore.qrz.                 //
// ************************************************************************** //

macro B = 4;    // base of the radix numbers
macro M = 3;    // number of digits used for x
macro N = 2;    // 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));
macro z(i,j) = (j==0 ? x[i] : r[j-1]); // should be z[j+i] of IntDivModSeq0

module IntDivRestore([M]nat{B} ?x,[N]nat{B} ?y,[M+1]nat{B} q,[N]nat{B} r,event !rdy) {
    [N+1]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

    //-------------------------------------------------------------------------
    // apply digit extension to x and put the extended digits in r
    //-------------------------------------------------------------------------
    for(j=0..N-1)
        r[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 r[N-1..0] = r[N-1..0] - 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 = sgn_y * (j==N-1 ? alpha(y[j]) : y[j]))
        let(sm = r[j] + c[j] - alpha(q[M]) * yj) {
        next(r[j]) = sm % B;
        c[j+1] = sm / B;
        }
    w0: pause;
    //-------------------------------------------------------------------------
    // computing the remaining quotient digits
    //-------------------------------------------------------------------------
    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(yj = sgn_y * (j==N-1 ? alpha(y[j]) : y[j]))
                let(sm = +s[j] - (yj + 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);
        
        // determine next partial remainder by possible restore step
        if(isNegative) { 
            // restore by adding next(r[N-1..0]) = y[N-1..0] + s[N-1..0]
            c[0] = 0;
            for(j=0..N-1)
                let(yj = sgn_y * (j==N-1 ? alpha(y[j]) : y[j]))
                let(sm = c[j] + s[j] + yj) {
                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;
        }
    }
    //-------------------------------------------------------------------------
    // 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(x,M) == intval(q,M+1) * intval(y,N) + intval(r,N));
        assert(intval(r,N) < abs(intval(y,N)));
    }
}