INT n = 4;
MODE V = [n]REAL;
MODE TEN = [n,n]REAL;
TEN m2;
[20]V ma;
first random(45);
PROC rx = REAL: 2*next random-1;
FOR j TO UPB ma DO FOR k TO n DO ma[j][k] := rx OD OD;
FOR i TO n DO FOR j TO n DO m2[i,j] := 0 OD OD;
FOR k TO UPB ma DO
  FOR i TO n DO FOR j TO n DO m2[i,j] +:= ma[k][i]*ma[k][j] OD OD OD;
pnt("moments", m2);
PROC pv = (V x)VOID:
    (FOR j TO n DO print(fixed(x[j], 12, 5)) OD; print(newline));
PROC pnv = (STRING n, V  x)VOID: (print((n, newline)); pv(x));
PROC pnt = (STRING m, TEN x)VOID: (print((m, newline));
         FOR j TO n DO pv(x[j,]) OD);
PROC ip = (V a, b)REAL: (REAL s:=0; FOR j TO n DO s +:= a[j]*b[j] OD; s);
first random(43);
PROC sum = (PROC (INT)REAL pr)REAL: (REAL sm := 0;
  FOR j TO n DO sm +:= pr(j) OD; sm);
TEN omet;
FOR i TO n DO FOR j TO n DO omet[i,j] := (i<j|rx|:i=j|0|-omet[j,i]) OD OD;
pnt("omet", omet);
(TEN momen; #From first principles#
  FOR i TO n DO FOR j TO n DO momen[i,j] := 0 OD OD;
  FOR k TO UPB ma DO V v; FOR i TO n DO v[i] :=
    (0=0|ip(omet[,i],ma[k][])| sum((INT a)REAL: (omet[a,i]*ma[k][a]))) OD;
    FOR i TO n DO FOR j TO n
      DO momen[i,j] +:= ma[k][i]*v[j] - ma[k][j]*v[i] OD OD OD;
  pnt("primitive L", momen));
([n,n,n,n]REAL it2;
FOR a TO n DO FOR b TO n DO FOR c TO n DO FOR d TO n DO
  it2[a,b,c,d] := (b=d|1|0)*m2[a,c] + (a=c|1|0)*m2[b,d] OD OD OD OD;
(TEN myl;
FOR j TO n DO FOR k TO n DO
   myl[j,k] := sum((INT i)REAL: sum((INT v)REAL: it2[i,v,j,k]*omet[i,v])) OD OD;
pnt("L in first form", myl));
(0<0|FOR i TO n DO FOR j TO n DO pnt((j=1|"---"|""), it2[i,,j,]) OD OD));
(TEN l2;
 FOR j TO n DO FOR k TO n DO
    l2[j,k] := ip(m2[j,], omet[,k]) + ip(m2[k,], omet[j,]) OD OD;
 pnt("L in second form", l2));
(TEN l3;
 FOR i TO n DO FOR j TO n DO
 l3[i,j] := ip(omet[i,], m2[j,]) - ip(omet[j,], m2[i,]) OD OD;
 pnt("L in third form", l3))