# Algol 68 #

MODE R = REAL;
MODE C = LONG COMPLEX;
MODE LR = LONG REAL;
R small = 32 * ABS SHORTEN long small real;
PROC (REF FILE) VOID nl = newline;

# determinant from ../determ/d.a #
PROC determ = (REF [,]C m)C: (INT k = 1 UPB m; 
IF k /= 2 UPB m THEN print("not square"); GOTO hell FI;
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] * determ(m[2:k, 2:k])) FI);
PROC det = ([,]C m)C: ([1 UPB m, 1 UPB m]C z := m; determ(z));

#  Matrix Inversion
Original at ../GR/inv.a #
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 := SHORTEN ABS am[i, i];
     FOR j FROM i+1 TO n DO R ae = SHORTEN ABS 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 FROM i 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;
   am[,n+1:2*n]);

# Matrix multiplication -- not now in use. #
PROC mq = ([,]C a, [,]C b)[,]C: (INT y = 1 UPB a; [y,y]C z;
  FOR j TO y DO FOR k TO y DO C f := 0; FOR i TO y DO f +:= a[j,i]*b[i,k] OD; 
    z[j,k] := f OD OD; z);

MODE P = [~] C; # polynomial - Ascending powers in sequence #
# Evaluate a polynomial #
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));

# Take the derivitive of a 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;

PROC sm = (C a)R: (R x = ABS SHORTEN re OF a, y = ABS SHORTEN im OF a; (x<y|y|x));

PROC rr = R: nextrandom - 0.5;

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

MODE DR = STRUCT(P q, P r);

#Synthetic division of two polynomials, 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 sm((r OF sol)[1]) > small THEN print(("Failure", nl))
        ELSE mp := q OF sol FI OD; rts);

PROC lagr = (INT n)[,]C: ([n,n]C p;
# Choose a hexagonal set of n points about the origin of the complex plane and fill array p with their powers, 0 thru n-1. #
LR w = longsqrt(3)/2; C tr = (1/2) I w; INT k := 2;
p[1,1] := 1; FOR j FROM 2 TO n DO p[1,j] := 0 OD; 
FOR j DO C st := j, head := (-1/2) I w;
  TO 6 DO (C s := st; TO j DO p[k,1] := 1;
    FOR j TO n-1 DO p[k,j+1] := p[k,j]*s OD;
    s +:= head; k +:= 1; (k>n|ex) OD); st *:= tr; head *:= tr OD OD; ~
EXIT ex: p);

# Compute the characteristic polynomial. #
PROC cp = ([,]C mat)[]C: (INT n = 1 UPB mat; INT m = n+1;
  [,]C l = lagr(m); [,]C li = inv(l);
  [m]C v; FOR k TO m DO [n,n]C r := mat;
    FOR j TO n DO r[j,j] -:= l[k,2] OD; v[k] := det(r) OD;
  [m]C pc; FOR y TO m DO C f := 0;
     FOR z TO m DO f +:= li[y,z]*v[z] OD; pc[y] := f OD;
  pc);

# Compute all eigenvalues #
PROC aeval = ([,]C m)[]C: allr(cp(m));

# Compute all eigenvalues and eigenvectors. #
PROC tot = ([,]C m, REF[]C vals, REF[,]C vecs) VOID: (INT n = 1 UPB m; vals := aeval(m);
  print(("M=", nl, m, nl, "eigenvalues=", nl, vals, nl));
  FOR j TO n DO [n,n]C mm := m; FOR k TO n DO mm[k,k] -:= 1.e-10+vals[j] OD; [,]C nv = inv(mm);
  print(("mm=", nl, mm, nl, "nv=", nl, nv, nl)); print(("prd=", nl, mq(mm, nv), nl));
  (INT i := 0; R mx := 0; FOR j TO n DO FOR k TO n DO
    (R nm = sm(nv[j,k]); (nm>mx|mx:=nm; i:=j)) OD OD;
  FOR k TO n DO vecs[j,k] := nv[i,k]/mx OD) OD);

# Test #
PROC evls = ([,]C m)VOID: print(("m:", nl, m, nl, "evals:", aeval(m), nl));
PROC tst = ([,]C m) VOID: (INT n = 1 UPB m; [n]C vals; [n,n]C vecs; tot(m, vals, vecs);
  FOR j TO n DO FOR k TO n DO C s := 0; FOR i TO n DO s +:= m[k,i]*vecs[j,i] OD;
    (ABS (vecs[j,k]*vals[j] - s) > 1.e-20 |
     print(("Bad", j, k, s, nl, m, nl, vals, nl, vecs, nl)); hell) OD OD);

(3|([,]C mx = lagr(5); print(mx);
  print(nl); [,]C m = inv(mx); print(m);
  C dmx = det(mx), dm = det(m);
  print((nl, dm, dmx, nl));
  print(dmx*dm)),
tst(((1, 0, 0), (0, 2, 0), (0, 0, 3))),
tst(((0, 1), (2, 0))),
tst(((3, 0), (0, 5))),
tst(((1, 0), (0, 0 I 1/10))),
tst(((0, 1), (-1, 0))),
tst(((3, 0, 0), (0, 3, 0), (0, 0, 3))),
print(("roots are ", allr((0, 0, 4))))
)
EXIT hell: print((nl, "The Bitter End"))