H := binomial(3*n,3*k)^2 * binomial(3*n,3*k+1);

# Corresponding sequence:
SEQ := [seq(add(H, k=0..n), n=0..20)];

with(LREtools): _Env_LRE_tau := Sn; _Env_LRE_x := n;
with(SumTools:-Hypergeometric):
_MAXORDER := 20;

phi := proc(f) global n,k; subs(k=n-k, f) end:

# Projections to M+ and M-
phi_p := proc(f) (f + phi(f))/2 end:
phi_m := proc(f) (f - phi(f))/2 end:


tau := proc(f) global n,k; subs(k=k+1/3, f) end:

# Projections to M_{tau-1} and M_{tau^2+tau+1}
tau_1 := proc(f) (f  +tau(f)+tau(tau(f)))/3 end:
tau_2 := proc(f) (2*f-tau(f)-tau(tau(f)))/3 end:

t0 := time():
L1 := Zeilberger(tau_1(phi_p(H)), n, k, Sn)[1]: degree(%, Sn);
L2 := Zeilberger(tau_2(phi_p(H)), n, k, Sn)[1]: degree(%, Sn);
L3 := Zeilberger(tau_1(phi_m(H)), n, k, Sn)[1]: degree(%, Sn);
L4 := Zeilberger(tau_2(phi_m(H)), n, k, Sn)[1]: degree(%, Sn);

lprint("Time to compute telescoper in factored form", time() - t0);

Lmin := LCLM(L1, L2): degree(%, Sn);

# Test that the sequence indeed satisfies L_min:
seq(add(eval(coeff(Lmin,Sn,i), n=j) * SEQ[1+i+j], i=0..degree(Lmin,Sn)), j=0..10);

# The full telescoper:
L := LCLM(L1, L2, L3, L4): degree(%, Sn);

# Compare:
t0 := time():
LZ := Zeilberger(H, n, k, Sn)[1]: degree(%, Sn);
lprint("Time to compute telescoper with Zeilbergers algorithm", time() - t0);
normal(L - LZ);