MODE C = LONG COMPLEX;
MODE RR = LONG REAL;
MODE R = REAL;
MODE P = [~] C;
PROC pp = (STRING m, P a)VOID: (print((m, newline));
  FOR k FROM LWB a TO UPB a DO print((whole(k, 4)+" ", a[k], newline)) OD);
PROC pm = (STRING m, [,]C mm)VOID: (print((m, newline));
  FOR j TO UPB mm DO pp("", mm[j,]) OD);
OP AB = (LONG REAL x)REAL: ABS SHORTEN x;
R small = 32 * AB long small real;
PROC rr = R: nextrandom - 0.5;

# derivative of polynomial #
PROC der = (P p)P: IF INT d = UPB p; d = 0 THEN ()
    ELSE [d-1] C v; FOR i TO d-1 DO v[i] := i*p[i+1] OD; v FI;

# p(x) for polynomial p and argument x #
PROC ev = (P p, C x)C: (INT d := UPB p; (d=0|0|
   C v := p[d]; WHILE (d -:= 1)>0 DO v := x*v+p[d] OD; v));

# Crude magnitude of complex number #
PROC cm = (C a)R: (R x = AB re OF a, y = AB im OF a; (x<y|y|x));

PRIO MX = 4;
OP MX = (REF R x, R y)VOID: (x<y|x:=y);

# Crude vector magnitude #
PROC cvm = (P x)R: (R m := 0; FOR j TO UPB x DO m MX cm(x[j]) OD; m);

# Crude Matrix magnitude #
PROC cmmg = ([,]C m)R: (R x := 0; FOR j TO UPB m DO x MX cvm(m[j,]) OD; x);

# Scalar times Vector #
OP * = (RR s, P v)P: (INT y = UPB v; [y]C a;
  FOR j TO y DO a[j] := (s*re OF v[j]) I (s*im OF v[j]) OD; a);

# Vector subtraction #
OP - = (P a, b)P: (INT n = UPB a; [n]C an; FOR j TO n DO an[j] := a[j] - b[j] OD; an);

PROC root = (P p)C: (C z := 0; P dv = der(p);
   l: (C e = ev(p, z); R h = cm(e);
   IF h < small THEN GOTO done FI;
   C f = ev(dv, z);
   IF f = 0 THEN z +:= rr I rr; GOTO l FI; C d := e/f; 
   m: (C tp = z - d;
   IF cm(ev(p, tp))>h
     THEN d *:= 0.5 I (0.1 * rr); GOTO m FI;
   z := tp; GOTO l);
   done: z));

#Tests#
# RC4: Pseudo Random Number Generator #
MODE RS = STRUCT([0:255]CHAR s, INT i, j);
OP MD = (INT j)INT: (j<256|j|MD (j-256));
PROC init = (STRING key)RS: ([0:255]CHAR s;
  FOR j FROM 0 TO 255 DO s[j] := REPR j OD;
  INT len = UPB key; INT j := 0;
  FOR i FROM 0 TO 255 DO CHAR t = s[i];
    j := MD (ABS t + j + ABS key[(i %* len) + 1]);
    s[i] := s[j]; s[j] := t OD; (s, 0, 0));
PROC nr = (REF RS s)INT: (i OF s := MD(i OF s+1);
    CHAR a = (s OF s)[i OF s]; INT ai = ABS a;
    j OF s := MD (j OF s + ai); CHAR b = (s OF s)[j OF s];
    (s OF s)[i OF s] := b; (s OF s)[j OF s] := a;
    ABS((s OF s)[MD(ai + ABS b)]));

RS st := init("Ambidexterously");
PROC nextrandom = REAL: (REAL x := 0; TO 7 DO x := 256*x + nr(st) OD; x*(2.0**-56));

# pfr(n) constructs a polynomual with n given complex roots. #
PROC pfr = (P rts)P: (INT d = UPB rts; [0:d]C r; r[0] := 1;
   FOR j TO d DO r[j] := 0 OD;
   FOR k TO d DO C rt := rts[k]; FOR j FROM d DOWNTO 1 DO
       r[j] := r[j-1] - rt*r[j] OD; r[0] *:= -rt OD; r[0:d]);

# Algol68 Normal Random Deviate Generator #
MODE RNDS = STRUCT(BOOL f, REAL v);
RNDS szx := (1<0, ~); 
PROC rnd = (REF RNDS s)REAL: (f OF s | f OF s := 1<0; v OF s
  | REAL x = nextrandom*2 - 1, y = nextrandom*2 - 1;
  REAL r2 = x*x + y*y;
    (r2 > 1.0 | rnd(s) | REAL q = sqrt(-2*ln(r2)/r2);
      s := (0<1, x*q); y*q));

# Complex rnd #
PROC crnd = C: rnd(szx) I rnd(szx);

MODE DR = STRUCT(P q, P r);
#Synthetic division of two polynimials, returns quotient and remainder.#
PROC div = (P n, d)DR: (
INT nn = UPB n, nd = UPB d - 1; [nn]C mn := n;
   [nn-nd]C q; C qr = 1/d[nd+1];
   FOR j FROM UPB q DOWNTO 1 DO C qs = qr*mn[j+nd]; q[j] := qs;
   FOR k TO nd+1 DO mn[k+j-1] -:= qs*d[k] OD OD; (q, mn[1:nd]));

# Routine allr finds all roots of a polynomial. #
PROC allr = (P p)P: (INT n = UPB p; [n-1]C rts; FLEX [n]C mp := p;
   FOR j TO n-1 DO C rt = root(mp); rts[j] := rt;
   DR sol = div(mp, (-rt, 1));
     IF cm((r OF sol)[1]) > small THEN print(("Failure", newline))
        ELSE mp := q OF sol FI OD; rts);

# Random Square Complex Matrix #
PROC rm = (INT k)[,]C: ([k,k]C x;
   FOR i TO k DO FOR j TO k DO x[i, j] := crnd OD OD; x);

# Matrix Inversion #
PROC inv = ([,]C m) [,]C: (INT n = 1 UPB m;
   IF n /= 2 UPB m THEN print("not square") FI;
   [n, 2*n]C am; FOR i TO n DO FOR j TO n DO
      am[i, j] := m[i, j]; am[i, j+n] := (i=j|1|0) OD OD;
   FOR i TO n DO INT im := 0; R mv := cm(am[i, i]);
     FOR j FROM i+1 TO n DO R ae = cm(am[j, i]);
       IF mv < ae THEN mv := ae; im := j FI OD;
     IF mv = 0.0 THEN print("Singular") FI;
     IF im > i THEN FOR n FROM i TO 2*n DO
        C s = am[i, n]; am[i, n] := am[im, n]; am[im, n] := s OD FI;
   (C r = 1.0 / am[i, i]; FOR j TO 2*n DO am[i, j] *:= r OD);
   FOR k FROM i+1 TO n DO C f = am[k, i]; 
     FOR j FROM i+1 TO 2*n DO am[k, j] -:= f * am[i, j] OD OD OD;
   FOR i TO n DO FOR j FROM i+1 TO n DO C p = am[i, j];
     FOR k FROM i+1 TO 2*n DO am[i, k] -:= am[j, k] * p OD OD OD;
   ([n, n]C a;
   FOR i TO n DO FOR j TO n DO a[i, j] := am[i, n+j] OD OD; a));

# Determinant #
PROC det = (REF [,]C m)C: (INT k = 1 UPB m; 
IF k = 1 THEN m[1,1]
ELSE FOR j TO k DO IF m[j,1] /= 0.0 THEN 
   IF j /= 1 THEN FOR n TO k DO
      C t = m[j, n]; m[j, n] := m[1, n]; m[1, n] := -t OD FI;
   GOTO zz FI OD; 0 EXIT
zz: (C t = 1.0 / m[1,1];
  FOR i FROM 2 TO k DO REF []C q = m[i,]; C a = t*q[1];
     FOR r FROM 2 TO k DO q[r] -:= a*m[1, r] OD OD;
  m[1,1] * det(m[2:k, 2:k])) FI);

# determanent of M-li #
PROC sh = ([,]C mt, C l)C: (INT n = UPB mt; [n,n]C m := mt;
  FOR j TO n DO m[j,j] -:= l OD; det(m));

# inverse of M-li #
PROC ish = ([,]C mt, C l)[,]C: (INT n = UPB mt; [n,n]C m := mt;
  FOR j TO n DO m[j,j] -:= l OD; inv(m));

# squared mangitude #
PROC m2 = (C x)RR: x[re]*x[re] + x[im]*x[im];

# squared vector magnitude #
PROC vm2 = (P v)RR: (RR s := 0; FOR j TO UPB v DO s +:= m2(v[j]) OD; s);

# Scalar times Vector #
OP * = (C s, P x)P: (INT y = UPB x; [y]C a; FOR j TO y DO a[j] := s*x[j] OD; a);

# Scalar times Matrix #
OP * = (RR s, [,]C mx)[,]C: (INT m = 1 UPB mx, n = 2 UPB mx; [m,n]C a;
  FOR i TO m DO FOR j TO n DO a[i,j] := s*mx[i,j] OD OD; a);

# Matrix times Vector #
OP * = ([,]C m, []C v)[]C: (INT n = 2 UPB m; [n]C a;
  FOR j TO n DO C s := 0; FOR k TO UPB m DO s +:= m[j,k]*v[k] OD;
     a[j] := s OD; a);

# Matrix times Matrix #
OP * = ([,]C a, [,]C b)[,]C: (INT i = UPB a, j = UPB b, k = 2 UPB b;
  [i,k]C an; FOR u TO i DO FOR v TO k DO C s := 0;
     FOR w TO j DO s +:= a[u,w]*b[w,v] OD; an[u,v] := s OD OD; an);

# Inner product of two vectors #
PROC ip = (P x, y)C: (C s := 0; FOR i TO UPB x DO s +:= x[i]*CONJ y[i] OD; s);

# Gram Schmidt #
PROC gs = ([,]C p)[,]C: (INT m = 1 UPB p, n = 2 UPB p;
  MODE R = LONG REAL;
  [m,n]C q := p; FOR j TO m DO FOR k TO j-1 DO 
      C f = ip(q[j,], q[k,]);
      FOR i TO n DO q[j,i] -:= f*q[k,i] OD OD;
    R r = 1/long sqrt(vm2(q[j,])); FOR k TO n DO q[j, k] *:= r OD OD; q);

# Compute Characteristic Polynomial #
PROC chrply = ([,]C m)P: (INT n = UPB m; INT s = ROUND sqrt(n+1);
  [n+1,n+1]C z; [n+1]C t; INT u := 1;
  FOR i FROM -s%2 DO FOR j FROM -s%2 TO s%2 DO C x = i I j;
     t[u] := sh(m, x);
     (C v := 1; FOR w TO n+1 DO z[w, u] := v; v *:= x OD);
      IF u = n+1 THEN GOTO ex FI; u+:= 1 OD OD; ~ EXIT ex: ([,]C p = inv(z);
     ([n+1]C an; FOR j TO n+1 DO C f := 0; FOR k TO n+1 DO f +:= t[k]*p[k,j] OD;
        an[j] := f OD; an)));

# Compute engenvalues of matrix #
PROC evals = ([,]C m)P: allr(chrply(m));

# Robust normalization #

PRIO NRM = 5;
#OP NRM = (P v2)P: (P v1 = C(LENG (1/cvm(v2)))*v2; (1/long sqrt(ip(v1, v1)))*v1);#

# Report quality of putative eigenvector. #
PROC qev = ([,]C m, P v2, C ev)R: (P
   v1 = LENG (1/cvm(v2))*v2,
   v = (1/long sqrt(vm2(v1)))*v1,
   u = m*v,
   w = ev*v;
   pp("v1", v1); pp("v", v); pp("u", u); pp("w", w);
   cvm(w-u));

# Test for Jordan Normal Form #
PROC jorp = ([,]C m)BOOL: (INT n = UPB m;
   FOR j TO n-2 DO FOR k FROM j+2 TO n DO (m[j,k] /= 0 | GOTO bad) OD OD;
   FOR j FROM 2 TO n DO FOR k TO j-1 DO (m[j,k] /= 0 | GOTO bad) OD OD;
   FOR j TO n-1 DO (m[j,j+1] /= 0 |(m[j,j+1] /= 1 | GOTO bad |
         (m[j,j] /= m[j+1,j+1] | GOTO bad))) OD; TRUE EXIT bad: FALSE);

PROC z = VOID: print("Out of Range");

# 12 diverse tests and demos #
##(1|
(PROC t = (INT n, x, y)VOID:
(IF x>n OR y>n OR x-1>y THEN print(("Foo", newline)) ELSE
[n,n]C a; FOR i TO n DO FOR j TO n DO a[i,j] := (i=j|3|:i+1=j|1|0) OD OD;
a[x,y] := 0.0001;
[n,n]C q; FOR i TO n DO FOR j TO n DO q[i,j] := (i=j|1|0) OD OD;
q[x,y-1] := -0.0001;
[,]C p = inv(q), c = q*a*p;
  pm("a=", a); pm("c = ", c); print(((jorp(c)|"is"|"isnt"), " JNF.")) FI);
#t(5, 1, 5)#
t(5, 2, 5)
),
  
(PROC tst = (P p)VOID:
##(1|(C r = root(p);
  print(("p=", p, newline, "r=", r, newline, "v=", ev(p, r), newline))),
  (P rts = allr(p); FOR j TO UPB rts DO print((rts[j], newline)) OD)|z);
##(7|tst((1, 3)),
tst((-4 I 3, 0.3 I -1, 1 I 0.1, 0.6)),
tst((4, 1, 1)),
tst((1, 2, 1)),
tst((1, 0, -1)),
tst((1, 0, 1)),
tst((1, 0, 0, 0, 0, 0, 0, 0, 1))|z)),

(INT d = 10; [d]C rts; FOR k TO d DO rts[k] := crnd OD; P r = pfr(rts);
  FOR j TO d DO print(("j=", j, ", |ev|=", cm(ev(r, rts[j])),
      ", rt:", newline, rts[j], newline)) OD;
  print((root(r)))),

(
# Choose random nth degree monic polynomial p, find its roots rts,
  compute polynomial with those roots q, compare p with q. #
PROC tst = (INT n)VOID: (
  [n+1]C p; FOR j TO n DO p[j] := crnd OD; p[n+1] := 1;
  []C q = pfr(allr(p)); FOR j TO n DO R e = cm(p[j] - q[j]);
  IF e > small THEN print(("bad", j, e, small, p[j], q[j], newline)) FI OD);

TO 1000 DO tst(6) OD),

print(det(LOC [2,2]C := ((1, 0), (0, 1 I 1)))),
print(inv([,]C((0,0,0,1), (0,0,1,0), (1,0,0,0), (0,1,0,0)))),

(INT b = 30; [,]C f = rm(b), n = inv(inv(f)); 
  FOR i TO b DO FOR j TO b DO IF cm(f[i,j] - n[i,j]) > small THEN
     print(("f[",i,",",j,"]=",f[i,j], newline,
            "n[",i,",",j,"]=",n[i,j], newline)) FI OD OD;
    print("Done Inv Test")),

FOR j FROM 2 TO 5 DO print((j, newline)); [j,j]C z;
  FOR t TO j DO FOR u TO j DO z[t, u] := (t=u|1|0) OD OD;
  pp("Coefficients", chrply(z)) OD,

FOR j FROM 2 TO 9 DO [,]C m = rm(j); P p = chrply(m); C y = crnd;
     print((ev(p, y), newline)); print((sh(m, y), newline)) OD,

(PROC t = ([,]C w)VOID: (INT y = UPB w; [,]C g = gs(w), h = gs(g);
  FOR j TO y DO FOR k TO y DO (cm(g[j,k] - h[j,k]) > small | print("George")) OD OD;
  (C d = det(LOC [y,y]C := g); print((d, m2(d)))));
t((2| rm(4), ((1, 0), (0, -1))))),

([,]C a = rm(5)[1:3,], b = gs(a);
FOR i TO 3 DO FOR j TO 3 DO print(ip(b[i,], b[j,])) OD OD),

(PROC t = (INT y)VOID: (P o = allr(chrply(gs(rm(y)))); pp("roots", o);
  FOR j TO UPB o DO (RR e = m2(o[j]) - 1; ABS e > small
     | print(("off ", e, newline))) OD);
  t(5); t(10); t(11)),

(PROC t = ([,]C m)VOID: (INT n = UPB m; P r = evals(m);
  FOR j TO n DO [,]C mx = ish(m,  r[j]), smx = LENG (1/cmmg(mx))*mx;
  print(("Eigenvalue =", newline, r[j], newline));
  pm("inverse =", mx);
  print((whole(j, 4) + ", det: ", newline, det(LOC [n,n]C := smx), newline));
    FOR i TO n DO pp(whole(i, 3) + " evec?", smx[i,]);
    IF cvm(smx[i,]) > 0.1 THEN print((newline, "qev = ",
       qev(m, smx[i,], r[j]), newline)) FI OD OD);
t((1|rm(3), ((0, -1), (1, 0)), ((1.1, 0), (0, 1))))),

(PROC t = (C d)VOID: ([2,2]C m := ((d, -1), (1, d)); pm("sy", inv(m));
pp("mm", ([,]C((0, -1), (1, 0))) * P(1, 0 I -1)));
t((0 I 1)+ 0.00001*(-1 I 0))),

print ((qev(((0, -1), (1, 0)), (1, 0 I -1), (0 I 1))))

|z)