MODE C = LONG COMPLEX;
MODE R = REAL;
MODE P = [~] C;
OP AB = (LONG REAL x)REAL: ABS SHORTEN x;
R small = 32 * AB long small real;
PROC rr = R: nextrandom - 0.5;
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 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));
PROC sm = (C a)R: (R x = AB re OF a, y = AB im OF a; (x<y|y|x));

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; 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));

#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 n = (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("Slee sought");
PROC nextrandom = REAL: (REAL x := 0; TO 7 DO x := 256*x + n(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 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));

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

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

# Three diverse tests and demos.#
##(1|(
PROC tst = (P p)VOID:
##(2|(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 C root = rts[j];
    print(("root=", root, "v=", ev(p, root), newline)) OD; print(newline))|z);
##(tst((1, 3));
tst((1, -3, 1, 1));
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)))),

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

TO 1000 DO tst(4) OD)|z)