MODE R = REAL;
MODE S = STRUCT(R c, R m);
PROC rr = (STRING m, R r)R: (print((m, r, newline)); r);
PROC rs = (STRING m, S r)S: (print((m, r, newline)); r);
OP - = (S a, b)S: (c OF a - c OF b, m OF a - m OF b);
PROC c = ([]R a)S: (INT n = UPB a;
  PROC r = (R b, INT j)S: (b<0 | (0, 0) | (j=0 | R t=b**n; (t, t*b) |
     rs("s0=", r(b, j-1)) - rs("s1=", r(b - a[j], j-1))));
  (R p:=1, b:=0; FOR j TO n DO p *:= a[j]*j; b +:=a[j]*a[j] OD;
      (S s=rs("s=", r(1, n));
      (c OF s/p, rr("q=", rr("m=", m OF s)/(sqrt(b)*p*(n+1)))))));

print(("S= ", c((1/3)) ))
 # (2, 2) → 1/8, 0.0147313912747197
    (4/3, 2/3) → 0.5, 0.121120348781239
    (2/3, 4/3) → 0.5, 0.121120348781239
    (1/3) → 1, 2.5
    (2) → 0.5, 0.125
    (R t = 2; (t, t, t)) → 1/48, 0.00150351632601465
    (R t = 1; (t, t, t)) → 1/6, 0.0240562612162344
    (R t = 2/3; (t, t, t)) → 0.5, 0.117274273429143
    #