MODE C = LONG COMPLEX;
MODE R = REAL;
MODE P = [~] C;
OP SH = (LONG REAL x)REAL: SHORTEN x;
R small = 16 * SH long small real; R s2 = small*small;
PROC ev = (P p, C x)C: (INT d := UPB p; C v := p[d];
   WHILE (d -:= 1)>0 DO v := x*v+p[d] OD; v);
PROC der = (P p)P: IF INT d = UPB p; d < 2 THEN print("Constant Polynomial"); ()
    ELSE [d-1] C v; FOR i TO d-1 DO v[i] := i*p[i+1] OD; v FI;
PROC m2 = (C a)R: (R x = SH re OF a, y = SH im OF a; x*x + y*y);
PROC root = (P p)C: (C z := 0; P dv = der(p); WHILE
   C e = ev(p, z); C f = ev(dv, z); C d := e/f; R h = m2(e);
   WHILE C tp = z - d; R i = m2(ev(p, tp)); (i<=h|z := tp; 1<0|0<1)
     DO d *:= 0.5 I (0.1 * nextrandom - 0.05) OD; m2(e) > s2 DO ~ OD; z);

#Tests#
# 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 s := (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));

# Three diverse tests and demos.#
##(1|(
PROC tst = (P p)VOID: (C r = root(p);
  print(("p=", p, newline, "r=", r, newline, "v=", ev(p, r), newline)));
##(5|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)),
~)),

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

(
#Synthetic division of two polynimials, returns quotient and remainder.#
MODE DR = STRUCT(P q, P r);
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 m2((r OF sol)[1]) > small THEN print(("Failure", newline))
        ELSE mp := q OF sol FI OD; rts);

# 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] := rnd(s) I rnd(s) OD; p[n+1] := 1;
  []C rts = allr(p), q = pfr(rts);
  FOR j TO n DO R e = m2(p[j] - q[j]); print(e);
  IF e > small THEN print(("bad", j, p[j], q[j], newline)) FI OD);

tst(5)))