// ************************************************************************** //
//                                                                            //
//    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 division of two n-bit radix-2 numbers [N]bool x,y can be reduced to    //
// the problem to compute the reciprocal z of y with n-bits, i.e., precision  //
// 2^{-n} which is a number z with 0 <= 1/y - z <= 2^{-n}. Given z, we can    //
// compute the quotient either as q=floor(x*z) or q=floor(x*z)+1, depending on//
// whether the remainder x-q*y satisfies x-q*y<y.                             //
//  To compute an approximation of the reciprocal, we apply Newton's iteration//
// to the function f(x) := 1/x - b so that x_{i+1} := x_i*(2-y*x_i) is        //
// obtained. For the computation, we first scale the operand y by shifting it //
// to the right so that the first bit on the right of the binary point is 1.  //
// This way, the new operand ys is in the interval 0.5 <= ys < 1.0 so that its//
// reciprocal is 2.0 >= 1/ys > 1.0. It is then sufficient to work with fixed  //
// point numbers with two bits on the left and N bits to the right of the     //
// binary point. The result z contains finally the N computed bits of the     //
// approximation to the reciprocal of 1/y, where only the N bits to the right //
// of the binary point are given.                                             //
// ************************************************************************** //

macro N = 8;
macro bval(x) = (x ? 1 : 0);
macro fval(x) = exp(2.0,-N) * sum(i=0..N+1) (bval(x{i}) * exp(2,i));
macro fmul(x,y) = (nat2bv(bv2nat(x) * bv2nat(y),2*(N+2)){2*N+1:N});
macro fsub(x,y) =  nat2bv(bv2nat(x) - bv2nat(y),N+2);
macro rsh(x) = false@(x{-1:1});
macro lsh(x) = x{-2:0}@false;

macro one = false @ true  @ {false::N};
macro two = true  @ false @ {false::N};


module ReciprocalNewton(bv{N} ?y,bv{N} z,event rdy) {
    bv{N+2} x_old,x_new,ys;
    nat{N+1} j;
    real r_new,r_ys;

    // scale the operand y to ys = y*2^-j into interval [0.5,1.0]
    // determine j such that y{j-1} is the leftmost bit that is true
    ys = false@false@y;
    j = N;
    while(!ys{-3}) {
        r_ys = fval(ys);
        next(ys) = lsh(ys);
        next(j) = j-1;
        wu: pause;
    }
    r_ys = fval(ys);

    // iteration to compute 1/ys starting from 1.5 with max. initial error 0.5
    // log(N) iterations are required to compute the N bits
    x_new = false@true@true@{false::N-1};
    
    do {
        r_new = fval(x_new);
        next(x_old) = x_new; 
        next(x_new) = fmul(x_new,fsub(two,fmul(ys,x_new)));
        wi: pause;
    } while(x_old != x_new);
    r_new = fval(x_new);

    // having 1/ys, scale it to obtain 1/y
    while(j>0) {
        next(x_new) = rsh(x_new);
        next(j) = j-1;
        wr: pause;
        r_new = fval(x_new);
    }

    r_new = fval(x_new);
    z = x_new{N-1:0};
    emit(rdy);

}
drivenby {
    y = nat2bv(128,N);
    await(rdy);
}