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);
OP + = ([]R x, []R y)[]R: (INT n = UPB x; [n]R t;
   (~(n=UPB y)|print(("Foul =", newline))); FOR j TO n DO t[j]:=x[j]+y[j] OD; t);
PROC na = (INT n, PROC (INT)R p)[]R: ([n]R t; FOR j TO n DO t[j]:= p(j) OD; t);
PROC c = ([]R a)S: (INT n = UPB a;
  S bt = (R p:=1; FOR j TO n DO p *:= a[j]*j OD;
             (1/p, na(n, (INT j)R: 1/(p*a[j]*(n+1)))));
  PROC r = (R b, INT j)S: (b<0 | (0, ~) | (j=0
      | R t = 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 | s0 |
         (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))))); r(1, n));

PROC u = (INT n, R v)[]R: na(n, (INT j)R: v);
PROC sv = ([]R x)R: (R s := 0; FOR j TO UPB x DO s +:= x[j] OD; s);
PROC sml = ([]R x)BOOL: UPB x=0 ORF (ABS x[1] < 0.00000000001 ANDF sml(x[2:]));
(0=1 | []STRUCT([]R t, S a) t = (
    ((2, 2), (1/8, (1/48, 1/48))),
    ((1, 2), (1/4, (1/12, 1/24))),
    ((1, 1), (1/2, (1/6, 1/6))),
    ((1/2, 1/2), (1, (1/2, 1/2))),
    ((4/3, 2/3), (1/2, (13/96, 5/24))),
    ((2/3, 4/3), (1/2, (5/24, 13/96))),
    ((2, 1, 3/2), (1/18, (1/144, 1/72, 1/108))),
    ((1, 3/2, 2), (1/18, (1/72, 1/108, 1/144))),
    ((3/2, 2, 1), (1/18, (1/108, 1/144, 1/72))),
    ((1/3), (1, (1/2, 8))),
    ((2), (1/2, (1/8, 11))),
    (u(3, 2), (1/48, u(3, 1/(3*128)))),
    (u(3, 1), (1/6, u(3, 1/24))),
    (u(3, 2/3), (1/2, u(3, 35/(3*64)))),
    (u(3, 2/5), (47/48, u(3, 185/384))),
    ((1/4.5, 1/5.1, 1/4.7, 1/6.2), (1, u(4, 1/2))),
    ((1/2, 1), (3/4, (1/3, 7/24))),
    ((1/7, 8/7), (13/16, (19/48, 127/384))));

 FOR j TO UPB t DO []R tc = t OF t[j];
   FOR k TO UPB tc DO print(tc[k]) OD; print((" case:", newline));
   (S s = c(tc); print((c OF s, " = content;  "));
     FOR k TO UPB m OF s DO print((m OF s)[k]/c OF s) OD; print((" = centroid.", 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;
 (PROC ip = ([]R x)[]R: (R f = 1/(sv(x)-1); na(UPB x, (INT j)R: x[j]*f));
   PROC ts = ([]R x)VOID: (S a = c(x), b = c(ip(x)); print(("x=", x, newline));
     INT n = UPB x;
     print(("total", c OF a, c OF b, c OF a + c OF b, newline));
     print((m OF a, newline)); print((m OF b, newline));
     []R mq = m OF a + u(n, c OF b) + -1.0* m OF b;
     (~sml(mq + -1.0*u(n, 1/2)) | print(("Foul cent: ", mq, newline))));
   FOR j TO UPB t DO IF c OF a OF t[j] < 1 THEN ts(t OF t[j]) FI OD;
   ts((3, 0.001))));
(0=1|print(("U: ", c((3, 0.7, 0.9, 0.001)), newline));
     print(("V: ", c((3, 0.7, 0.9)), newline)));
(0=0| PROC slv = (R top, PROC(R)BOOL t)R: (R x := 0, b := top;
     WHILE R q = x+b; q>x DO (t(q)| x := q); b *:= 0.5 OD; x);
  IF 1=1 THEN
  FOR j TO 6 DO print((j, "D   den           summit            centroid elev", newline));
    R r = sqrt(j);
    FOR z TO 6 DO R den = (6 - z)/10, dn = (den>1.e-13|den|1.e-13);
      R sl = slv(1, (REAL y)BOOL: c OF c(u(j, 1/(j*y)))<1-dn);
      S s = c(u(j, 1/(j*sl)));
      print((dn, " ", sl*r, r*(sl-(m OF s)[1]/c OF s), newline)) OD OD
  ELSE
  FOR z TO 6 DO R den = (6 - z)/10, dn = (den>1.e-13|den|1.e-13);
    print(("      dim        summit            centroid elev  den=", dn, newline));
    FOR j TO 6 DO R r = sqrt(j);
      R sl = slv(1, (REAL y)BOOL: c OF c(u(j, 1/(j*y)))<1-dn);
      S s = c(u(j, 1/(j*sl)));
      print((j, " ", sl*r, r*(sl-(m OF s)[1]/c OF s), newline)) OD OD FI)