{a,b,c,d,e}::Indices;
\delta{#}::KroneckerDelta;
ex:=\delta^{a}_{b}^{c}_{d};
expand_delta(_);


{a,b,c,d,e,m,n,p,q,r,s,t}::Indices;
A_{a b c}::AntiSymmetric;
\delta{#}::KroneckerDelta;
ex:=\delta^{a}_{b}^{c}_{d} A_{a c e};
expand_delta(_);

# Lovelock 2nd order (Gauss-Bonnet)
{a,b,c,d,e,m,n,p,q,r,s,t}::Indices;
\delta{#}::KroneckerDelta;
R_{m n p q}::RiemannTensor;
ex:=\delta_{a m b n c p d q} R_{a b m n} R_{c d p q};
expand_delta(_);
distribute(_);
eliminate_kronecker(_);
meld(_);
substitute(_, $R_{m n m n}=R, R_{m c p c} = R_{m p}, R_{c m c p} = R_{m p}$);

# Lovelock 3rd order
{a#, m#}::Indices;
\delta{#}::KroneckerDelta;
R_{m1 m2 m3 m4}::RiemannTensor;
ex:=\delta_{a1 m1 a2 m2 a3 m3 a4 m4 a5 m5 a6 m6} R_{a1 a2 m1 m2} R_{a3 a4 m3 m4} R_{a5 a6 m5 m6};
expand_delta(_);
distribute(_);
eliminate_kronecker(_);
substitute(_, $R_{m1 m2 m1 m2}=R, R_{m1 a1 m1 a2} = R_{a1 a2}, R_{a1 m1 m1 a2} = -R_{a1 a2}, R_{m1 a1 a2 m1} = -R_{a1 a2}, R_{a1 m1 a2 m1} = R_{a1 a2}$, repeat=True);
sort_product(_)
meld(_);
print(len(ex))

# Lovelock 4th order    
{a#, m#}::Indices;
\delta{#}::KroneckerDelta;
R_{m1 m2 m3 m4}::RiemannTensor;
ex:=\delta_{a1 m1 a2 m2 a3 m3 a4 m4 a5 m5 a6 m6 a7 m7 a8 m8} R_{a1 a2 m1 m2} R_{a3 a4 m3 m4} R_{a5 a6 m5 m6} R_{a7 a8 m7 m8};
expand_delta(_);
distribute(_);
eliminate_kronecker(_);
print(len(ex))
# 2520 terms
substitute(_, $R_{m1 m2 m1 m2}=R, R_{m1 a1 m1 a2} = R_{a1 a2}, R_{a1 m1 m1 a2} = -R_{a1 a2}, R_{m1 a1 a2 m1} = -R_{a1 a2}, R_{a1 m1 a2 m1} = R_{a1 a2}$, repeat=True);
sort_product(_)
print(len(ex))
# 1489 terms
canonicalise(_);
rename_dummies(_);
print(len(ex))
# 25 terms
meld(_);
print(len(ex))  
# 20 terms

# Lovelock 5th order    
{a#, m#}::Indices;
\delta{#}::KroneckerDelta;
R_{m1 m2 m3 m4}::RiemannTensor;
ex:=\delta_{a1 m1 a2 m2 a3 m3 a4 m4 a5 m5 a6 m6 a7 m7 a8 m8 a9 m9 a10 m10} R_{a1 a2 m1 m2} R_{a3 a4 m3 m4} R_{a5 a6 m5 m6} R_{a7 a8 m7 m8} R_{a9 a10 m9 m10};
expand_delta(_)
distribute(_) # need l_applied_no_new_dummies; slow now
eliminate_kronecker(_) # ditto; slow now
# 113400 terms
print(len(ex))
canonicalise(_);  # takes a long time
# 15814 terms
print(len(ex))
rename_dummies(_);
# 452 terms
print(len(ex))
meld(_);
#  terms
print(len(ex))  
substitute(_, $R_{m1 m2 m1 m2}=R, R_{m1 a1 m1 a2} = R_{a1 a2}, R_{a1 m1 m1 a2} = -R_{a1 a2}, R_{m1 a1 a2 m1} = -R_{a1 a2}, R_{a1 m1 a2 m1} = R_{a1 a2}$, repeat=True);


{m,n,p,q,r,s,t,u,x,y,s,v,z,w,a,b,c,d,e,f,i,j,k,l#}::Indices(vector).
{m,n,p,q,r,s,t,u,x,y,s,v,z,w,a,b,c,d,e,f,i,j,k,l#}::Integer(0..D-1).
g_{i j}::Metric.
g^{i j}::InverseMetric.
g_{i j}::Symmetric.
S_{i j}::Symmetric.
W_{m n p q}::WeylTensor.
\delta{#}::KroneckerDelta;
R_{m n p q}::RiemannTensor;

ex5:=\delta^{i}_{a}^{j}_{b}^{k}_{c}^{l}_{d}^{m}_{e}^{n}_{r}^{t}_{p}^{x}_{q}^{y}_{z}^{s}_{w}
R^{a b}_{i j}R^{c d}_{k l}R^{e r}_{m n}R^{p q}_{t x}R^{z w}_{y s};
expand_delta(_);