// ************************************************************************** //
//                                                                            //
//    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 by reducing the  //
// problem to the computation of the eigenvalues of the companion matrix of   //
// the polynomial. The eigenvalues are then computed by the QR-iteration.     //
// Since the companion matrix is a Hessenberg matrix (it has zeros everywhere //
// below the diagonal except for the first subdiagonal), also the matrices    //
// that appear in the QR-decomposition are Hessenberg matrices, and so is the //
// result of the procedure EigenvaluesQR. The algorithm below has however the //
// drawback that convergence is only guaranteed if the absolute values of the //
// eigenvalues are pairwise distinct.                                         //
// ************************************************************************** //

package Algorithms.Polynomials;
import  Algorithms.Matrices.*;

macro N = 5;

// numeric precision
macro eps = 1.0e-8;
macro rabs(x) = (x<0 ? -x : +x);

module PolynomialRootsQR([N+1]real ?c,[N](real * real) root,event !rdy) {
[N][N](real * real) a;

// determine the companion matrix a
for(i=0..N-1)
for(j=0..N-1) {
let(j0=N-1-j)
a[i][j].0 = (i==0 ? -c[j0]/c[N] : (i-1==j ? 1.0 : 0.0));
a[i][j].1 = 0.0;
}

// perform QR-iteration to compute its eigenvalues
EigenvaluesShiftQRC(a,a,root,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);
}