// ************************************************************************** // // // // 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]); }