(REAL t = 1; REAL s =  1/(t*t), s1 = -2/(t*t*t);
OP * = ([]REAL a, []REAL b)REAL: (REAL s := 0; FOR j TO 4 DO s +:= a[j]*b[j] OD; s);
[4, 4] REAL g, gu;
FOR i TO 4 DO FOR j TO 4 DO g[i, j] := IF i=j THEN IF i=1 THEN -s
  ELSE s FI ELSE 0 FI OD OD;
FOR i TO 4 DO FOR j TO 4 DO gu[i, j] := IF i=j THEN 1/g[i, j] ELSE 0 FI OD OD;
[4, 4, 4] REAL dg, ch, chu;
FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO dg[i, j, k] := IF i > 1 THEN 0 ELSE
     IF j=1 ANDF k=1 THEN -s1 ELSE IF j=k THEN s1 ELSE 0 FI FI FI OD OD OD;
FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO ch[i, j, k] :=
  0.5 * (dg[j, k, i] + dg[k, i, j] - dg[i, j, k]) OD OD OD;
PROC p4 = (STRING m, [,,,]REAL a)VOID: (print(("p4: "+m, newline));
  FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO FOR l TO 4 DO
  IF REAL r = a[i, j, k, l]; r /= 0 THEN
    print((i-1, j-1, k-1, l-1, "  ", r, newline)) FI OD OD OD OD);
PROC p2 = (STRING m, [,]REAL a)VOID: (print(("p2: "+m, newline));
  FOR i TO 4 DO FOR j TO 4 DO print(a[i, j]) OD; print(newline) OD);
FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO chu[i, j, k] := gu[i, ]*ch[ , j, k] OD OD OD;
(0=0 | print(("Christoffel", newline));
  FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO print(chu[i, j, k]) OD;
    print(newline) OD; print(newline) OD);
(0=1 | p2("gu", gu));
PROC dchu = (INT i, INT j, INT k, INT l)REAL: IF i>0 THEN 0
  ELSE -chu[j, k, l]/t FI;
[4, 4, 4, 4] REAL rie, riel;
FOR m TO 4 DO FOR s TO 4 DO FOR a TO 4 DO FOR b TO 4 DO
  rie[m, s, a, b] := - dchu(b, m, s, a) + dchu(a, m, s, b)
      + chu[m, , a]*chu[ , s, b] - chu[m, , b]*chu[ , s, a] OD OD OD OD;
FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO FOR l TO 4 DO riel[i, j, k, l] :=
  g[i, ]*rie[ , j, k, l] OD OD OD OD;
[4, 4] REAL rie2;
FOR i TO 4 DO FOR j TO 4 DO rie2[i, j] :=
   (REAL s := 0; FOR k TO 4 DO s +:= rie[k, i, j, k] OD; s) OD OD;
p2("Rij = ", rie2);
REAL rie0 = (REAL s := 0; FOR j TO 4 DO s +:= gu[j, ]*rie2[j, ] OD; s);
print(("R = ", rie0)); p4("riel", riel);
[4, 4, 4, 4] REAL c; # Weyl #
FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO FOR l TO 4 DO c[i, j, k, l] := riel[i, j, k, l]
  + (g[i, l]*rie2[k, j] + g[j, k]*rie2[i, l] - g[i, k]*rie2[l, j] - g[j, l]*rie2[k, i])/(4-2)
  + (g[i, k]*g[l, j] - g[i, l]*g[k, j])*rie0/((4-1)*(4-2)) OD OD OD OD;
p4("c", c);

(0=1| # Decomposition of Riemann tensor #
[4, 4, 4, 4] REAL su;
FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO FOR l TO 4 DO su[i, j, k, l] :=
   (rie0/12)*(g[i, k]*g[j, l] - g[i, l]*g[j, k]) OD OD OD OD;
(0=1|p4("su", su));

[4, 4] REAL su2;
FOR i TO 4 DO FOR j TO 4 DO su2[i, j] := rie2[i, j] - rie0*g[i, j]/4 OD OD;
p2("su2", su2);

[4, 4, 4, 4] REAL e;
FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO FOR l TO 4 DO e[i, j, k, l] := (1/(4-2))*
   (g[i, k]*su2[j, l] - g[i, l]*su2[j, k] + g[j, l]*su2[i, k] - g[j, k]*su2[i, l])
OD OD OD OD;
p4("e", e);

[4, 4, 4, 4] REAL c;
  FOR i TO 4 DO FOR j TO 4 DO FOR k TO 4 DO FOR l TO 4 DO c[i, j, k, l] :=
    riel[i, j, k, l] - su[i, j, k, l] - e[i, j, k, l] OD OD OD OD;
p4("wiki c", c))
)