MODE R = REAL;
MODE S = STRUCT(R c, [~]R m);
OP * = (R s, []R r)[]R: (
      INT n = UPB r; [n]R t; FOR j TO n DO t[j] := s*r[j] OD; t);
PROC pow = (R x, INT k)R: (R t := 1; TO k DO t *:= x OD; t);
PROC c = ([]R a)S: (INT n = UPB a;
  S bt = (R p:=1, b := 0; FOR j TO n DO p *:= a[j]*j; b +:= a[j]*a[j] OD;
             (1/p, ([n]R t; FOR j TO n DO t[j] := 1/(p*a[j]*(n+1)) OD; t)));
  PROC r = (R b, INT j)S: (b<0 | (0, ~) | (j=0
      | (R t = pow(b, n); (t*c OF bt, b*t*m OF bt))
      | (S s0 = r(b, j-1), s1 = r(b - a[j], j-1); (c OF s1 > 0 |
         (c OF s0 - c OF s1, ([n]R t; FOR k TO n
             DO t[k] := (m OF s0)[k] - (m OF s1)[k] OD;
             t[j] -:= (c OF s1); t)) | s0)))); r(1, n));

PROC cn = ([]R a)S: (INT n = UPB a; R p := 0, q := 0;
  FOR j TO n DO (a[j] < 0 | q | p) +:= a[j] OD;
  IF [n]R na; p < 1 THEN (1, (FOR j TO n DO na[j] := 1/2 OD; na))
  ELSE BOOL b = p + q > 2;
     FOR j TO n DO
       na[j] := (b = (a[j] > 0) | a[j]/(a[j] - 1) | a[j]/(1 - a[j])) OD; (42, na) FI);
  
 (0=0 | print(("w=", cn((- 1/6, 4/3))))
 
 | (0=0 | print(("S= ", c( (2, -2) ))) |
 
     []STRUCT([]R t, S a) t = (
    ((2, 2), (1/8, (1/48, 1/48))),
    ((2, -2), (1/8, (1/48, 5/48))),
    ((-2, 2), (1/8, (5/48, 1/48))),
    ((-2, -2), (1/8, (0, 0))));

FOR j TO UPB t DO []R tc = t OF t[j];
   FOR k TO UPB tc DO print(tc[k]) OD; print(newline);
     (S s = cn(tc); print(("c = ", c OF s));
     FOR k TO UPB m OF s DO print((m OF s)[k]/c OF s) OD; print(newline);
     IF (ABS (c OF s - c OF a OF  t[j])) > 1.e-14
       THEN print(("Foul c: ", j, c OF s,  c OF a OF  t[j], newline)) FI;
     FOR k TO UPB m OF s DO IF ABS ((m OF s)[k] - (m OF a OF t[j])[k]) > 1.e-13
       THEN print(("Foul m: ", j, k, (m OF s)[k], (m OF a OF t[j])[k], newline))FI OD) OD))