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