```// ************************************************************************** //
//                                                                            //
//    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 improves module is a version of EigenvaluesShiftQR that//
// works on complex numbers.                                                  //
// ************************************************************************** //

macro N = 5;

// -------------------------------------------------------------------------- //
// numeric precision
// -------------------------------------------------------------------------- //
macro eps = 1.0e-8;
// -------------------------------------------------------------------------- //
// macros on real numbers
// -------------------------------------------------------------------------- //
macro rsign(x) = (x<0.0 ? -1.0 : +1.0);
macro rabs(x)  = (x<0.0 ? -x : +x);
macro rroot(x) = exp(x,0.5);
// -------------------------------------------------------------------------- //
// 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;
macro unit = (1.0,0.0);
macro zero = (0.0,0.0);
// -------------------------------------------------------------------------- //

module EigenvaluesShiftQRC([N][N](real * real) ?m,a,[N](real * real) e,event !rdy) {
[N][N](real * real) q,r;                    // matrices q,r of QR-decomposition
(real * real) trc2,det,dis,rtdis,e1,e2,p;   // to determine shift parameter p
real dist1,dist2;   // distances of eigenvectors for Wilkinson shift
nat{N+1} n;         // current size of the matrix (due to deflation)
event QRrdy;        // whether QR-decomposition terminated

// initialize
n = N;
for(i=0..N-1)
for(j=0..N-1)
a[i][j] = m[i][j];

do {
do {

// determine Wilkinson shift p: the idea is to compute the eigenvalues
// of the matrix [[a[N-2][N-2],a[N-2][N-1]],[a[N-1][N-2],a[N-1][N-1]]]
// and to select that eigenvalue that is closer to a[N-1][N-1] as shift p
// note that this is the only step where complex numbers may be introduced
// for real matrices, otherwise the QR-iteration would fail for real matrices
// that have complex eigenvalues

let(a00 = a[n-2][n-2])
let(a01 = a[n-2][n-1])
let(a10 = a[n-1][n-2])
let(a11 = a[n-1][n-1]) {
det = csub(cmul(a00,a11),cmul(a01,a10));
dis = csub(cmul(trc2,trc2),det);
// calling croot directly caused inaccuracies, so the following is better
rtdis = (dis.1==0.0 ? (dis.0<0.0 ? (0.0,rroot(-dis.0)) : (rroot(dis.0),0.0))
: croot(dis));
e2 = csub(trc2,rtdis);
dist1 = sqabs(csub(e1,a11));
dist2 = sqabs(csub(e2,a11));
p = (dist1<dist2 ? e1 : e2);
}

// shift matrix a
for(i=0..N-1)
next(a[i][i]) = csub(a[i][i],p);
pause;

// decompose shifted matrix A_k-p*E into factors Q_k and R_k
DecompQRC(a,q,r,QRrdy);

// compute A_{k+1} := R_k*Q_k+p_k*E and check for convergence
for(i=0..N-1) {
for(j=0..N-1) {
let(p0 = (i==j ? p : zero))
let(new_aij_re = Re(p0) + sum(k=0..N-1) Re(cmul(r[i][k],q[k][j])))
let(new_aij_im = Im(p0) + sum(k=0..N-1) Im(cmul(r[i][k],q[k][j])))
next(a[i][j]) = (new_aij_re,new_aij_im);
}
}
pause;
} while(sqabs(a[n-1][n-2])>eps*eps);
next(n) = n-1;
// eigenvalue e[n-1] has been found, we can deflate the matrix
e[n-1] = a[n-1][n-1];
for(i=0..N-1) {
next(a[i][n-1]) = zero;
next(a[n-1][i]) = zero;
}
pause;
} while(n>1);

e[n-1] = a[n-1][n-1];
emit(rdy);
}
drivenby m0 {
m = [(1.0,0.0),(2.0,0.0),(3.0,0.0),(4.0,0.0),(5.0,0.0)];
m = [(2.0,0.0),(3.0,0.0),(4.0,0.0),(5.0,0.0),(1.0,0.0)];
m = [(3.0,0.0),(4.0,0.0),(5.0,0.0),(1.0,0.0),(2.0,0.0)];
m = [(4.0,0.0),(5.0,0.0),(1.0,0.0),(2.0,0.0),(3.0,0.0)];
m = [(5.0,0.0),(1.0,0.0),(2.0,0.0),(3.0,0.0),(4.0,0.0)];
await(rdy);
}
drivenby m1 {
// note: x^5 - 15x^4 + 85x^3 - 225x^2 + 274x - 120 = (x-1)(x-2)(x-3)(x-4)(x-5)
m = [ (15.0,0.0), (-85.0,0.0), (225.0,0.0), (-274.0,0.0), (120.0,0.0)];
m = [ unit,  zero,  zero,  zero,  zero];
m = [ zero,  unit,  zero,  zero,  zero];
m = [ zero,  zero,  unit,  zero,  zero];
m = [ zero,  zero,  zero,  unit,  zero];
await(rdy);
}
drivenby m2 {
// note: x^5-2x^4-3x^3+6x^2+2x-4 = (x-2)(x-1)(x+1)(x-sqrt(2))(x+sqrt(2))
m = [ (2.0,0.0), (3.0,0.0), (-6.0,0.0), (-2.0,0.0), (4.0,0.0)];
m = [ unit,  zero,  zero,  zero,  zero];
m = [ zero,  unit,  zero,  zero,  zero];
m = [ zero,  zero,  unit,  zero,  zero];
m = [ zero,  zero,  zero,  unit,  zero];
await(rdy);
}
drivenby m3 {
// note: x^5+6x^4+14x^3+24x^2+33x+18 = (x+3)(x+2)(x+1)(x^2+3)
m = [ (-6.0,0.0), (-14.0,0.0), (-24.0,0.0), (-33.0,0.0), (-18.0,0.0)];
m = [ unit,  zero,  zero,  zero,  zero];
m = [ zero,  unit,  zero,  zero,  zero];
m = [ zero,  zero,  unit,  zero,  zero];
m = [ zero,  zero,  zero,  unit,  zero];
await(rdy);
}
drivenby m4 {
// note: x^5-2x^4+2x^3-4x^2-3x+6 = (x-2)(x-1)(x+1)(x^2+3)
m = [ (2.0,0.0), (-2.0,0.0), (4.0,0.0), (3.0,0.0), (-6.0,0.0)];
m = [ unit,  zero,  zero,  zero,  zero];
m = [ zero,  unit,  zero,  zero,  zero];
m = [ zero,  zero,  unit,  zero,  zero];
m = [ zero,  zero,  zero,  unit,  zero];
await(rdy);
}
```