MODE R = REAL;
PROC eigen = (REF[,] R m, REF[]R l, REF[,] R vc) VOID: (
  INT n = 1 UPB m; [n,n] R m2 := m;
  FOR i TO n DO FOR j TO n DO vc[i, j] := (i=j|1|0) OD OD;
  WHILE R mod := 0; INT i, j; FOR k TO n DO FOR m FROM k+1 TO n DO R q = ABS m2[k,m];
    IF q > mod THEN mod:=q; i:=k; j:=m FI OD OD; mod > 1.e-14 DO 
    R th = 0.5*(R d = m2[i,i] - m2[j,j]; (d=0|pi/2|atan(2.0*m2[i,j]/d)));
    R c = cos(th), s = sin(th);
    PROC twst = (REF[,]R m)VOID: FOR k TO n DO R t = c*m[i,k]+s*m[j,k];
      m[j,k] := -s*m[i,k]+c*m[j,k]; m[i,k]:=t OD;
    FOR k TO n DO R t = c*m2[k,i] + s*m2[k,j];
      m2[k,j] := -s*m2[k,i] + c*m2[k,j]; m2[k,i] := t OD;
    twst(m2); twst(vc) OD;
  FOR j TO n DO l[j] := m2[j,j] OD);

INT s = 38;
[s,s]R m;
FOR j TO s DO FOR k FROM j TO s DO m[j, k] := next random; m[ k, j] := m[j, k] OD OD;
[s]R evl; [s, s]R evc;
eigen(m, evl, evc);
R me := 0;
# below: i is which eigenvector we consider;
j is which component of that vector we consider;
k is which summand of that component. #
FOR i TO s DO FOR j TO s DO R t := 0;
  FOR k TO s DO t +:= m[j, k] * evc[i, k] OD;
  R er = ABS(evl[i]*evc[i, j] - t);
  IF er > me THEN me := er FI OD OD;
  print((me, newline, evl))