MODE R = REAL;
MODE S = STRUCT(R c, R m);
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) |
     r(b, j-1) - 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=r(1, n); (c OF s/p, m OF s/(sqrt(b)*p*(n+1))))));
       print(("S= ", c((2, 2)) ))
 # (2, 2) → 1/8, 0.0147313912747197
    (4/3, 2/3) → 1/2, 0.121120348781239
    (2/3, 4/3) → 1/2, 0.121120348781239
    (1/3) → 1, 2.5
    (2) → 1/2, 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)) → 1/2, 0.117274273429143
    #