// ************************************************************************** //
//                                                                            //
//    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                                             //
//                                                                            //
// ************************************************************************** //
// Given numbers a[0..N-1] and b[0..N-1] as well as a further initial value   //
// y0, the task is to compute the numbers y[0..N-1] that are recursively      //
// defined as y[0] := a[0]*y0+b[0], and y[i+1] := a[i+1]*y[i]+b[i+1].         //
// A clever idea of [OkNT95] is to reduce the problem to the parallel prefix  //
// computation as follows: We define vectors x[i] = [y[i],1] and matrices     //
// M[i] := [[a[i],b[i]],[0,1]] so that x[i+1] = M[i] * x[i] follows. By       //
// induction it also follows that x[i+1] = M[i]*...*M[0]*x[0] holds, so       //
// that we can compute all x[i+1] from the matrices M[i]*...*M[0] and the     //
// initial value x[0]. Moreover, it is easy to prove that the following holds //
// for 2x2 matrices:                                                          //
//  | a b |   | c d |   | a*c a*d+b |                                         //
//  | 0 1 | * | 0 1 | = |  0    1   |                                         //
// Thus, we can ignore the second rows, since these remain constant.          //
// The algorithm below computes therefore from the given values a[0..N-1] and //
// b[0..N-1] the first rows of the matrices M[i]*...*M[0] for i=0,...,N-1.    //
// Since matrix multiplication is associative, we can apply parallel prefix   //
// computation to obtain these matrices in time O(log(N)) using O(N) procs.   //
// and O(N) work.                                                             //
// ************************************************************************** //

macro N = 8;


module FirstOrderLinearRec(int ?y0,[N]int ?a,?b,!y,event !rdy){
    [N][2]int m;

    for(i=0..N-1) {
        m[i][0] = a[i];
        m[i][1] = b[i];
    }
    for(i=0..log(N)-1) {
        for(j=exp(2,i)..N-1) {
            next(m[j][0]) = m[j][0] * m[j-exp(2,i)][0];
            next(m[j][1]) = m[j][0] * m[j-exp(2,i)][1] + m[j][1];
        }
        pause;
    }
    // having computed the matrix products, we compute the y[i]
    for(i=0..N-1) {
        y[i] = m[i][0] * y0 + m[i][1];
    }
    emit(rdy);
}
drivenby {
    a = [+1,+2,-2,-1,+3,-3,+0,+1];
    b = [+3,+1,+0,+1,-2,-1,+3,-4];
    y0 = +1;
    await(rdy);
    for(i=0..N-1)
        assert(y[i] == a[i] * (i==0 ? y0 : y[i-1]) + b[i]);
}