// ************************************************************************** //
//                                                                            //
//    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 algorithm below solves a linear equation system given as a MxN matrix  //
// a in that it applies the Gauss Jordan algorithm.                           //
// ************************************************************************** //

macro M = 3;
macro N = 4;

// define arithmetic over the rationals
macro add(a,b) = (a.0*b.1+a.1*b.0, a.1*b.1);    // addition of rationals
macro sub(a,b) = (a.0*b.1-a.1*b.0, a.1*b.1);    // subtraction of rationals
macro mul(a,b) = (a.0*b.0,a.1*b.1);             // multiplication of rationals
macro div(a,b) = (a.0*b.1,a.1*b.0);             // division of rationals


module LinearEquGaussJordanRat([M][N](int * int) ?a,b,event !rdy) {

    // store input matrix a in matrix b
    for(i=0..M-1)
        for(j=0..N-1)
            b[i][j] = a[i][j];

    // perform Gaussian-Jordan algorithm using arithmetic of rationals
    for(k=0..M-1) {
        for(i=0..M-1) {
            let(lambda = div(b[i][k],b[k][k])) {
            for(j=0..N-1)
                if(i==k)
                    next(b[i][j]) = div(b[i][j],b[k][k]);
                else
                    next(b[i][j]) = sub(b[i][j],mul(lambda,b[k][j]));
            }
        }
        pause;
    }
    // cancel by gcd
    for(i=0..M-1)
        for(j=0..N-1) {
            int x,y;
            x = b[i][j].0;
            y = b[i][j].1;
            while(y != 0) {
                next(x) = y;
                next(y) = x % y;
                pause;
            }
            next(b[i][j].0) = b[i][j].0 / x;
            next(b[i][j].1) = b[i][j].1 / x;
        }
    pause;
    emit(rdy);
}
drivenby {
    a[0][0].0 = 3;  a[0][0].1 = 1;
    a[0][1].0 = 2;  a[0][1].1 = 1;
    a[0][2].0 = 1;  a[0][2].1 = 1;
    a[0][3].0 = 0;  a[0][3].1 = 1;
    a[1][0].0 = 4;  a[1][0].1 = 1;
    a[1][1].0 = 3;  a[1][1].1 = 1;
    a[1][2].0 = 1;  a[1][2].1 = 1;
    a[1][3].0 = 1;  a[1][3].1 = 1;
    a[2][0].0 = 9;  a[2][0].1 = 1;
    a[2][1].0 = 3;  a[2][1].1 = 1;
    a[2][2].0 = 1;  a[2][2].1 = 1;
    a[2][3].0 = 3;  a[2][3].1 = 1;
    await(rdy);
}