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