// ************************************************************************** //
//                                                                            //
//    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 module below computes all roots of a given polynomial. The algorithm   //
// applies Newton-Raphson iteration to find an approximation root[i] of a root//
// and divides then the polynomial by the monomial (x-root[i]). To this end,  //
// the algorithm works with complex numbers to that we know that there are N  //
// complex roots. As soon as all roots have been computed with precision eps, //
// we polish them by another round of Newton-Raphson iterations starting with //
// the computed roots, but this time with the original polynomial.            //
//  The algorithm below is not recommended for use in scientific computing;   //
// to that purpose, it is better to use methods that reduce the problem to the//
// computation of the eigenvalues of the companion matrix of the polynomial.  //
// ************************************************************************** //


macro N = 5;

// macros on real numbers
macro rroot(x) = exp(x,0.5);
macro rsign(x) = (x<0 ? -1 : +1);

// arithmetic on complex numbers
macro cadd(x,y) = (x.0+y.0,x.1+y.1);
macro csub(x,y) = (x.0-y.0,x.1-y.1);
macro cmul(x,y) = (x.0*y.0-x.1*y.1,x.1*y.0+x.0*y.1);
macro cdiv(x,y) = ((x.0*y.0+x.1*y.1)/(y.0*y.0+y.1*y.1),
                   (x.1*y.0-x.0*y.1)/(y.0*y.0+y.1*y.1));
macro sqabs(x) = (x.1*x.1 + x.0*x.0);
macro cabs(x)  = rroot(x.1*x.1 + x.0*x.0);
macro croot(x) = (rsign(x.1) * rroot(0.5 * (cabs(x) + x.0)),
                               rroot(0.5 * (cabs(x) - x.0)));
macro Re(x) = x.0;
macro Im(x) = x.1;

// numeric precision
macro eps = 1.0e-10;


module PolynomialRootsNewton([N+1]real ?c,[N](real*real) root,event !rdy) {
    [N+1](real*real) cc,dd; // coefficients as complex numbers
    (real*real) pp,dp;      // value of p and its derivative
    (real*real) x;          // value used in Newton iteration

    // initialization: convert given coefficients to complex numbers
    for(i=0..N) {
        cc[i].0 = c[i];
        cc[i].1 = 0.0;
    }
   
    for(i=0..N-1) {
        // apply Newton-Raphson iteration to find a complex root
        do {
            // evaluate p(x) and its derivative at x in O(N) time, which can be
            // improved to log(N) time but does not matter for small polynomials
            next(pp) = cc[N];
            next(dp) = (0.0,0.0);
            pause;
            for(j=0..N-1) {
                next(dp) = cadd(cmul(dp,root[i]),pp);
                next(pp) = cadd(cmul(pp,root[i]),cc[N-1-j]);
                pause;
            }
            // to avoid division by zero, we cheat a bit
            let(dp1 = (dp==(0.0,0.0) ? (eps,eps) : dp)) 
            next(root[i]) = csub(root[i],cdiv(pp,dp1));
            next(x) = root[i];
            pause;
        } while(sqabs(csub(x,root[i])) > eps);        
        
        // deflation step: divide polynomial cc by (x-root[i]) to continue with
        // the computation of the remaining roots using a polynomial of lower degree
        dd[N] = (0.0,0.0);
        for(j=1..N)
            dd[j-1] = (j==N-i ? cc[N-i] 
                              : cadd(cc[j],cmul(root[i],dd[j])));
        pause;
        // now assign cc := dd to proceed with the computation of remaining roots
        for(j=0..N)
            cc[j] = dd[j];
    }

    pause;
    // At this point, we have computed all the N roots. However, due to the deflation
    // the precision of the computed roots might not be that good. For this reason,
    // we polish them by another round of Newton-Raphson iterations starting with the
    // so far obtained roots.

    for(i=0..N) {
        cc[i].0 = c[i];
        cc[i].1 = 0.0;
    }

    for(i=0..N-1) {
        // apply Newton-Raphson iteration starting with the current approximations
        do {
            // evaluate p(x) and its derivative at x in O(N) time, which can be
            // improved to log(N) time but does not matter for small polynomials
            next(pp) = cc[N];
            next(dp) = (0.0,0.0);
            pause;
            for(j=0..N-1) {
                next(dp) = cadd(cmul(dp,root[i]),pp);
                next(pp) = cadd(cmul(pp,root[i]),cc[N-1-j]);
                pause;
            }
            // to avoid division by zero, we cheat a bit
            let(dp1 = (dp==(0.0,0.0) ? (eps*eps,eps*eps) : dp)) 
            next(root[i]) = csub(root[i],cdiv(pp,dp1));
            next(x) = root[i];
            pause;
        } while(sqabs(csub(x,root[i])) > eps*eps);        
    }

    emit(rdy);
}
drivenby p1 {
    // note: x^5 - 15x^4 + 85x^3 - 225x^2 + 274x - 120 = (x-1)(x-2)(x-3)(x-4)(x-5)
    c = [-120.0,+274.0,-225.0,+85.0,-15.0,1.0];
    await(rdy);
}
drivenby p2 {
    // note: x^5-2x^4-3x^3+6x^2+2x-4 = (x-2)(x-1)(x+1)(x-sqrt(2))(x+sqrt(2))
    c = [-4.0,+2.0,+6.0,-3.0,-2.0,1.0];
    await(rdy);
}
drivenby p3 {
    // note: x^5+6x^4+14x^3+24x^2+33x+18 = (x+3)(x+2)(x+1)(x^2+3)
    c = [+18.0,33.0,24.0,14.0,6.0,1.0];
    await(rdy);
}
drivenby p4 {
    // note: x^5-2x^4+2x^3-4x^2-3x+6 = (x-2)(x-1)(x+1)(x^2+3)
    c = [+6.0,-3.0,-4.0,+2.0,-2.0,1.0];
    await(rdy);
}