This Mathematica notebook accompanies the paper by Christoph Koutschan and Thotsaporn Thanatipanonda with the title

A curious family of binomial determinants that count rhombus tilings of a holey hexagon

The following Mathematica packages are required. They are part of the RISCErgoSum bundle, which can be downloaded from
    http://www.risc.jku.at/research/combinat/software/ergosum/installation.html
For the download a password is required. It can be obtained by sending an e-mail to Peter Paule (ppaule@risc.uni-linz.ac.at).

<<RISC`HolonomicFunctions`;
<<RISC`Guess`;
<<RISC`LinearSystemSolver`;
SetDirectory[NotebookDirectory[]];

HolonomicFunctions Package version 1.7.3 (21-Mar-2017)
written by Christoph Koutschan
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Guess Package version 0.52
written by Manuel Kauers
Copyright Research Institute for Symbolic Computation (RISC),
Johannes Kepler University, Linz, Austria

Section 3: Related Determinants

Global Definitions

(* An auxiliary procedure to simplify certain expressions. *)
MySimp[expr_]:=expr/.
(* Simplify Floors *)
Floor[a_]:>Floor[Expand[a]]//.Floor[Optional[a_Integer]*(v:(i|j|k|n))+b_.]->a*v+Floor[b]/.
(* Simplify powers *)
(2^a_):>2^Expand[a]/.m1^a_:>Simplify[(-1)^a,Element[{i,j,k,n},Integers]]/.
(* Expand arguments of Pochhammers. *)
Pochhammer[a_,b_]:>Pochhammer@@Expand[{a,b}]//.
(* Simplify quotients of conjugate Pochhammers to rational functions. *)
Pochhammer[a1_,b1_]^c1_.*Pochhammer[a2_,b2_]^c2_./;IntegerQ[Expand[a1-a2]]&&IntegerQ[Expand[b1-b2]]&&c1>0&&c2<0:>With[{m=Min[c1,-c2]},FunctionExpand[(Pochhammer[a1,b1]/Pochhammer[a2,b2])^m]*Pochhammer[a1,b1]^(c1-m)*Pochhammer[a2,b2]^(c2+m)];

(* Definitions of matrix and determinant D_{s,t}(n) *)
DstMat[s_,t_,n_,mu_]:=Table[FunctionExpand[KroneckerDelta[i,j]+Binomial[mu+i+j-2,j]],{i,s,n+s-1},{j,t,n+t-1}];
dst[s_,t_,n_,i_,j_,mu_]:=FunctionExpand[KroneckerDelta[i+s,j+t]+Binomial[mu+i+s+j+t-4,j+t-1]];
dst[s_,t_,n_,i_,j_]:=dst[s,t,n,i,j,mu];
DstMat[s_,t_,n_,mu_]:=Table[dst[s,t,n,i,j,mu],{i,n},{j,n}];
DstMat[s_,t_,n_]:=DstMat[s,t,n,mu];
Dst[s_,t_,n_,mu_]:=Det[DstMat[s,t,n,mu]];
Dst[s_,t_,n_]:=Det[DstMat[s,t,n]];

(* Definitions of rational functions in lemmas and propositions of Section 2 *)
R00e[n_]:=Pochhammer[mu+2n,n]*Pochhammer[mu/2+2n+1/2,n-1]/Pochhammer[n,n]/Pochhammer[mu/2+n+1/2,n-1];
R00o[n_]:=Pochhammer[mu+2n-2,n-1]*Pochhammer[mu/2+2n-1/2,n]/Pochhammer[n,n]/Pochhammer[mu/2+n-1/2,n-1];
R00[n_]:=If[EvenQ[n],R00e[n/2],R00o[(n+1)/2]];
R10[n_]:=(1-2n)*Pochhammer[mu+2n,n]^2*Pochhammer[mu/2+2n+1/2,n-1]^2/(mu+2n)/Pochhammer[n,n]^2/Pochhammer[mu/2+n+1/2,n-1]^2;
R10[n_]:=-Pochhammer[mu+2n,n]*Pochhammer[mu+2n+1,n-1]*Pochhammer[mu/2+2n+1/2,n-1]^2/Pochhammer[n,n]/Pochhammer[n,n-1]/Pochhammer[mu/2+n+1/2,n-1]^2;
R01[n_]:=-Pochhammer[mu+2n-2,n+2]*Pochhammer[mu+2n+1,n-1]*Pochhammer[mu/2+2n+1/2,n-1]^2/Pochhammer[n,n-1]/Pochhammer[n,n+2]/Pochhammer[mu/2+n+1/2,n-1]^2;
R20[n_]:=Pochhammer[mu/2+2n+1/2,n-1]*Pochhammer[mu+2n+1,n-1]/Pochhammer[n,n-1]/Pochhammer[mu/2+n+1/2,n-1];
R02[n_]:=(2n-1)*Pochhammer[mu/2+2n+1/2,n-1]*Pochhammer[mu+2n-2,n+2]/(mu+2n)/Pochhammer[n,n+2]/Pochhammer[mu/2+n+1/2,n-1];

(* Tests *)
Table[Together[{
Dst[0,0,2n]/Dst[0,0,2n-1]-R00[2n-1],
Dst[1,0,2n+1]/Dst[1,0,2n-1]-R10[n],
Dst[0,1,2n+1]/Dst[0,1,2n-1]-R01[n],
Dst[2,0,2n]/Dst[2,0,2n-1]-R20[n],
Dst[0,2,2n]/Dst[0,2,2n-1]-R02[n],
Dst[-1,1,2n+1]/Dst[-1,1,2n]-R02[n]
}],{n,5}]
Table[Together[{
Dst[0,0,n]-2*Product[R00[i],{i,n-1}],
Dst[1,0,n]-If[EvenQ[n],0,Product[R10[i],{i,(n-1)/2}]],
Dst[0,1,n]-If[EvenQ[n],0,(mu-1)*Product[R01[i],{i,(n-1)/2}]]
}],{n,10}]

(* Alternative (old) definitions. The above formulas, however, are nicer. *)
R00a[n_]:=If[EvenQ[n],
2^(n/2)*Pochhammer[mu/2+n/2,Floor[(n+2)/4]]*Pochhammer[mu/2+n+1/2,(n-2)/2]/
Pochhammer[n/2,n/2]/Pochhammer[mu/2+Floor[3/4n]+1/2,Floor[(n-2)/4]],
2^((n-1)/2)*Pochhammer[mu/2+n/2-1/2,Floor[(n+1)/4]]*Pochhammer[mu/2+n+1/2,(n+1)/2]/
Pochhammer[(n+1)/2,(n+1)/2]/Pochhammer[mu/2+Floor[3/4(n-1)]+1/2,Floor[(n+1)/4]]];
R10a[n_]:=-(mu+2n)*Pochhammer[mu/2+2*n+1/2,n-1]^2*Pochhammer[mu/2+n+1,Floor[(n-1)/2]]^2/Pochhammer[mu/2+Floor[3/2n]+1/2,Floor[(n-1)/2]]^2/Pochhammer[n,n]/Pochhammer[1/2,n-1];
R01a[n_]:=(-2)*Pochhammer[mu+2*n-2,3]*Pochhammer[mu/2+2*n+1/2,n-1]^2*Pochhammer[mu/2+n+1,Floor[(n-1)/2]]^2/Pochhammer[mu/2+Floor[3/2n]+1/2,Floor[(n-1)/2]]^2/n/Pochhammer[n+2,n+1]/Pochhammer[1/2,n-1];
R20a[n_]:=Pochhammer[mu/2+2n+1/2,n-1]*Pochhammer[mu/2+n+1,Floor[(n-1)/2]]/2^(n-1)/Pochhammer[1/2,n-1]/Pochhammer[mu/2+Floor[3n/2]+1/2,Floor[(n-1)/2]];
R02a[n_]:=2^(n-1)*(2n-1)*Pochhammer[mu+2n-2,2]*Pochhammer[mu/2+2n+1/2,n-1]*Pochhammer[mu/2+n+1,Floor[(n-1)/2]]/Pochhammer[mu/2+Floor[3n/2]+1/2,Floor[(n-1)/2]]/Pochhammer[n,n+2];

(* Tests *)
Table[Together[{R00[n]/R00a[n],R10[n]/R10a[n],R01[n]/R01a[n],R20[n]/R20a[n],R02[n]/R02a[n]}],{n,20}]

Theorem 2:  (Subscript[D, 1,1])(2n)/(Subscript[D, 1,1])(2n-1)

The proof of this theorem is not part of the present paper. Here we verify it only on a few instances.

Table[Together[
(Dst[1,1,2n]/Dst[1,1,2n-1])/(Pochhammer[mu+2n,n]*Pochhammer[mu/2+2n+1/2,n-1]/Pochhammer[n,n]/Pochhammer[mu/2+n+1/2,n-1])],{n,5}]

Lemma 3:  (Subscript[D, 1,0])(2n)=0

We generate some data for the bivariate sequence (Subscript[c, n,j]). For 1≤n≤30 we compute the nullspace of the matrix of (Subscript[D, 1,0])(2n) using the command LinSolveUniv from the LinearSystemSolver package. For each n this nullspace is spanned by a single vector. We divide all entries of this vector by its last entry, therefore normalizing it such that the last entry is 1.We fill the obtained triangular array with zeros, since the guessing procedure requires a rectangular array as input.

Timing[data=PadRight[Table[ns=LinSolveUniv[DstMat[1,0,2n,mu],mu][[1]];Together[ns/ns[[-1]]],{n,30}]];]

We now use this data to guess recurrences. The settings (structure set and degree) have been found by trial and error. The Constraints option ensures that only the relevant data is used for the guessing, and not the padded zeros. The recurrences are converted into Ore polynomials (a data structure provided by the HolonomicFunctions package that represents recurrences as operators). The command NormalizeCoefficients removes any denominators and common factors among the coefficients of the recurrences. We verify that the guessed operators form a left Groebner basis, by feeding them into OreGroebnerBasis command and check that the output essentially equals the input (up to sign, this is achieved by GBEqual).

(* Directly guess a Groebner basis of annihilating operators for c[n,j]. *)
Timing[
ann=NormalizeCoefficients/@ToOrePolynomial[Join[
GuessMultRE[data,{c[n,j],c[n,j+1],c[n+1,j],c[n,j+2]},{n,j},7,StartPoint->{1,1},Constraints->j≤2n],
GuessMultRE[data,{c[n,j],c[n,j+1],c[n+1,j],c[n+1,j+1]},{n,j},8,StartPoint->{1,1},Constraints->j≤2n],
GuessMultRE[data,{c[n,j],c[n,j+1],c[n+1,j],c[n+2,j]},{n,j},17,StartPoint->{1,1},Constraints->j≤2n]],c[n,j]];
GBEqual[ann,OreGroebnerBasis[ann]]
]

Alternatively, we can guess recurrences of higher orders (and lower degree), which requires fewer data points and is therefore faster. The corresponding operators do not form a Groebner basis, so we have to call Buchberger’s algorithm (provided by the OreGroebnerBasis command) to obtain it. We verify that we get exactly the same result as before.

(* A faster alternative. Don't guess the Groebner basis directly, but some higher-order recurrences. *)
(* Then apply Buchberger. *)
Timing[
ann1=OreGroebnerBasis[NormalizeCoefficients/@ToOrePolynomial[
GuessMultRE[data,{c[n,j],c[n,j+1],c[n+1,j],c[n,j+2],c[n+1,j+1],c[n+2,j]},{n,j},7,StartPoint->{1,1},Constraints->j≤2n],c[n,j]]];
GBEqual[ann,ann1]
]

So far,  we have constructed a set ann of operators that generate an annihilating ideal in the Ore algebra defined by the shift operators (Subscript[S, n]) and (Subscript[S, j]). Together with a few initial values, these operators uniquely define a bivariate sequence (Subscript[c, n,j]). We now show that this sequence has the desired property, namely that the vector (Subscript[c, n,j]) with 1≤j≤2n is in the nullspace of the matrix, exhibiting that (Subscript[D, 1,0])(2n)=0. Note that for the proof of Lemma 3, it is not relevant how the set ann of operators was obtained. So in some sense, the proof of Lemma 3 starts here. We can load the precomputed set ann.

(*Put[ann,"ann_c_1_0.m"];*)
ann=Get["ann_c_1_0.m"];
UnderTheStaircase[ann]
ByteCount[ann]
Exponent[#,{n,j}]&/@ann

We first check whether the recurrence operators in ann allow to uniquely define a bivariate sequence. The command AnnihilatorSingularities gives all points (n,j) that cannot be computed using the recurrences (either because they lie under the staircase or because all applicable recurrence have a vanishing leading coefficient at the point). The command does not finish in reasonable time for symbolic μ, but for concrete μ it does. This is sufficient to conclude that there cannot be infinitely many singularities for symbolic μ.

AnnihilatorSingularities[ann,{0,0}]

AnnihilatorSingularities[ann/.mu->67,{0,0}]

(* Some tests *)
(* 1. Test that the annihilator is valid on the given data. *)
test=ApplyOreOperator[ann,c[n,j]];
Union[Flatten[Together[Table[test,{n,28},{j,28}]/.c[n_,j_]:>data[[n,j]]]]]
(* 2. Verify the identity Sum[d[2n,i,j]*c[n,j],{j,1,2n}]==0 for all 1≤i≤2n. *)
Union[Flatten[Table[Together[Sum[dst[1,0,2n,i,j,mu]*data[[n,j]],{j,1,2n}]],{n,30},{i,2n}]]]
(* 3. Same as 2., but with explicit form of matrix entries. This test can be viewed as initial value check. *)
Union[Flatten[Table[Together[Sum[FunctionExpand[Binomial[mu+i+j-3,j-1]]*data[[n,j]],{j,2n}]+data[[n,i+1]]],{n,29},{i,2n}]]]

Hence we want to prove the following identity, which shows that the vector (Subscript[c, n,j]), 1≤j≤2n, is in the nullspace of the matrix of the determinant (Subscript[D, 1,0])(2n).

TraditionalForm[HoldForm[Sum[Binomial[mu+i+j-3,j-1]*c[n,j],{j,1,2n}]==-c[n,i+1]]]

The sum has natural boundaries, since for j≤0 the binomial coefficient vanishes, and for j>2n we have (Subscript[c, n,j])=0. Thus creative telescoping yields an annihilator for the sum on the left hand side. Before calling the creative telescoping command, we have to compute the annihilator for the product of the binomial coefficient and (Subscript[c, n,j]), which is achieved by the command DFiniteTimes. It requires that both annihilators are given with respect to the same Ore algebra. Therefore, the Annihilator command for the binomial also includes the shift operator S[n], and to the annihilator ann, we have to add the shift operator S[i].

Timing[ct=FindCreativeTelescoping[DFiniteTimes[ToOrePolynomial[Append[ann,S[i]-1],OreAlgebra[S[n],S[j],S[i]]],Annihilator[Binomial[mu+i+j-3,j-1],{S[n],S[j],S[i]}]],S[j]-1];]

Since this computation takes some time, we can load the precomputed result.

ct=<<"ct_1_0.m";

It remains to verify that the sum satisfies the same system of recurrences as the (Subscript[c, n,i+1]) on the right hand side. Indeed, we find that the two annihilators agree. In addition, we have to compare a few initial values, which has already been done before. The proof is complete.

GBEqual[DFiniteSubstitute[ann,{j->i+1},Algebra->OreAlgebra[S[n],S[i]]],ct[[1]]]

Lemma 4:  (Subscript[D, 0,1])(2n)=0

The proof of Lemma 4 is very similar to the proof of Lemma 3. Therefore we give only few comments here and refer to the proof of Lemma 3 for more details.

(* Directly guess a Groebner basis of annihilating operators for c[n,j]. *)
Timing[
data=PadRight[Table[ns=LinSolveUniv[DstMat[0,1,2n,mu],mu][[1]];Together[ns/ns[[-1]]],{n,15}]];
ann=NormalizeCoefficients/@ToOrePolynomial[Join[
GuessMultRE[data,{c[n,j],c[n,j+1],c[n+1,j]},{n,j},9,StartPoint->{1,1},Constraints->j≤2n],
GuessMultRE[data,{c[n,j],c[n,j+1],c[n,j+2]},{n,j},6,StartPoint->{1,1},Constraints->j≤2n]
],c[n,j]];
GBEqual[ann,OreGroebnerBasis[ann]]
]

(*Put[ann,"ann_c_0_1.m"];*)
ann=Get["ann_c_0_1.m"];
UnderTheStaircase[ann]
ByteCount[ann]
Exponent[#,{n,j}]&/@ann

We check whether the recurrence operators in ann allow to uniquely define a bivariate sequence. We find that there are infinitely many singularities: the recurrences in ann are not applicable to compute the values (Subscript[c, n,2n]).

AnnihilatorSingularities[ann/.mu->67,{0,0}]

However, we can find another recurrence that is a consequence of the previous ones. This recurrence allows us to compute the values for (Subscript[c, n,2n]) (they are constant =1).

Factor[LeadingCoefficient[dfs=DFiniteSubstitute[ann,{j->2n},Algebra->OreAlgebra[S[n]]][[1]]]]

To see that the constant sequence is a solution of the above recurrence, we reduce its operator with the operator (Subscript[S, n])-1 and obtain 0. Hence (Subscript[S, n])-1 is a right factor of the recurrence dfs.

OreReduce[dfs,ToOrePolynomial[{S[n]-1}]]

(* Ee want to prove the following: *)
TraditionalForm[HoldForm[Sum[Binomial[mu+i+j-3,j]*c[n,j],{j,1,2n}]==-c[n,i-1]]]

(* Test initial values. *)
Union[Flatten[Table[Together[Sum[FunctionExpand[Binomial[mu+i+j-3,j]]*data[[n,j]],{j,2n}]+If[i===1,0,data[[n,i-1]]]],{n,14},{i,2n}]]]

(* The sum has natural boundaries. Thus creative telescoping gives its annihilator. *)
Timing[ct=FindCreativeTelescoping[DFiniteTimes[ToOrePolynomial[Append[ann,S[i]-1],OreAlgebra[S[n],S[j],S[i]]],Annihilator[Binomial[mu+i+j-3,j],{S[n],S[j],S[i]}]],S[j]-1];
]

ct=<<"ct_0_1.m";

GBEqual[DFiniteSubstitute[ann,{j->i-1},Algebra->OreAlgebra[S[n],S[i]]],ct[[1]]]

Lemma 5:  (Subscript[D, 0,0])(2n)/(Subscript[D, 0,0])(2n-1)

The determinant evaluation of (Subscript[D, 0,0])(n) was first proven by George E. Andrews. We include our computer proof just for sake of illustration. It follows the very same strategy as the proof of Lemma 6 (where also more detailed explanations can be found, see below).

(* Directly guess a Groebner basis of annihilating operators for c[n,j]. *)
Timing[
data=PadRight[Table[ns=LinSolveUniv[Most[DstMat[0,0,2n,mu]],mu][[1]];Together[ns/ns[[-1]]],{n,15}]];
ann=NormalizeCoefficients/@ToOrePolynomial[Join[
GuessMultRE[data,{c[n,j],c[n,j+1],c[n+1,j]},{n,j},8,StartPoint->{1,1},Constraints->j≤2n],
GuessMultRE[data,{c[n,j],c[n,j+1],c[n,j+2]},{n,j},4,StartPoint->{1,1},Constraints->j≤2n]
],c[n,j]];
GBEqual[ann,OreGroebnerBasis[ann]]
]

(*Put[ann,"ann_c_0_0.m"];*)
ann=Get["ann_c_0_0.m"];
UnderTheStaircase[ann]
ByteCount[ann]
Exponent[#,{n,j}]&/@ann

(* Identity (1): Show that c[n,2n]=1 for all n. *)
diag=DFiniteSubstitute[ann,{j->2n},Algebra->OreAlgebra[S[n]]][[1]];
Factor[LeadingCoefficient[diag]]
OreReduce[diag,{ToOrePolynomial[S[n]-1]}]

(* Identity (2) *)
TraditionalForm[HoldForm[Sum[Binomial[μ+i+j-4,j-1]*c[n,j],{j,1,2n}]==-c[n,i]"    "(1≤i≤2n-1)]]

(* Identity (2): numerical check (= initial values) *)
Union[Flatten[Table[Together[Sum[dst[0,0,2n,i,j,mu]*data[[n,j]],{j,1,2n}]],{n,15},{i,2n-1}]]]
Union[Flatten[Table[Together[Sum[FunctionExpand[Binomial[mu+i+j-4,j-1]]*data[[n,j]],{j,1,2n}]+data[[n,i]]],{n,15},{i,2n-1}]]]

(* The sum has natural boundaries. Thus creative telescoping gives its annihilator. *)
Timing[ct=FindCreativeTelescoping[DFiniteTimes[ToOrePolynomial[Append[ann,S[i]-1],OreAlgebra[S[n],S[j],S[i]]],Annihilator[Binomial[mu+i+j-4,j-1],{S[n],S[j],S[i]}]],S[j]-1];
]

ct=<<"ct_0_0_A.m";

GBEqual[DFiniteSubstitute[ann,{j->i},Algebra->OreAlgebra[S[n],S[i]]],ct[[1]]]

(* Identity (3) *)
TraditionalForm[HoldForm[Sum[(KroneckerDelta[2n,j]+Binomial[μ+2n+j-4,j-1])*c[n,j],{j,1,2n}]==(Subscript[D, 0,0])[2n]/(Subscript[D, 0,0])[2n-1]]]

(* Numerical check of Identity (3) (= initial value check) *)
Table[Together[Sum[FunctionExpand[Binomial[mu+2n+j-4,j-1]]*data[[n,j]],{j,1,2n}]+data[[n,2n]]-R00[2n-1]],{n,10}]

(* The sum has natural boundaries. Thus creative telescoping gives its annihilator. *)
Timing[
ct=FindCreativeTelescoping[DFiniteTimes[ann,Annihilator[Binomial[mu+2n+j-4,j-1],{S[n],S[j]}]],S[j]-1];
]

ct=<<"ct_0_0_B.m";

(* Add the part coming from the KroneckerDelta. *)
ann3=DFinitePlus[ct[[1]],DFiniteSubstitute[ann,{j->2n}]];

(* The recurrence we get is not too unhandy... *)
Factor[ann3]

(* Plug our expression into the recurrence *)
test=ApplyOreOperator[ann3[[1]],R00o[n]];
(* Divide by the (non-zero!) expression itself, to make simplification easier. *)
test=If[Head[test]===Plus,#/R00o[n]&/@test,test/R00o[n]];
(* Simplify *)
Together[MySimp[test]]

Lemma 6:  (Subscript[D, 2,0])(2n)/(Subscript[D, 2,0])(2n-1)

For 1≤n≤15  we compute the nullspace of the (2n-1)×(2n) matrix that is obtained by removing the last row from the matrix of (Subscript[D, 2,0])(2n). In all cases we encounter a one-dimensional nullspace, hence it is spanned by a single vector. We normalize this vector by dividing it by its last entry. For each n we obtain a vector of length 2n. The triangular array formed by these vectors is filled with 0’s such that we get a rectangular array.

Timing[
data=PadRight[Table[ns=LinSolveUniv[Most[DstMat[2,0,2n,mu]],mu][[1]];Together[ns/ns[[-1]]],{n,15}]];
]

As described in the proof of Lemma 3, we use the GuessMultRE command to construct recurrence equations that are satisfied by the computed (finite amount of) data. The recurrences are converted into operators and a left Groebner basis of the left ideal generated by these operators is computed.

(* Don't guess the Groebner basis directly, but some higher-order recurrences. Then apply Buchberger. *)
Timing[
ann=OreGroebnerBasis[NormalizeCoefficients/@ToOrePolynomial[
GuessMultRE[data,{c[n,j],c[n,j+1],c[n+1,j],c[n,j+2],c[n+1,j+1]},{n,j},7,StartPoint->{1,1},Constraints->j≤2n],c[n,j]]];
]

In the following we prove certain properties and identities concerning the infinite bivariate sequence (Subscript[c, n,j]) defined by the recurrences (ann) and initial values (data).

(*Put[ann,"ann_c_2_0.m"];*)
ann=Get["ann_c_2_0.m"];
UnderTheStaircase[ann]
ByteCount[ann]
Exponent[#,{n,j}]&/@ann

The first identity is (Subscript[c, n,2n])=1. From the annihilator of (Subscript[c, n,j]) we compute a recurrence operator for the diagonal sequence (Subscript[c, n,2n]). In order to study the solutions of this univariate recurrence, we have to look at its singularities (integer roots of its leading coefficient). By listing all linear factors of degree 1 in n, we see that there are no positive integer roots (μ is considered a symbolic parameter, and hence we need not care about special values of μ). When we reduce the diag operator with the operator (Subscript[S, n])-1, we obtain 0, which tells us that (Subscript[S, n])-1 is a right factor, and therefore that any constant sequence is a solution of diag. Since the initial values for (Subscript[c, n,2n]) are 1 (by construction), we conclude that (Subscript[c, n,2n])=1 for all n.

(* Identity (1): Show that c[n,2n]=1 for all n. *)
diag=DFiniteSubstitute[ann,{j->2n},Algebra->OreAlgebra[S[n]]][[1]];
Select[Factor[LeadingCoefficient[diag]],Exponent[#,n]===1||Head[#]=!=Plus&]
OreReduce[diag,{ToOrePolynomial[S[n]-1]}]

The second identity that we want to prove is displayed below. We first perform a numerical check with the available data. This check also serves as comparison of initial values that we will have to do at the end.

(* Identity (2) *)
TraditionalForm[HoldForm[Sum[Binomial[μ+i+j-2,j-1]*c[n,j],{j,1,2n}]==-c[n,i+2]"    "(1≤i≤2n-1)]]

(* Identity (2): numerical check (= initial values) *)
Union[Flatten[Table[Together[Sum[dst[2,0,2n,i,j,mu]*data[[n,j]],{j,1,2n}]],{n,15},{i,2n-1}]]]
Union[Flatten[Table[Together[Sum[FunctionExpand[Binomial[mu+i+j-2,j-1]]*data[[n,j]],{j,1,2n}]+data[[n,i+2]]],{n,14},{i,2n-1}]]]

The sum has natural boundaries (meaning that the summand is zero outside the given summation range), since for j<1 the binomial coefficient is 0 and for j>2n we have (Subscript[c, n,j])=0. Thus creative telescoping gives the annihilator of the sum (ct[[1]]), provided that the certificate (ct[[2]]) does not have poles inside the summation range.

(* The sum has natural boundaries. Thus creative telescoping gives its annihilator. *)
Timing[ct=FindCreativeTelescoping[DFiniteTimes[ToOrePolynomial[Append[ann,S[i]-1],OreAlgebra[S[n],S[j],S[i]]],Annihilator[Binomial[mu+i+j-2,j-1],{S[n],S[j],S[i]}]],S[j]-1];
]

(* Load the precomputed expression *)
ct=<<"ct_2_0_A.m";

(* Look at the singularities of the certificate: no pole for 1 ≤ j ≤ 2n. *)
Factor[PolynomialLCM@@Denominator[Together[Flatten[OrePolynomialListCoefficients/@Flatten[ct[[2]]]]]]]

We see that the recurrences satisfied by the sum and constructed by creative telescoping are identical with the recurrences satisfied by (Subscript[c, n,i+2]). Hence by comparing initial values (already done), we conclude that both expressions agree, thereby establishing identity (2).

GBEqual[DFiniteSubstitute[ann,{j->i+2},Algebra->OreAlgebra[S[n],S[i]]],ct[[1]]]

Identitites (1) and (2) together imply that (Subscript[c, n,j]) is the cofactor of the Laplace expansion of the matrix of (Subscript[D, 2,0])(2n) with respect to the last row, divided by (Subscript[D, 2,0])(2n-1). Hence, by proving the following identity (3), we establish the statement of Lemma 6.

(* Identity (3) *)
TraditionalForm[HoldForm[Sum[(KroneckerDelta[2n+1,j-1]+Binomial[μ+2n+j-2,j-1])*c[n,j],{j,1,2n}]==(Subscript[D, 2,0])[2n]/(Subscript[D, 2,0])[2n-1]==r20]/.r20->R20[n]/.mu->μ]

(* Numerical check of Identity (3) (= initial value check) *)
Table[Together[Sum[FunctionExpand[Binomial[mu+2n+j-2,j-1]]*data[[n,j]],{j,1,2n}]-R20[n]],{n,10}]

As before, we find that the sum in identity (3) has natural boundaries. Note that the Kronecker delta can be omitted since it does not contribute to the sum. Using the DFiniteTimes command, we construct an annihilator for the product of the binomial times (Subscript[c, n,j]).

(* The sum has natural boundaries. Thus creative telescoping gives its annihilator. *)
Timing[
ct=CreativeTelescoping[DFiniteTimes[ann,Annihilator[Binomial[mu+2n+j-2,j-1],{S[n],S[j]}]],S[j]-1];
]

ct=<<"ct_2_0_B.m";

(* Look at the singularities of the certificate: no pole for 1 ≤ j ≤ 2n. *)
Factor[PolynomialLCM@@Denominator[Together[OrePolynomialListCoefficients[Flatten[ct][[2]]]]]]

As a result, we get a nice recurrence (in n) for the sum on the left-hand side of identity (3):

Factor[ct[[1]]]

It remains to show that the expression (Subscript[R, 2,0])(n) also satisfies this recurrence (initial values have already been checked). This can, for example, be done by plugging (Subscript[R, 2,0])(n) into the recurrence and by simplifying:

(* Plug our expression R_{2,0} into the recurrence *)
test=ApplyOreOperator[ct[[1,1]],R20[n]];
(* Divide by the (non-zero!) expression itself, to make simplification easier. *)
test=If[Head[test]===Plus,#/R20[n]&/@test,test/R20[n]];
(* Simplify *)
Together[MySimp[test]]

Alternatively, we can compute an annihilator for (Subscript[R, 2,0])(n) and show that it is a right factor of the recurrence operator ct[[1,1]]; the latter is done via the OreReduce command that gives a zero remainder.

OreReduce[ct[[1,1]],Annihilator[R20[n],S[n]]]

Lemma 7:  (Subscript[D, 0,2])(2n)/(Subscript[D, 0,2])(2n-1)

The proof of Lemma 7 follows the very same strategy as the proof of Lemma 6. There, detailed explanations have been given that apply equally to the following calculations. Hence, we list only the Mathematica commands for the computations, and refer to Lemma 6 for the explanations.

Timing[
data=PadRight[Table[ns=LinSolveUniv[Most[DstMat[0,2,2n,mu]],mu][[1]];Together[ns/ns[[-1]]],{n,23}]];
]

Timing[
rec1=First[GuessMultRE[data,{c[n,j],c[n,j+1],c[n+1,j]},{n,j},{10,11},StartPoint->{1,1},Constraints->j≤2n&&n≥1,AdditionalEquations->Infinity]];
]

Timing[
rec2=First[GuessMultRE[data,{c[n,j],c[n,j+1],c[n,j+2]},{n,j},12,StartPoint->{1,1},Constraints->j≤2n&&n≥1,AdditionalEquations->Infinity]];
]

(* CAVEAT: these recurrences do not hold for n=j=1. *)
TableForm[Map[If[#=!={0,0},"*",0]&,Together[Table[{rec1,rec2},{n,10},{j,2n}]/.c[n_,j_]:>data[[n,j]]],{2}]]

Together[{rec1,rec2}/.{j->1,n->1}/.c[n_,j_]:>data[[n,j]]]

ann=OreGroebnerBasis[NormalizeCoefficients/@ToOrePolynomial[{rec1,rec2},c[n,j]]];
Support[ann]

(* Put[ann,"ann_c_0_2.m"]; *)
ann=Get["ann_c_0_2.m"];
UnderTheStaircase[ann]
ByteCount[ann]
Exponent[#,{n,j}]&/@ann

(* Identity (1): Show that c[n,2n]=1 for all n. *)
diag=DFiniteSubstitute[ann,{j->2n},Algebra->OreAlgebra[S[n]]][[1]];
Select[Factor[LeadingCoefficient[diag]],Exponent[#,n]===1||Head[#]=!=Plus&]
OreReduce[diag,{ToOrePolynomial[S[n]-1]}]

(* Identity (2) *)
TraditionalForm[HoldForm[Sum[Binomial[μ+i+j-2,j+1]*c[n,j],{j,1,2n}]==-c[n,i-2]"    "(1≤i≤2n-1)]]

(* Identity (2): numerical check (= initial values). Also the special case n=1 is covered. *)
Union[Flatten[Table[Together[Sum[dst[0,2,2n,i,j,mu]*data[[n,j]],{j,1,2n}]],{n,15},{i,2n-1}]]]
Union[Flatten[Table[Together[Sum[FunctionExpand[Binomial[mu+i+j-2,j+1]]*data[[n,j]],{j,1,2n}]+If[i<3,0,data[[n,i-2]]]],{n,14},{i,2n-1}]]]

(* The sum has natural boundaries. Thus creative telescoping gives its annihilator. *)
Timing[ct=FindCreativeTelescoping[DFiniteTimes[ToOrePolynomial[Append[ann,S[i]-1],OreAlgebra[S[n],S[j],S[i]]],Annihilator[Binomial[mu+i+j-2,j+1],{S[n],S[j],S[i]}]],S[j]-1];
]

ct=<<"ct_0_2_A.m";

(* Look at the singularities of the certificate: no pole for 1 ≤ j ≤ 2n. *)
Factor[PolynomialLCM@@Denominator[Together[Flatten[OrePolynomialListCoefficients/@Flatten[ct[[2]]]]]]]

GBEqual[DFiniteSubstitute[ann,{j->i-2},Algebra->OreAlgebra[S[n],S[i]]],ct[[1]]]

(* Identity (3) *)
TraditionalForm[HoldForm[Sum[(KroneckerDelta[2n-1,j+1]+Binomial[μ+2n+j-2,j+1])*c[n,j],{j,1,2n}]==(Subscript[D, 0,2])[2n]/(Subscript[D, 0,2])[2n-1]==r02]/.r02->R02[n]/.mu->μ]

(* Numerical check of Identity (3) (= initial value check) *)
Table[Together[Sum[FunctionExpand[Binomial[mu+2n+j-2,j+1]]*data[[n,j]],{j,1,2n}]+If[n>1,data[[n,2n-2]],0]-R02[n]],{n,10}]

(* The sum has natural boundaries. Thus creative telescoping gives its annihilator. *)
Timing[
ct=FindCreativeTelescoping[DFiniteTimes[ann,Annihilator[Binomial[mu+2n+j-2,j+1],{S[n],S[j]}]],S[j]-1];
]

(*Put[ct,"ct_0_2_B.m"];*)
ct=Get["ct_0_2_B.m"];

(* Look at the singularities of the certificate: no pole for 1 ≤ j ≤ 2n. *)
Factor[PolynomialLCM@@Denominator[Together[OrePolynomialListCoefficients[Flatten[ct][[2]]]]]]

In contrast to Lemma 6, here the Kronecker delta in identity (3) plays a role (it gives 1 for j=2n-2). Hence we have to add the corresponding term (Subscript[c, n,2n-2]). The command DFinitePlus yields a recurrence operator (of order 2) for the sum (without Kronecker delta) plus (Subscript[c, n,2n-2]).

annSum=DFinitePlus[ct[[1]],DFiniteSubstitute[ann,{j->2n-2}]];
Support[annSum]

(* Plug our expression into the recurrence. *)
test=ApplyOreOperator[annSum[[1]],R02[n]];
(* Divide by the (non-zero!) expression itself, to make simplification easier. *)
test=If[Head[test]===Plus,#/R02[n]&/@test,test/R02[n]];
(* Simplify *)
Together[MySimp[test]]

(* Alternative: reduce annSum with the annihilator of R02 *)
OreReduce[annSum[[1]],Annihilator[R02[n],S[n]]]

Section 4: Nice Formula for (Subscript[D, 1,1])(n)

Simple formula for (∏_(j=1)^(k-1))(((Subscript[R, 1,0])(j) (Subscript[R, 0,1])(j)) / ((Subscript[R, 0,0])(2j-1) (Subscript[R, 0,0])(2j)))

(* The factor inside the product *)
fac=MySimp[R10[j]*R01[j]/R00o[j]/R00e[j]]/.Pochhammer[j,j+2]->Pochhammer[j,j-1]*(2j-1)*(2j)*(2j+1)

(* Test *)
Table[Together[(R10[j]*R01[j]/R00[2j-1]/R00[2j])/fac],{j,10}]

(* Final form for the product (use prod to avoid evaluation) *)
PR0110=(Pochhammer[mu,3k-3]/Pochhammer[mu/2+k-1/2,k-1]/(2k-1)!)*
prod[Pochhammer[mu/2+2j+1/2,j-1]*Pochhammer[mu+2j+1,j-1]/Pochhammer[j,j-1]/Pochhammer[mu/2+j+1/2,j-1],{j,1,k-1}]^2;
TraditionalForm[(HoldForm@@{PR0110})/.{prod->Product,mu->μ}]

(* Test *)
Table[Together[Product[R10[j]*R01[j]/R00[2j-1]/R00[2j],{j,1,k-1}]/(PR0110/.prod->Product)],{k,10}]

Simple formula for (∏_(i=2k)^n)(Subscript[R, 0,0])(i)

(* Split the product into even and odd instances of R00 *)
TraditionalForm[Product[R[i],{i,2k,n}]==Product[R[2j],{j,k,Floor[n/2]}]*Product[R[2j+1],{j,k,Floor[(n-1)/2]}]]

(* Derivation *)
PR00=prod[R00e[j],{j,k,Floor[n/2]}]*prod[R00o[j+1],{j,k,Floor[(n-1)/2]}];
PR00=PR00/.Pochhammer[a_,b_]:>Pochhammer[Expand[a],Expand[b]];
PR00e=PR00/.Floor[a_]:>FullSimplify[Floor[a],Element[n/2,Integers]];
PR00o=PR00/.{Floor[n/2]->(n-1)/2,Floor[(n-1)/2]->(n-1)/2}/.
prod[a_,c_]*prod[b_,c_]->prod[a*b,c]/.
Pochhammer[mu/2+2j+1/2,j-1]->Pochhammer[mu/2+2j-1/2,j]/(mu/2+2j-1/2)/.
Pochhammer[mu/2+j+1/2,j-1]->Pochhammer[mu/2+j+1/2,j]/(mu/2+2j-1/2);
{PR00e,PR00o}={PR00e,PR00o}/.prod[a_,b_]:>prod[a,Expand[b]]

(* Tests *)
Flatten[Table[Together[(PR00o/.prod->Product)/Product[R00[i],{i,2k,n}]],{n,1,9,2},{k,0,Floor[(n+1)/2]}]]
Flatten[Table[Together[(PR00e/.prod->Product)/Product[R00[i],{i,2k,n}]],{n,2,10,2},{k,0,Floor[(n+1)/2]}]]

Final formula for (Subscript[D, 1,1])(n)

(* The case distinction (k=0 vs. k>0) *)
Together[PR0110*(2k-1)!/.k->0/.prod->Product]

(* n even *)
FFe=sum[(mu-1)/2*If[k==0,4(mu-2),1/(2k-1)!]*PR00e*(PR0110*(2k-1)!),{k,0,n/2}];
(* n odd *)
FFo=sum[(mu-1)/2*If[k==0,4(mu-2),1/(2k-1)!]*PR00o*(PR0110*(2k-1)!),{k,0,(n+1)/2}];
{FFe,FFo}={FFe,FFo}/.(mu-1)*Pochhammer[mu,3k-3]->Pochhammer[mu-1,3k-2];

TraditionalForm[HoldForm[FF]/.FF->FFe/.{prod->Product,sum->Sum,mu->μ}]

TraditionalForm[HoldForm[FF]/.FF->FFo/.{prod->Product,sum->Sum,mu->μ}]

(* Test *)
Table[Together[(If[EvenQ[n],FFe,FFo]/.sum->Sum/.prod->Product)/Dst[1,1,n]],{n,10}]

Section 5: Proof of the Monstrous Conjecture

The Conjecture

Clear[myC,myE,myF,myT,myS1,myS2,myP1,myP2,myG,D34];
myC[n_Integer?Positive]:=((-1)^n+3)/2*Product[Floor[i/2]!/i!,{i,1,n}];
myE[n_Integer?Positive,mu_]:=Pochhammer[mu+1,n]*Product[(mu+2i+6)^(2*Floor[(i+2)/3]),{i,1,Floor[3/2*Floor[(n-1)/2]-2]}]*
Product[(mu+2i+2*Floor[3/2*Floor[n/2+1]]-1)^(2*Floor[Floor[n/2]/2-(i-1)/3]-1),{i,1,Floor[3/2*Floor[n/2]-2]}];
myF[m:(0|1),n_Integer?Positive,mu_]:=
Product[(mu+2i+n+m)^(1-2i-m),{i,1,Floor[(n-1)/4]}]*Product[(mu-2i+2n-2m+1)^(1-2i-m),{i,1,Floor[n/4-1]}];
myF[n_Integer?Positive,mu_]:=Which[
EvenQ[n],myE[n,mu]*myF[0,n,mu],
OddQ[n],myE[n,mu]*myF[1,n,mu]*Product[mu+2n+2i-1,{i,1,(n-5)/2}]];
myT[k_,mu_]:=-12+84*k+288*k^2-5856*k^3+20352*k^4-41472*k^5+55296*k^6+10*mu+76*k*mu-2176*k^2*mu+9888*k^3*mu-25344*k^4*mu+41472*k^5*mu+10*mu^2-261*k*mu^2+1676*k^2*mu^2-5472*k^3*mu^2+11520*k^4*mu^2-10*mu^3+115*k*mu^3-488*k^2*mu^3+1440*k^3*mu^3+2*mu^4-15*k*mu^4+76*k^2*mu^4+k*mu^5;
myS1[n_Integer?Positive,mu_]:=Sum[myT[k,mu]*2^(6k)*(mu+8 k-1)*Pochhammer[1/2,2k-1]^2*Pochhammer[(mu+4k+2)/2,2n-2k-2]*Pochhammer[(mu+5)/2,2k-3]*Pochhammer[(mu+4k+2)/2,k-2]/
((2k)!*Pochhammer[(mu+6k-3)/2,3k+4]),{k,1,n-1}];
myS2[n_Integer?Positive,mu_]:=
Sum[myT[k+1/2,mu]*2^(6k)*(mu+8 k+3)*Pochhammer[1/2,2k]^2*Pochhammer[(mu+4k+4)/2,2n-2k-2]*Pochhammer[(mu+5)/2,2k-2]*Pochhammer[(mu+4k+4)/2,k-2]/
((2k+1)!*Pochhammer[(mu+6k+1)/2,3k+5]),{k,1,n-1}];
myP1[n_Integer?Positive,mu_]:=Together[
2^(3n-1)*Pochhammer[(mu+6n-3)/2,3n-2]/Pochhammer[(mu+5)/2,2n-3]*
(2^(-13)*mu*(mu-1)*myS1[n,mu]+1/(mu+3)^2*Pochhammer[(mu+2)/2,2n-2])];
myP2[n_Integer?Positive,mu_]:=Together[
2^(3n-1)*Pochhammer[(mu+6n+1)/2,3n-1]/Pochhammer[(mu+5)/2,2n-2]*
(2^(-9)mu*(mu-1)*myS2[n,mu]+(mu+14)/((mu+7)*(mu+9))*Pochhammer[(mu+4)/2,2n-2])];
myG[n_Integer?Positive,mu_]:=If[EvenQ[n],myP2[n/2,mu],myP1[(n+1)/2,mu]];
D34[n_Integer?Positive,mu_]:=myC[n]*myF[n,mu]*myG[Floor[(n+1)/2],mu];

(* even n *)
Table[Together[(FFe/.sum->Sum/.prod->Product)/(myC[n]*myE[n,mu]*myF[0,n,mu]*myP2[n/4,mu])],{n,4,12,4}]
Table[Together[(FFe/.sum->Sum/.prod->Product)/(myC[n]*myE[n,mu]*myF[0,n,mu]*myP1[(n+2)/4,mu])],{n,2,10,4}]
(* odd n *)
Table[Together[(FFo/.sum->Sum/.prod->Product)/(myC[n]*myE[n,mu]*myF[1,n,mu]*Product[mu+2n+2i-1,{i,1,(n-5)/2}]*myP1[(n+3)/4,mu])],{n,1,9,4}]
Table[Together[(FFo/.sum->Sum/.prod->Product)/(myC[n]*myE[n,mu]*myF[1,n,mu]*Product[mu+2n+2i-1,{i,1,(n-5)/2}]*myP2[(n+1)/4,mu])],{n,3,11,4}]

(* These are the recurrence operators that we used to define P_1 and P_2. *)
rec1=NormalizeCoefficients[ToOrePolynomial[-myP1fast[n,mu]+(-8*(-10+mu+4*n)*(-8+mu+4*n)*(-14+mu+6*n)*(-12+mu+6*n)*(-10+mu+6*n)*(-9+mu+8*n)*(-31+mu+12*n)*(-29+mu+12*n)*(-27+mu+12*n)*(-25+mu+12*n)*(-23+mu+12*n)*(-21+mu+12*n)*(81-72*n+16*n^2)*(49-56*n+16*n^2)*(123168-78946*mu+18939*mu^2-2053*mu^3+93*mu^4-mu^5+(-9+mu)*(70956-30208*mu+3989*mu^2-158*mu^3+mu^4)*n+4*(346032-149656*mu+21803*mu^2-1202*mu^3+19*mu^4)*n^2+96*(-9+mu)*(1861-402*mu+15*mu^2)*n^3+384*(2753-606*mu+30*mu^2)*n^4+41472*(-9+mu)*n^5+55296*n^6)*myP1fast[-2+n,mu]+(-13+mu+8*n)*(-6*(-15+mu)*(-11+mu)*(-9+mu)*(77272834343040-90508623095808*mu+46786094223720*mu^2-14041912717156*mu^3+2707887452266*mu^4-350541498059*mu^5+30888280625*mu^6-1838952303*mu^7+72032193*mu^8-1778033*mu^9+26555*mu^10-241*mu^11+mu^12)+(-13+mu)*(615591764176296960-787691318438414592*mu+453271146257615040*mu^2-154970921382725880*mu^3+35030740197791460*mu^4-5511255715119386*mu^5+618465455797003*mu^6-49890145667170*mu^7+2877469024970*mu^8-116576723262*mu^9+3218550024*mu^10-58094110*mu^11+655730*mu^12-4400*mu^13+13*mu^14)*n+(43790163197061415680-55769554581921674496*mu+32100807569482408752*mu^2-11046065343390418896*mu^3+2532539665806086200*mu^4-408068212472225048*mu^5+47486735062736003*mu^6-4036853597489641*mu^7+250606824181572*mu^8-11237476473228*mu^9+356071800098*mu^10-7704642502*mu^11+108621484*mu^12-941780*mu^13+4611*mu^14-9*mu^15)*n^2+2*(-13+mu)*(5765368315087296000-6423796647403130880*mu+3186986194272026736*mu^2-928737086880929008*mu^3+176577512806080224*mu^4-23002876518214396*mu^5+2097912117891133*mu^6-134465197774532*mu^7+5992468266728*mu^8-181075265324*mu^9+3560096842*mu^10-42928700*mu^11+293696*mu^12-1000*mu^13+mu^14)*n^3+8*(44967647815472773440-49875119477893931904*mu+24771543294236452512*mu^2-7277588373063623552*mu^3+1407087781066080464*mu^4-188436568279081716*mu^5+17910169812661579*mu^6-1217322600443922*mu^7+58827888448174*mu^8-1983671151898*mu^9+45113742796*mu^10-655655046*mu^11+5605730*mu^12-24666*mu^13+41*mu^14)*n^4+32*(-13+mu)*(1545137447830050528-1468846207754989056*mu+613359955784013384*mu^2-148046294338567160*mu^3+22867645137091796*mu^4-2363768523778396*mu^5+166104951524749*mu^6-7900529853234*mu^7+248588564859*mu^8-4947975304*mu^9+57722923*mu^10-345266*mu^11+785*mu^12)*n^5+128*(6923436910786740816-6551979917272781760*mu+2741775205145125620*mu^2-668624737408815316*mu^3+105402483452844020*mu^4-11258804752461004*mu^5+830334150499955*mu^6-42256983681030*mu^7+1457399275653*mu^8-32763679904*mu^9+447520681*mu^10-3258554*mu^11+9319*mu^12)*n^6+1024*(-13+mu)*(72414477952775604-57105723925009800*mu+19399742350341207*mu^2-3719307354992416*mu^3+442850412559382*mu^4-33955375237500*mu^5+1681820711178*mu^6-52507834704*mu^7+974233650*mu^8-9518828*mu^9+36355*mu^10)*n^7+4096*(204759442490425380-160746724570083012*mu+54801297077548677*mu^2-10648677530738482*mu^3+1300829127395384*mu^4-103865351431818*mu^5+5455145057379*mu^6-184594947228*mu^7+3811103508*mu^8-42749540*mu^9+194248*mu^10)*n^8+49152*(-13+mu)*(920215916156142-577914239846832*mu+151701784373213*mu^2-21614250577806*mu^3+1815722558519*mu^4-91353917016*mu^5+2663224490*mu^6-40669644*mu^7+245586*mu^8)*n^9+196608*(1693595159851230-1058822980698432*mu+279542833819585*mu^2-40572445515984*mu^3+3526446267001*mu^4-187021320840*mu^5+5872755784*mu^6-99020958*mu^7+679074*mu^8)*n^10+21233664*(-13+mu)*(550446775412-258091315032*mu+47985773125*mu^2-4496668860*mu^3+222288724*mu^4-5456352*mu^5+51547*mu^6)*n^11+28311552*(1958821138060-914306594496*mu+171668385371*mu^2-16540689390*mu^3+859090262*mu^4-22689546*mu^5+236549*mu^6)*n^12+21403533312*(-13+mu)*(57395792-17859456*mu+1964631*mu^2-89610*mu^3+1425*mu^4)*n^13+12230590464*(290157464-89880912*mu+10081119*mu^2-483594*mu^3+8337*mu^4)*n^14+3522410053632*(-13+mu)*(12823-1986*mu+69*mu^2)*n^15+3522410053632*(19340-2982*mu+111*mu^2)*n^16+380420285792256*(-13+mu)*n^17+169075682574336*n^18)*myP1fast[-1+n,mu])/((-1+n)*(-3+2*n)*(-6+mu+4*n)*(-5+mu+4*n)*(-4+mu+4*n)*(-3+mu+4*n)*(-9+mu+6*n)*(-7+mu+6*n)*(-5+mu+6*n)*(-17+mu+8*n)*(-2*(-2619750+910279*mu-117666*mu^2+6856*mu^3-168*mu^4+mu^5)+(-17+mu)*(862188-199648*mu+14213*mu^2-302*mu^3+mu^4)*n+4*(4278168-996880*mu+77747*mu^2-2282*mu^3+19*mu^4)*n^2+96*(-17+mu)*(6541-762*mu+15*mu^2)*n^3+384*(9773-1146*mu+30*mu^2)*n^4+41472*(-17+mu)*n^5+55296*n^6)),myP1fast[n,mu],OreAlgebra[S[n]]]];
rec2=NormalizeCoefficients[ToOrePolynomial[-myP2fast[n,mu]+(-8*(-8+mu+4*n)*(-6+mu+4*n)*(-12+mu+6*n)*(-10+mu+6*n)*(-8+mu+6*n)*(-5+mu+8*n)*(-25+mu+12*n)*(-23+mu+12*n)*(-21+mu+12*n)*(-19+mu+12*n)*(-17+mu+12*n)*(-15+mu+12*n)*(49-56*n+16*n^2)*(25-40*n+16*n^2)*(-((-3+mu)*(2788-2196*mu+577*mu^2-54*mu^3+mu^4))+2*(-5+mu)*(7620-5536*mu+1253*mu^2-86*mu^3+mu^4)*n+8*(35820-26716*mu+6791*mu^2-662*mu^3+19*mu^4)*n^2+192*(-5+mu)*(601-222*mu+15*mu^2)*n^3+768*(863-336*mu+30*mu^2)*n^4+82944*(-5+mu)*n^5+110592*n^6)*myP2fast[-2+n,mu]+(-9+mu+8*n)*(-((-11+mu)*(-7+mu)*(-5+mu)*(-3+mu)*(-941137562880+1369543037568*mu-856059425680*mu^2+301467356208*mu^3-65925560840*mu^4+9300152544*mu^5-851420265*mu^6+49707939*mu^7-1788230*mu^8+38538*mu^9-505*mu^10+3*mu^11))+(-9+mu)*(2174231624313600-4271307638939136*mu+3746500640981808*mu^2-1938172937860384*mu^3+658024132807528*mu^4-154336161708664*mu^5+25631896940311*mu^6-3038647883536*mu^7+255911958856*mu^8-15059474264*mu^9+601933862*mu^10-15728672*mu^11+258008*mu^12-2528*mu^13+11*mu^14)*n-4*(-41108205131322624+79558217840920896*mu-69190984849287408*mu^2+35769692404688632*mu^3-12252335726377252*mu^4+2933722316842738*mu^5-504752475079572*mu^6+63145862893203*mu^7-5745369671196*mu^8+376356342416*mu^9-17384266580*mu^10+548066954*mu^11-11274720*mu^12+143142*mu^13-1032*mu^14+3*mu^15)*n^2+4*(-9+mu)*(23845345590072960-40269695568954624*mu+30117142128190992*mu^2-13158415762916400*mu^3+3730778330679232*mu^4-721067843021868*mu^5+97108500711985*mu^6-9153045269192*mu^7+597928404668*mu^8-26432573136*mu^9+759984806*mu^10-13416552*mu^11+134684*mu^12-676*mu^13+mu^14)*n^3+16*(195243602402676096-325694901477820032*mu+242221596032134128*mu^2-106098978486724128*mu^3+30459989915673992*mu^4-6033975669037412*mu^5+845417566861997*mu^6-84452424919988*mu^7+5983741160080*mu^8-295319349276*mu^9+9821158066*mu^10-208695676*mu^11+2610160*mu^12-16816*mu^13+41*mu^14)*n^4+64*(-9+mu)*(14622810947299008-20883005872697088*mu+13042640269010160*mu^2-4687978533249048*mu^3+1073821472622084*mu^4-163965505744412*mu^5+16961587465549*mu^6-1184203363074*mu^7+54575767659*mu^8-1588856808*mu^9+27087171*mu^10-236578*mu^11+785*mu^12)*n^5+256*(68159047060299744-96275531839385520*mu+59912949582646848*mu^2-21651596638546640*mu^3+5041402403618604*mu^4-793018597591700*mu^5+85896040596299*mu^6-6405365947182*mu^7+323083532589*mu^8-10605978520*mu^9+211271829*mu^10-2240614*mu^11+9319*mu^12)*n^6+2048*(-9+mu)*(1541341241341668-1813373921002968*mu+915446118884163*mu^2-259814716685092*mu^3+45629241741242*mu^4-5143009129752*mu^5+373337413062*mu^6-17038328436*mu^7+461072406*mu^8-6556280*mu^9+36355*mu^10)*n^7+8192*(4500207031276008-5239264901634576*mu+2640189965261667*mu^2-755987780804488*mu^3+135697154047598*mu^4-15878627119200*mu^5+1219280284095*mu^6-60190646760*mu^7+1809241320*mu^8-29487896*mu^9+194248*mu^10)*n^8+98304*(-9+mu)*(43421763841182-40431075715248*mu+15676365711905*mu^2-3287037266982*mu^3+404944404503*mu^4-29779385976*mu^5+1265065310*mu^6-28070508*mu^7+245586*mu^8)*n^9+393216*(81960492523446-75530247171240*mu+29304166747543*mu^2-6232480720254*mu^3+791707261321*mu^4-61212289536*mu^5+2795677186*mu^6-68402040*mu^7+679074*mu^8)*n^10+42467328*(-9+mu)*(56814548324-39259013448*mu+10716187369*mu^2-1468655040*mu^3+105783400*mu^4-3770148*mu^5+51547*mu^6)*n^11+56623104*(206001269260-140889461280*mu+38650036817*mu^2-5426533428*mu^3+409650908*mu^4-15687096*mu^5+236549*mu^6)*n^12+42807066624*(-9+mu)*(12790352-5830560*mu+935607*mu^2-61962*mu^3+1425*mu^4)*n^13+24461180928*(65459144-29536992*mu+4812879*mu^2-334554*mu^3+8337*mu^4)*n^14+7044820107264*(-9+mu)*(6091-1374*mu+69*mu^2)*n^15+7044820107264*(9242-2064*mu+111*mu^2)*n^16+760840571584512*(-9+mu)*n^17+338151365148672*n^18)*myP2fast[-1+n,mu])/((-1+n)*(-1+2*n)*(-4+mu+4*n)*(-3+mu+4*n)*(-2+mu+4*n)*(-1+mu+4*n)*(-5+mu+6*n)*(-3+mu+6*n)*(-1+mu+6*n)*(-13+mu+8*n)*(2136180-963208*mu+161921*mu^2-12281*mu^3+391*mu^4-3*mu^5+2*(-13+mu)*(298788-89728*mu+8309*mu^2-230*mu^3+mu^4)*n+8*(1475028-447124*mu+45455*mu^2-1742*mu^3+19*mu^4)*n^2+192*(-13+mu)*(3841-582*mu+15*mu^2)*n^3+768*(5723-876*mu+30*mu^2)*n^4+82944*(-13+mu)*n^5+110592*n^6)),myP2fast[n,mu],OreAlgebra[S[n]]]];

Preparing the stage

(* Alternative closed form for D_{1,1}(n) *)
{FFeAlt,FFoAlt}={FFe,FFo}/.Pochhammer[j,j-1]->Pochhammer[1/2,j-1]*2^(2j-2)/.
Pochhammer[mu+2j+1,j-1]/Pochhammer[mu/2+j+1/2,j-1]->2^(2j-3)*Pochhammer[mu/2+j+1,j-2]/Pochhammer[mu+3j,j-2]/.
Pochhammer[mu+2j,j]^2/Pochhammer[mu/2+j+1/2,j]^2->2^(2j)*Pochhammer[mu/2+j,Floor[(j+1)/2]]^2/Pochhammer[mu/2+Floor[3j/2]+1/2,Floor[(j+1)/2]]^2/.
Pochhammer[mu+2j,j]/Pochhammer[mu/2+j+1/2,j]->2^j*Pochhammer[mu/2+j,Floor[(j+1)/2]]/Pochhammer[mu/2+Floor[3j/2]+1/2,Floor[(j+1)/2]]/.
Pochhammer[mu+2j,j]/Pochhammer[mu/2+j+1/2,j-1]->2^j*Pochhammer[mu/2+j,Floor[(j+1)/2]]/Pochhammer[mu/2+Floor[3j/2]+1/2,Floor[(j-1)/2]]/.
prod[2^a_*b_,c_]:>Product[2^a,c]*prod[b,c]/.
prod[a_/2,{j,1,k-1}]^2*If[k==0,b_,c_]->2^(2-2k)*If[k==0,b/4,c]*prod[a,{j,1,k-1}]^2/.
2^a_:>2^Expand[a];

(* Test *)
Table[Together[(FFeAlt/FFe)/.{sum->Sum,prod->Product}],{n,2,10,2}]
Table[Together[(FFoAlt/FFo)/.{sum->Sum,prod->Product}],{n,1,9,2}]

even n

FFe1=FFeAlt/.prod[a_,{j,k,b_}]:>prod[a,{j,1,b}]/prod[a,{j,1,k-1}]*if[k==0,a/.j->0,1]//.prod[a1_,b_]^c1_.*prod[a2_,b_]^c2_.->prod[a1^c1*a2^c2,b]//.sum[(a_/;FreeQ[a,k])*b_,c_]->a*sum[b,c]//.
prod[a_*b_,{j,1,n/2+c_.}]->prod[a,{j,1,n/2+c}]*prod[b,{j,1,n/2+c}];
Check0=FFe1/(
(* myC[n] *)
((-1)^n+3)/2*prod[Floor[i/2]!/i!,{i,1,n}]*
(* myE[n,mu] *)
Pochhammer[mu+1,n]*prod[(mu+2i+6)^(2*Floor[(i+2)/3]),{i,1,Floor[3/2*Floor[(n-1)/2]-2]}]*
prod[(mu+2i+2*Floor[3/2*Floor[n/2+1]]-1)^(2*Floor[Floor[n/2]/2-(i-1)/3]-1),{i,1,Floor[3/2*Floor[n/2]-2]}]*
(* myF[0,n,mu] *)
prod[(mu+2i+n)^(1-2i),{i,1,Floor[(n-1)/4]}]*prod[(mu-2i+2n+1)^(1-2i),{i,1,Floor[n/4-1]}]);
Check0=Check0/.{(-1)^n->1,Floor[n/2]->n/2,Floor[(n-1)/2]->n/2-1}/.Floor[a_]:>Floor[Together[a]];

n=0 mod 4

Check00=Check0/.{Floor[3/4*(n-2)]->3/4n-2,Floor[n/4]->n/4,Floor[3/4n]->3/4n,Floor[(n-1)/4]->n/4-1,Floor[3/4(n+2)]->3/4n+1,Floor[(3n-4i+4)/12]->n/4+Floor[(1-i)/3]}/.
prod[Pochhammer[j+(mu/2),Floor[((1+j)/2)]],{j,1,(n/2)}]*prod[Pochhammer[j+(mu/2),Floor[((1+j)/2)]],{j,1,-1+(n/2)}]/(prod[(((6+2 i+mu))^(2 Floor[((2+i)/3)])),{i,1,-4+((3 n)/4)}]*prod[(((2 i+mu+n))^(1-2 i)),{i,1,-1+(n/4)}])->Pochhammer[mu/2+1,n/2]*Pochhammer[mu/2+1,n/2-1]/2^(n^2/8-3/4n+1)/.
prod[Pochhammer[(1/2)+2 j+(mu/2),-1+j],{j,1,(n/2)}] prod[Pochhammer[(3/2)+2 j+(mu/2),1+j],{j,1,-1+(n/2)}]prod[(1/(Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[(1/2) (-1+j)]])),{j,1,(n/2)}] prod[(1/(Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[((1+j)/2)]])),{j,1,-1+(n/2)}]/(prod[(((-1+2 i+mu+2 (1+((3 n)/4))))^(-1+2 ((n/4)+Floor[((1-i)/3)]))),{i,1,-2+((3 n)/4)}] prod[(((1-2 i+mu+2 n))^(1-2 i)),{i,1,-1+(n/4)}])->Pochhammer[mu/2+3/4n+1/2,3/4n-1]/(mu+3)/2^(n^2/8-n/2)/.
prod[(1/(Pochhammer[j,j])),{j,1,(n/2)}] prod[(1/(Pochhammer[1+j,1+j])),{j,1,-1+(n/2)}] /prod[((Floor[(i/2)]!)/(i!)),{i,1,n}]->2^(n/2)/.
(2^a_):>2^Expand[a]/.
If->if//.if[a_,b1_,c1_]*if[a_,b2_,c2_]:>if[a,Together[b1*b2],Together[c1*c2]]/.
Pochhammer[1+(mu/2),(n/2)]/Pochhammer[1+mu,n]->1/2^n/Pochhammer[mu/2+1/2,n/2]/.
(* Now the product expression inside the sum *)
(* We first rewrite this Pochhammer to separate even and odd factors *)
Pochhammer[3 j+mu,-2+j]->Pochhammer[mu/2+Floor[3/2j+1/2],Floor[(j-2)/2]]*Pochhammer[mu/2+Floor[3/2j]+1/2,Floor[(j-1)/2]]*2^(j-2)/.
prod[a_Times,b_]:>(prod[#,b]&/@a)/.
prod[(1/((Pochhammer[j+(mu/2),Floor[((1+j)/2)]])^2)),{j,1,-1+k}] prod[((Pochhammer[1+j+(mu/2),-2+j])^2),{j,1,-1+k}] prod[(1/((Pochhammer[(mu/2)+Floor[(1/2)+((3 j)/2)],Floor[(1/2) (-2+j)]])^2)),{j,1,-1+k}]->if[k==0,4/mu^2,1]/Pochhammer[mu/2+1,k-1]^2/.
prod[Pochhammer[(1/2)+2 j+(mu/2),-1+j],{j,1,-1+k}] prod[(1/(Pochhammer[(3/2)+2 j+(mu/2),1+j])),{j,1,-1+k}] prod[(1/(Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[(1/2) (-1+j)]])),{j,1,-1+k}] prod[Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[((1+j)/2)]],{j,1,-1+k}]->if[k==0,1,(mu+3)/2]/Pochhammer[mu/2+2k-1/2,k]/.
prod[(1/((Pochhammer[(1/2),-1+j])^2)),{j,1,-1+k}] prod[Pochhammer[j,j],{j,1,-1+k}] prod[Pochhammer[1+j,1+j],{j,1,-1+k}]->if[k==0,1/8,1]*2^(2k(k-1))*Pochhammer[3/2,k-1]*Pochhammer[1/2,k-1]^2/.
prod[2^a_,{j,1,k-1}]:>With[{cf=Product[2^a,{j,1,k-1}]},cf*if[k==0,1/(cf/.k->0),1]]//.
if[k==0,a1_,b1_]*if[k==0,a2_,b2_]->if[k==0,a1*a2,b1*b2]/.
if[k==0,a_*b_,a_*c_]->a*if[k==0,b,c]/.
a_*sum[b_,c_]->sum[a*b,c]/.
Pochhammer[mu/2+1,n/2-1]->Pochhammer[mu/2+1,k-1]*Pochhammer[mu/2+k,n/2-k]/.
Pochhammer[3/2,k-1]*if[k==0,a_,b_]->if[k==0,a,b*(2k-1)!/(2k-2)!!]/2^(k-1)/.(2^a_):>2^FullSimplify[a]/.if[k==0,a_,b_]->4*if[k==0,a/4,b/4]

Table[Together[(Check00/.{sum->Sum,prod->Product,if->If})/myP2[n/4,mu]],{n,4,20,4}]

TraditionalForm[HoldForm@@{Check00}/.{sum->Sum,prod->Product,if->If}]

(* Our expression fits the recurrence (this is the initial value check). *)
test=ApplyOreOperator[rec2,f[n]];
Together[Table[test,{n,0,4}]/.f[nn_]:>(Check00/.n->4nn/.sum->Sum/.if->If)]

(* the smnd for k=1 to n/2 *)
smnd=ExpandAll[Check00[[1]]/.n->4n/.if[k==0,_,a_]->a]

Factor[{{op},{cert}}=CreativeTelescoping[smnd,S[k]-1,S[n]]]

(* We use the following identity in order to construct an annihilating operator for the sum *)
TraditionalForm[HoldForm[
0==Sum[(p1[n]S[n]+p0[n])*f[n,k],{k,1,2n}]==
(p1[n]S[n]+p0[n])*Sum[f[n,k],{k,1,2n}]
-p1[n]*(f[n+1,2n+1]+f[n+1,2n+2])
]]
test=ApplyOreOperator[op,Check00/.n->4n/.sum[a_,{k,0,b_}]->sum[a,{k,1,b}]]-LeadingCoefficient[op]*((smnd/.n->n+1/.k->2n+1)+(smnd/.n->n+1/.k->2n+2));
Together[Table[test/.sum->Sum/.if->If,{n,0,4}]]

inh=Factor[LeadingCoefficient[op]*(1+FullSimplify[(smnd/.n->n+1/.k->2n+2)/(smnd/.n->n+1/.k->2n+1)])]*(smnd/.n->n+1/.k->2n+1);
rec00=Annihilator[inh,S[n]][[1]]**op;

smnd0=ExpandAll[Together[Check00[[1]]/.n->4n/.if[k==0,a_,_]->a/.k->0]]

(* smnd0 is also in the inhomogeneous part, but is already annihilated by rec2a. *)
OreReduce[rec00,Annihilator[smnd0,S[n]]]

GBEqual[{rec2},{rec00}]

n=2 mod 4

Check01=Check0/.Floor[a_]:>(FullSimplify[Floor[a/.n->4l-2],Element[l,Integers]]/.Floor[l+b_]:>Floor[Together[b]]+l/.l->(n+2)/4)/.
prod[a_^b_,{i,1,c_}]:>prod[Expand[a]^Simplify[b],{i,1,Expand[c]}]/.
prod[Pochhammer[j+(mu/2),Floor[((1+j)/2)]],{j,1,-1+(n/2)}]*prod[Pochhammer[j+(mu/2),Floor[((1+j)/2)]],{j,1,(n/2)}]/(prod[(((6+2 i+mu))^(2 Floor[((2+i)/3)])),{i,1,-(7/2)+((3 n)/4)}]prod[(((2 i+mu+n))^(1-2 i)),{i,1,-(1/2)+(n/4)}] )->
Pochhammer[mu/2+1,n/2]*Pochhammer[mu/2+1,n/2-1]/2^(n^2/8-3/4n+1)/.
prod[Pochhammer[(1/2)+2 j+(mu/2),-1+j],{j,1,(n/2)}] prod[Pochhammer[(3/2)+2 j+(mu/2),1+j],{j,1,-1+(n/2)}] prod[(1/(Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[(1/2) (-1+j)]])),{j,1,(n/2)}] prod[(1/(Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[((1+j)/2)]])),{j,1,-1+(n/2)}]/(prod[(((2+2 i+mu+((3 n)/2)))^((1/2) (n-4 Ceiling[(1/6)+(i/3)]))),{i,1,-(5/2)+((3 n)/4)}] prod[(((1-2 i+mu+2 n))^(1-2 i)),{i,1,-(3/2)+(n/4)}])->
Pochhammer[mu/2+3/4n,3/4n-1/2]/(mu+3)/2^(n^2/8-n/2-1/2)/.
prod[(1/(Pochhammer[j,j])),{j,1,(n/2)}] prod[(1/(Pochhammer[1+j,1+j])),{j,1,-1+(n/2)}] /prod[((Floor[(i/2)]!)/(i!)),{i,1,n}]->2^(n/2)/.
(2^a_):>2^Expand[a]/.
If->if//.if[a_,b1_,c1_]*if[a_,b2_,c2_]:>if[a,Together[b1*b2],Together[c1*c2]]/.
Pochhammer[1+(mu/2),(n/2)]/Pochhammer[1+mu,n]->1/2^n/Pochhammer[mu/2+1/2,n/2]/.
(* Now the product expression inside the sum *)
(* We first rewrite this Pochhammer to separate even and odd factors *)
Pochhammer[3 j+mu,-2+j]->Pochhammer[mu/2+Floor[3/2j+1/2],Floor[(j-2)/2]]*Pochhammer[mu/2+Floor[3/2j]+1/2,Floor[(j-1)/2]]*2^(j-2)/.
prod[a_Times,b_]:>(prod[#,b]&/@a)/.
prod[(1/((Pochhammer[j+(mu/2),Floor[((1+j)/2)]])^2)),{j,1,-1+k}] prod[((Pochhammer[1+j+(mu/2),-2+j])^2),{j,1,-1+k}] prod[(1/((Pochhammer[(mu/2)+Floor[(1/2)+((3 j)/2)],Floor[(1/2) (-2+j)]])^2)),{j,1,-1+k}]->if[k==0,4/mu^2,1]/Pochhammer[mu/2+1,k-1]^2/.
prod[Pochhammer[(1/2)+2 j+(mu/2),-1+j],{j,1,-1+k}] prod[(1/(Pochhammer[(3/2)+2 j+(mu/2),1+j])),{j,1,-1+k}] prod[(1/(Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[(1/2) (-1+j)]])),{j,1,-1+k}] prod[Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[((1+j)/2)]],{j,1,-1+k}]->if[k==0,1,(mu+3)/2]/Pochhammer[mu/2+2k-1/2,k]/.
prod[(1/((Pochhammer[(1/2),-1+j])^2)),{j,1,-1+k}] prod[Pochhammer[j,j],{j,1,-1+k}] prod[Pochhammer[1+j,1+j],{j,1,-1+k}]->if[k==0,1/8,1]*2^(2k(k-1))*Pochhammer[3/2,k-1]*Pochhammer[1/2,k-1]^2/.
prod[2^a_,{j,1,k-1}]:>With[{cf=Product[2^a,{j,1,k-1}]},cf*if[k==0,1/(cf/.k->0),1]]//.
if[k==0,a1_,b1_]*if[k==0,a2_,b2_]->if[k==0,a1*a2,b1*b2]/.
if[k==0,a_*b_,a_*c_]->a*if[k==0,b,c]/.
a_*sum[b_,c_]->sum[a*b,c]/.
Pochhammer[mu/2+1,n/2-1]->Pochhammer[mu/2+1,k-1]*Pochhammer[mu/2+k,n/2-k]/.
Pochhammer[3/2,k-1]*if[k==0,a_,b_]->if[k==0,a,b*(2k-1)!/(2k-2)!!]/2^(k-1)/.(2^a_):>2^FullSimplify[a]/.if[k==0,a_,b_]->4*if[k==0,a/4,b/4]

Table[Together[(Check01/.{sum->Sum,prod->Product,if->If})/myP1[n/4+1/2,mu]],{n,2,18,4}]

TraditionalForm[HoldForm@@{Check01}/.{sum->Sum,prod->Product,if->If}]

(* Our expression fits the recurrence (this is the initial value check). *)
test=ApplyOreOperator[rec1,f[n]];
Together[Table[test,{n,5}]/.f[nn_]:>(Check01/.n->4nn-2/.sum->Sum/.if->If)]

(* the smnd for k>=1 *)
smnd=Check01[[1]]/.n->4n-2/.if[k==0,_,a_]->a

Factor[{{op},{cert}}=CreativeTelescoping[smnd,S[k]-1,S[n]]]

(* We use the following identity in order to construct an annihilating operator for the sum *)
TraditionalForm[HoldForm[
0==Sum[(p1[n]S[n]+p0[n])*f[n,k],{k,1,2n-1}]==
(p1[n]S[n]+p0[n])*Sum[f[n,k],{k,1,2n-1}]
-p1[n]*(f[n+1,2n]+f[n+1,2n+1])
]]
test=ApplyOreOperator[op,Check01/.n->4n-2/.sum[a_,{k,0,b_}]->sum[a,{k,1,b}]]-LeadingCoefficient[op]*((smnd/.n->n+1/.k->2n)+(smnd/.n->n+1/.k->2n+1));
Together[Table[test/.sum->Sum/.if->If,{n,5}]]

inh=Factor[LeadingCoefficient[op]*(1+FullSimplify[(smnd/.n->n+1/.k->2n+1)/(smnd/.n->n+1/.k->2n)])]*(smnd/.n->n+1/.k->2n);
rec01=Annihilator[inh,S[n]][[1]]**op;

smnd0=Together[ExpandAll[Check01[[1]]/.n->4n-2/.if[k==0,a_,_]->a/.k->0]]

(* smnd0 is also in the inhomogeneous part, but is already annihilated by rec01. *)
OreReduce[rec01,Annihilator[smnd0,S[n]]]

GBEqual[{rec1},rec01]

odd n

FFo1=FFoAlt/.prod[a_,{j,k,b_}]:>prod[a,{j,1,b}]/prod[a,{j,1,k-1}]*if[k==0,a/.j->0,1]//.prod[a1_,b_]^c1_.*prod[a2_,b_]^c2_.->prod[a1^c1*a2^c2,b]//.sum[(a_/;FreeQ[a,k])*b_,c_]->a*sum[b,c]//.
prod[a_*b_,{j,1,n/2+c_.}]->prod[a,{j,1,n/2+c}]*prod[b,{j,1,n/2+c}];
Check1=FFo1/(
(* myC[n] *)
((-1)^n+3)/2*prod[Floor[i/2]!/i!,{i,1,n}]*
(* myE[n,mu] *)
Pochhammer[mu+1,n]*prod[(mu+2i+6)^(2*Floor[(i+2)/3]),{i,1,Floor[3/2*Floor[(n-1)/2]-2]}]*
prod[(mu+2i+2*Floor[3/2*Floor[n/2+1]]-1)^(2*Floor[Floor[n/2]/2-(i-1)/3]-1),{i,1,Floor[3/2*Floor[n/2]-2]}]*
(* myF[1,n,mu] *)
prod[(mu+2i+n+1)^(-2i),{i,1,Floor[(n-1)/4]}]*prod[(mu-2i+2n-1)^(-2i),{i,1,Floor[n/4-1]}]*
(* remaining factor in myF[n,mu] *)
Product[mu+2n+2i-1,{i,1,(n-5)/2}]);
Check1=Check1/.{(-1)^n->-1,Floor[n/2]->(n-1)/2,Floor[(n-1)/2]->(n-1)/2}/.Floor[a_]:>Floor[Together[a]];

(* n=1 and n=3 don't work because of the "remaining factor in myF[n,mu]". *)
Table[Together[(Check1/.{sum->Sum,prod->Product,if->If})/myP1[(n+3)/4,mu]],{n,1,19,4}]
Table[Together[(Check1/.{sum->Sum,prod->Product,if->If})/myP2[(n+1)/4,mu]],{n,3,21,4}]

n=1 mod 4

Check10=Check1/.Floor[a_]:>(FullSimplify[Floor[a/.n->4l-3],Element[l,Integers]]/.Floor[l+b_]:>Floor[Together[b]]+l/.l->(n+3)/4)/.
prod[a_^b_,{i,1,c_}]:>prod[Expand[a]^Simplify[b],{i,1,Expand[c]}]/.
prod[((Pochhammer[j+(mu/2),Floor[((1+j)/2)]])^2),{j,1,-(1/2)+(n/2)}]/prod[(((6+2 i+mu))^(2 Floor[((2+i)/3)])),{i,1,-(11/4)+((3 n)/4)}]/prod[(((1+2 i+mu+n))^(-2 i)),{i,1,-(1/4)+(n/4)}]->Pochhammer[mu/2+1,(n-1)/2]^2/2^((n^2-6n+5)/8)/.
prod[Pochhammer[-(1/2)+2 j+(mu/2),j],{j,1,-(1/2)+(n/2)}]prod[Pochhammer[(3/2)+2 j+(mu/2),1+j],{j,1,-(1/2)+(n/2)}]prod[(1/((Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[((1+j)/2)]])^2)),{j,1,-(1/2)+(n/2)}]/prod[(((-(1/2)+2 i+mu+((3 n)/2)))^((1/2) (1+n-4 Ceiling[((2+i)/3)]))),{i,1,-(11/4)+((3 n)/4)}]/prod[(((-1-2 i+mu+2 n))^(-2 i)),{i,1,-(5/4)+(n/4)}]->Pochhammer[mu/2+3/4(n-1)+3/2,3/4(n-1)-2]*Pochhammer[mu/2+(n-1)+3/2,(n-1)/2+1]/(mu+3)/2^((n^2-8n+15)/8)/.
prod[(1/(Pochhammer[j,j])),{j,1,-(1/2)+(n/2)}]prod[(1/(Pochhammer[1+j,1+j])),{j,1,-(1/2)+(n/2)}]/prod[((Floor[(i/2)]!)/(i!)),{i,1,n}]->2^((n-1)/2)/.
(2^a_):>2^Expand[a]/.
If->if//.if[a_,b1_,c1_]*if[a_,b2_,c2_]:>if[a,Together[b1*b2],Together[c1*c2]]/.
Pochhammer[(1/2)+(mu/2)+n,1+(1/2) (-1+n)]Pochhammer[(3/2)+(mu/2)+(3/4) (-1+n),-2+(3/4) (-1+n)]/Pochhammer[(1/2) (1+mu+2 n),(1/2) (-5+n)]->Pochhammer[(3/2)+(mu/2)+(3/4) (-1+n),(3n+1)/4]/.
Pochhammer[1+(mu/2),((n-1)/2)]^2/Pochhammer[1+mu,n]->Pochhammer[mu/2+1,(n-1)/2]/2^n/Pochhammer[mu/2+1/2,(n+1)/2]/.
(* Now the product expression inside the sum *)
(* We first rewrite this Pochhammer to separate even and odd factors *)
Pochhammer[3 j+mu,-2+j]->Pochhammer[mu/2+Floor[3/2j+1/2],Floor[(j-2)/2]]*Pochhammer[mu/2+Floor[3/2j]+1/2,Floor[(j-1)/2]]*2^(j-2)/.
prod[a_Times,b_]:>(prod[#,b]&/@a)/.
prod[(1/((Pochhammer[j+(mu/2),Floor[((1+j)/2)]])^2)),{j,1,-1+k}] prod[((Pochhammer[1+j+(mu/2),-2+j])^2),{j,1,-1+k}] prod[(1/((Pochhammer[(mu/2)+Floor[(1/2)+((3 j)/2)],Floor[(1/2) (-2+j)]])^2)),{j,1,-1+k}]->if[k==0,4/mu^2,1]/Pochhammer[mu/2+1,k-1]^2/.
prod[(1/(Pochhammer[-(1/2)+2 j+(mu/2),j])),{j,1,-1+k}] prod[((Pochhammer[(1/2)+2 j+(mu/2),-1+j])^2),{j,1,-1+k}] prod[(1/(Pochhammer[(3/2)+2 j+(mu/2),1+j])),{j,1,-1+k}] prod[(1/((Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[(1/2) (-1+j)]])^2)),{j,1,-1+k}] prod[((Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[((1+j)/2)]])^2),{j,1,-1+k}]->if[k==0,1,(mu+3)/2]/Pochhammer[mu/2+2k-1/2,k]/.
prod[(1/((Pochhammer[(1/2),-1+j])^2)),{j,1,-1+k}] prod[Pochhammer[j,j],{j,1,-1+k}] prod[Pochhammer[1+j,1+j],{j,1,-1+k}]->if[k==0,1/8,1]*2^(2k(k-1))*Pochhammer[3/2,k-1]*Pochhammer[1/2,k-1]^2/.
prod[2^a_,{j,1,k-1}]:>With[{cf=Product[2^a,{j,1,k-1}]},cf*if[k==0,1/(cf/.k->0),1]]//.
if[k==0,a1_,b1_]*if[k==0,a2_,b2_]->if[k==0,a1*a2,b1*b2]/.
if[k==0,a_*b_,a_*c_]->a*if[k==0,b,c]/.
a_*sum[b_,c_]->sum[a*b,c]/.
Pochhammer[mu/2+1,(n-1)/2]->Pochhammer[mu/2+1,k-1]*Pochhammer[mu/2+k,(n+1)/2-k]/.
Pochhammer[3/2,k-1]*if[k==0,a_,b_]->if[k==0,a,b*(2k-1)!/(2k-2)!!]/2^(k-1)/.(2^a_):>2^FullSimplify[a]/.Pochhammer[(3/2)+(mu/2)+(3/4) (-1+n),(1/4) (1+3 n)]->Pochhammer[(3/4)+(mu/2)+((3 n)/4),1/4+3 n/4]/.if[k==0,a_,b_]->4*if[k==0,a/4,b/4]

Table[Factor[(Check10/.{sum->Sum,prod->Product,if->If})/myP1[(n+3)/4,mu]],{n,1,19,4}]

TraditionalForm[HoldForm@@{Check10}/.{sum->Sum,prod->Product,if->If}]

(* Our expression fits the recurrence (this is the initial value check). *)
test=ApplyOreOperator[rec1,f[n]];
Together[Table[test,{n,5}]/.f[nn_]:>(Check10/.n->4nn-3/.sum->Sum/.if->If)]

(* the smnd for k>=1 *)
smnd=Check10[[1]]/.n->4n-3/.if[k==0,_,a_]->a/.Pochhammer[a_,b_]:>Pochhammer@@Expand[{a,b}]

Factor[{{op},{cert}}=CreativeTelescoping[smnd,S[k]-1,S[n]]]

(* We use the following identity in order to construct an annihilating operator for the sum *)
TraditionalForm[HoldForm[
0==Sum[(p1[n]S[n]+p0[n])*f[n,k],{k,1,2n-1}]==
(p1[n]S[n]+p0[n])*Sum[f[n,k],{k,1,2n-1}]
-p1[n]*(f[n+1,2n]+f[n+1,2n+1])
]]
test=ApplyOreOperator[op,Check10/.n->4n-3/.sum[a_,{k,0,b_}]->sum[a,{k,1,b}]]-LeadingCoefficient[op]*((smnd/.n->n+1/.k->2n)+(smnd/.n->n+1/.k->2n+1));
Together[Table[test/.sum->Sum/.if->If,{n,5}]]

inh=Factor[LeadingCoefficient[op]*(1+FullSimplify[(smnd/.n->n+1/.k->2n+1)/(smnd/.n->n+1/.k->2n)])]*(smnd/.n->n+1/.k->2n);
rec10=Annihilator[inh,S[n]][[1]]**op;

smnd0=Together[ExpandAll[Check10[[1]]/.n->4n-3/.if[k==0,a_,_]->a/.k->0]]

(* smnd0 is also in the inhomogeneous part, but is already annihilated by rec01. *)
OreReduce[rec10,Annihilator[smnd0,S[n]]]

GBEqual[{rec1},rec10]

n=3 mod 4

Check11=Check1/.Floor[a_]:>(FullSimplify[Floor[a/.n->4l-1],Element[l,Integers]]/.Floor[l+b_]:>Floor[Together[b]]+l/.l->(n+1)/4)/.
prod[a_^b_,{i,1,c_}]:>prod[Expand[a]^Simplify[b],{i,1,Expand[c]}]/.
prod[((Pochhammer[j+(mu/2),Floor[((1+j)/2)]])^2),{j,1,-(1/2)+(n/2)}]/prod[(((6+2 i+mu))^(2 Floor[((2+i)/3)])),{i,1,-(13/4)+((3 n)/4)}] /prod[(((1+2 i+mu+n))^(-2 i)),{i,1,-(3/4)+(n/4)}]->Pochhammer[mu/2+1,(n-1)/2]^2/2^((n^2-6n+9)/8)/.
prod[Pochhammer[-(1/2)+2 j+(mu/2),j],{j,1,-(1/2)+(n/2)}] prod[Pochhammer[(3/2)+2 j+(mu/2),1+j],{j,1,-(1/2)+(n/2)}] prod[(1/((Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[((1+j)/2)]])^2)),{j,1,-(1/2)+(n/2)}]/prod[((((1/2)+2 i+mu+((3 n)/2)))^((1/2) (-1+n-4 Ceiling[(1/6)+(i/3)]))),{i,1,-(13/4)+((3 n)/4)}] /prod[(((-1-2 i+mu+2 n))^(-2 i)),{i,1,-(7/4)+(n/4)}]->Pochhammer[mu/2+3/4n+5/4,(3n-13)/4]Pochhammer[mu/2+n+1/2,(n+1)/2]/(mu+3)/2^((n^2-8n+15)/8)/.
prod[(1/(Pochhammer[j,j])),{j,1,-(1/2)+(n/2)}]prod[(1/(Pochhammer[1+j,1+j])),{j,1,-(1/2)+(n/2)}]/prod[((Floor[(i/2)]!)/(i!)),{i,1,n}]->2^((n-1)/2)/.
(2^a_):>2^Expand[a]/.
If->if//.if[a_,b1_,c1_]*if[a_,b2_,c2_]:>if[a,Together[b1*b2],Together[c1*c2]]/.
 Pochhammer[(1/2)+(mu/2)+n,((1+n)/2)]Pochhammer[(5/4)+(mu/2)+((3 n)/4),(1/4) (-13+3 n)]/Pochhammer[(1/2) (1+mu+2 n),(1/2) (-5+n)]->Pochhammer[(5/4)+(mu/2)+((3 n)/4),3n/4-1/4]/.
Pochhammer[1+(mu/2),((n-1)/2)]^2/Pochhammer[1+mu,n]->Pochhammer[mu/2+1,(n-1)/2]/2^n/Pochhammer[mu/2+1/2,(n+1)/2]/.
(* Now the product expression inside the sum *)
(* We first rewrite this Pochhammer to separate even and odd factors *)
Pochhammer[3 j+mu,-2+j]->Pochhammer[mu/2+Floor[3/2j+1/2],Floor[(j-2)/2]]*Pochhammer[mu/2+Floor[3/2j]+1/2,Floor[(j-1)/2]]*2^(j-2)/.
prod[a_Times,b_]:>(prod[#,b]&/@a)/.
prod[(1/((Pochhammer[j+(mu/2),Floor[((1+j)/2)]])^2)),{j,1,-1+k}] prod[((Pochhammer[1+j+(mu/2),-2+j])^2),{j,1,-1+k}] prod[(1/((Pochhammer[(mu/2)+Floor[(1/2)+((3 j)/2)],Floor[(1/2) (-2+j)]])^2)),{j,1,-1+k}]->if[k==0,4/mu^2,1]/Pochhammer[mu/2+1,k-1]^2/.
prod[(1/(Pochhammer[-(1/2)+2 j+(mu/2),j])),{j,1,-1+k}] prod[((Pochhammer[(1/2)+2 j+(mu/2),-1+j])^2),{j,1,-1+k}] prod[(1/(Pochhammer[(3/2)+2 j+(mu/2),1+j])),{j,1,-1+k}] prod[(1/((Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[(1/2) (-1+j)]])^2)),{j,1,-1+k}] prod[((Pochhammer[(1/2)+(mu/2)+Floor[((3 j)/2)],Floor[((1+j)/2)]])^2),{j,1,-1+k}]->if[k==0,1,(mu+3)/2]/Pochhammer[mu/2+2k-1/2,k]/.
prod[(1/((Pochhammer[(1/2),-1+j])^2)),{j,1,-1+k}] prod[Pochhammer[j,j],{j,1,-1+k}] prod[Pochhammer[1+j,1+j],{j,1,-1+k}]->if[k==0,1/8,1]*2^(2k(k-1))*Pochhammer[3/2,k-1]*Pochhammer[1/2,k-1]^2/.
prod[2^a_,{j,1,k-1}]:>With[{cf=Product[2^a,{j,1,k-1}]},cf*if[k==0,1/(cf/.k->0),1]]//.
if[k==0,a1_,b1_]*if[k==0,a2_,b2_]->if[k==0,a1*a2,b1*b2]/.
if[k==0,a_*b_,a_*c_]->a*if[k==0,b,c]/.
a_*sum[b_,c_]->sum[a*b,c]/.
Pochhammer[mu/2+1,(n-1)/2]->Pochhammer[mu/2+1,k-1]*Pochhammer[mu/2+k,(n+1)/2-k]/.
Pochhammer[3/2,k-1]*if[k==0,a_,b_]->if[k==0,a,b*(2k-1)!/(2k-2)!!]/2^(k-1)/.(2^a_):>2^FullSimplify[a]/.if[k==0,a_,b_]->4*if[k==0,a/4,b/4]

Table[Factor[(Check11/.{sum->Sum,prod->Product,if->If})/myP2[(n+1)/4,mu]],{n,3,23,4}]

TraditionalForm[HoldForm@@{Check11}/.{sum->Sum,prod->Product,if->If}]

(* Our expression fits the recurrence (this is the initial value check). *)
test=ApplyOreOperator[rec2,f[n]];
Together[Table[test,{n,5}]/.f[nn_]:>(Check11/.n->4nn-1/.sum->Sum/.if->If)]

(* the smnd for k>=1 *)
smnd=Check11[[1]]/.n->4n-1/.if[k==0,_,a_]->a/.Pochhammer[a_,b_]:>Pochhammer@@Expand[{a,b}]

Factor[{{op},{cert}}=CreativeTelescoping[smnd,S[k]-1,S[n]]]

(* We use the following identity in order to construct an annihilating operator for the sum *)
TraditionalForm[HoldForm[
0==Sum[(p1[n]S[n]+p0[n])*f[n,k],{k,1,2n}]==
(p1[n]S[n]+p0[n])*Sum[f[n,k],{k,1,2n}]
-p1[n]*(f[n+1,2n+1]+f[n+1,2n+2])
]]
test=ApplyOreOperator[op,Check11/.n->4n-1/.sum[a_,{k,0,b_}]->sum[a,{k,1,b}]]-LeadingCoefficient[op]*((smnd/.n->n+1/.k->2n+1)+(smnd/.n->n+1/.k->2n+2));
Together[Table[test/.sum->Sum/.if->If,{n,5}]]

inh=Factor[LeadingCoefficient[op]*(1+FullSimplify[(smnd/.n->n+1/.k->2n+2)/(smnd/.n->n+1/.k->2n+1)])]*(smnd/.n->n+1/.k->2n+1);
rec11=Annihilator[inh,S[n]][[1]]**op;

smnd0=Together[ExpandAll[Check11[[1]]/.n->4n-1/.if[k==0,a_,_]->a/.k->0]]

(* smnd0 is also in the inhomogeneous part, but is already annihilated by rec01. *)
OreReduce[rec11,Annihilator[smnd0,S[n]]]

GBEqual[{rec2},rec11]

Summarize

{Check00,Check01,Check10,Check11}/.(if[k==0,(mu-2)/mu^2,1/8/(2k-2)!!]*Pochhammer[1/2,k-1]^2*Pochhammer[mu-1,3k-2]/Pochhammer[mu/2+1,k-1]/Pochhammer[mu/2+k-1/2,k-1]/Pochhammer[mu/2+2k-1/2,k])->(h[n,k]/2^k)

Section 6: The General Determinant

Overview

Do[Fam[i]={},{i,0,6}];
NiceQ[x_]:=Max[First/@TimeConstrained[FactorInteger[x],0.1,{{10^6,0}}]]<100;
Do[
tt=Table[Det[DstMat[s,t,n,37]],{n,17,19}];
Which[
MatchQ[tt,{(0)..}],AppendTo[Fam[0],{s,t}],
MatchQ[tt,{_?NiceQ,0,_?NiceQ}],AppendTo[Fam[2],{s,t}],
MatchQ[tt,{(_?NiceQ)..}],AppendTo[Fam[1],{s,t}],
NiceQ[tt[[2]]]&&s<0,AppendTo[Fam[4],{s,t}],
NiceQ[tt[[2]]],AppendTo[Fam[3],{s,t}],
tt=Rest[tt]/Most[tt];
NiceQ[tt[[1]]],AppendTo[Fam[5],{s,t}],
NiceQ[tt[[2]]],AppendTo[Fam[6],{s,t}]
];
,{s,-3,6},{t,-3,6}];
ListPlot[Table[Append[Fam[i],{i,-6}],{i,0,6}],
PlotStyle->{Black,Red,Blue,Green,Purple},PlotMarkers->{Automatic,Large},ImageSize->600]

Corollary 15:  (Subscript[D, -r,-r])(n)

TableForm[DstMat[-2,-2,10,37]]

TableForm[DstMat[-3,-3,10,37]]

Union[Flatten[Table[Factor[Dst[r,r,n+1]/Dst[r+1,r+1,n]],{r,-5,-1},{n,10}]]]

Table[Together[Dst[-r,-r,n]/If[r<n,Dst[0,0,n-r],1]],{r,0,5},{n,10}]

Proposition 16:  (Subscript[D, s,t])(n)=0  for  t≤-1  and  s≥t+1

TableForm[DstMat[0,-1,8,31]]

Union[Flatten[Table[Dst[s,t,n],{t,-1,-5,-1},{s,t+1,6},{n,10}]]]

Theorem 17: switching s and t

Table[Together[Dst[s,t,n]-Product[Pochhammer[mu+i+s-1,n]/Pochhammer[i+s+1,n],{i,0,t-s-1}]*Dst[t,s,n]],{s,0,5},{t,s,5},{n,8}]

(* Base case n=1 (we may assume t>s, hence Kronecker delta is not present at all). *)
(Dst[s,t,1]/Dst[t,s,1]/.KroneckerDelta[__]->0)/Product[Pochhammer[mu+i+s-1,1]/Pochhammer[i+s+1,1],{i,0,t-s-1}]

FullSimplify[%]

(* Base case n=2 (we may assume t>s). We start by assuming t=s+1. *)
Simplify[Dst[s,t,2]/Dst[t,s,2]/.t->s+1]/Product[Pochhammer[mu+i+s-1,2]/Pochhammer[i+s+1,2],{i,0,s+1-s-1}]

FullSimplify[%]

(* Base case n=2; now assume t≥s+2, hence Kronecker delta is not present at all. *)
Simplify[Dst[s,t,2]/Dst[t,s,2]/.KroneckerDelta[__]->0]/Product[Pochhammer[mu+i+s-1,2]/Pochhammer[i+s+1,2],{i,0,t-s-1}]

FullSimplify[%]

Theorem 18  (Family A):  (Subscript[D, 2r,0])(n)  and  (Subscript[D, 0,2r])(n)

quoAe=Pochhammer[mu+2n+4r,n-r]*Pochhammer[mu/2+2n+r+1/2,n-r-1]/Pochhammer[n-r,n-r]/Pochhammer[mu/2+n+2r+1/2,n-r-1];
quoAo=Pochhammer[mu+2n+4r-2,n-r-1]*Pochhammer[mu/2+2n+r-1/2,n-r]/Pochhammer[n-r,n-r]/Pochhammer[mu/2+n+2r-1/2,n-r-1];
Table[Together[{Dst[2r,0,2n+1]/Dst[2r,0,2n]/quoAe,Dst[2r,0,2n]/Dst[2r,0,2n-1]/quoAo}],{n,5},{r,0,n-1}]

RA[n_]:=If[EvenQ[n],
Pochhammer[mu+n+4r,n/2-r]*Pochhammer[mu/2+n+r+1/2,n/2-r-1]/Pochhammer[n/2-r,n/2-r]/Pochhammer[mu/2+n/2+2r+1/2,n/2-r-1],
Pochhammer[mu+n+4r-1,(n+1)/2-r-1]*Pochhammer[mu/2+n+r+1/2,(n+1)/2-r]/Pochhammer[(n+1)/2-r,(n+1)/2-r]/Pochhammer[mu/2+(n+1)/2+2r-1/2,(n+1)/2-r-1]];
Table[Together[Dst[2r,0,n+1]/Dst[2r,0,n]/RA[n]],{n,10},{r,0,(n-1)/2}]

{quoAei,quoAoi}={quoAe/.n->i/2,quoAo/.n->(i+1)/2};
Table[Together[Dst[2r,0,n]/If[n≤2r,1,2*Product[If[Mod[i,2]==0,quoAei,quoAoi],{i,2r+1,n-1}]]],{n,8},{r,8}]

Table[Dst[2r,0,n],{r,0,4},{n,2r}]

Union[Flatten[Table[Dst[2r,0,n,mu]-Dst[0,0,n-2r,mu+6r],{r,4},{n,2r+1,2r+5},{mu,-2,6}]]]

Table[Together[Dst[0,2r,n]/(Product[Pochhammer[mu+i-1,2r]/Pochhammer[i+1,2r],{i,0,n-1}]*Dst[2r,0,n])],{n,8},{r,8}]

(* Special case: Proposition 8 *)
(quoAe/.r->0)/R00e[n]

(* Special case: Proposition 8 *)
(quoAo/.r->0)/R00o[n]

(* Special case: Lemma 6 *)
(quoAo/.r->1)/R20[n]

FullSimplify[%]

(* Special case: Lemma 7 *)
((quoAo/.r->1)*(mu+2n-2)*(mu+2n-1)/(2n)/(2n+1))/R02[n]

FullSimplify[%]

Theorem 19  (Family B):  (Subscript[D, 2r-1,0])(n)  and  (Subscript[D, 0,2r-1])(n)

quoB=-Pochhammer[mu+2n+4r-4,n-r+1]*Pochhammer[mu+2n+4r-3,n-r]*Pochhammer[mu/2+r+2n-1/2,n-r]^2/Pochhammer[n-r+1,n-r]/Pochhammer[n-r+1,n-r+1]/Pochhammer[mu/2+n+2r-3/2,n-r]^2;
Table[Factor[(Dst[2r-1,0,2n+1]/Dst[2r-1,0,2n-1])/quoB],{n,5},{r,n}]//TableForm

Table[Together[Dst[2r-1,0,n]/Product[quoB,{n,r,(n-1)/2}]],{n,1,9,2},{r,5}]

TableForm[Table[Dst[2r-1,0,n],{n,2,10,2},{r,5}]]

Table[Together[Dst[2r-1,0,n,mu]-Dst[1,0,n-2r+2,mu+6r-6]],{r,5},{n,2r,10}]

Table[Dst[0,2r-1,n,37]/
(Product[Pochhammer[mu+i-1,2r-1]/Pochhammer[i+1,2r-1],{i,0,n-1}]*Product[quoB,{n,r,(n-1)/2}]/.mu->37),{n,1,9,2},{r,1,10}]

(* Special case: Proposition 9 *)
(quoB/.r->1)/R10[n]

(* Special case: Proposition 10 *)
((quoB/.r->1)*(mu+2n-1)*(mu+2n-2)/(2n+1)/(2n))/R01[n]

FullSimplify[%]

Conjecture 20  (Family C):  (Subscript[D, 2r,1])(n) and (Subscript[D, 1,2r])(n)

quoC=-(2n+2r)*(mu+2n+2r-1)*(mu+2n+2r)*Pochhammer[mu+2n+4r,n-r]^2*Pochhammer[mu/2+2n+r+3/2,n-r+1]^2/Pochhammer[n-r+1,n-r+1]^2/Pochhammer[mu/2+n+2r+1/2,n-r]^2/(mu+2n+1)/(2n+1)/(2n+2);
Table[Factor[(Dst[2r,1,2n+2]/Dst[2r,1,2n])/quoC],{n,6},{r,n}]

TableForm[Table[Together[Dst[2r,1,2n]/(Pochhammer[mu+2r-1,2n]/(2n)!)],{r,5},{n,r}]]

Table[Together[Dst[2r,1,2r]/((mu-1)*Pochhammer[mu+2r,2r-1]/(2r)!)],{r,5}]

Table[Together[Dst[2r,1,2n]/
If[n<r,Pochhammer[mu+2r-1,2n]/(2n)!,(mu-1)*Pochhammer[mu+2r,2r-1]/(2r)!*Product[quoC,{n,r,n-1}]]],{n,5},{r,5}]

Table[Together[Dst[1,2r,n]/(Product[Pochhammer[mu+i,2r-1]/Pochhammer[i+2,2r-1],{i,0,n-1}]*Dst[2r,1,n])],{n,8},{r,8}]

Conjecture 21  (Family D):  (Subscript[D, -1,2r])(n)

quoD=-Pochhammer[mu+2n-3,2r+1]*Pochhammer[mu+2n-1,2r]*Pochhammer[mu+2n+4r,n-r]^2*Pochhammer[mu/2+2n+r+1/2,n-r-1]^2/Pochhammer[n-r,n-r]^2/Pochhammer[2n+2,2r+1]/Pochhammer[2n+1,2r]/Pochhammer[mu/2+n+2r+1/2,n-r-1]^2;
Table[Together[(Dst[-1,2r,2n+2]/Dst[-1,2r,2n])/quoD],{n,5},{r,0,n-1}]

Table[Together[Dst[-1,2r,2n]/(Product[quoD,{n,r+1,n-1}]*Dst[-1,2r,2r+2])],{n,5},{r,0,n-1}]

Table[Together[(Dst[-1,2r,2n+2]/Dst[-1,2r,2n])/
(Pochhammer[mu+2n-2,2r]*Pochhammer[mu+2n-1,2r]/Pochhammer[2n+1,2r]/Pochhammer[2n+2,2r])],{r,0,5},{n,r-1}]

Table[Together[Dst[-1,2r,2r]/(Product[Pochhammer[mu+2i-2,2r]*Pochhammer[mu+2i-1,2r]/Pochhammer[2i+1,2r]/Pochhammer[2i+2,2r],{i,0,r-1}])],{r,4}]

Table[Together[Dst[-1,2r,2n]/
(Product[quoD,{n,r+1,n-1}]*(Dst[-1,2r,2r+2]/Dst[-1,2r,2r])*Product[Pochhammer[mu+2i-2,2r]*Pochhammer[mu+2i-1,2r]/Pochhammer[2i+1,2r]/Pochhammer[2i+2,2r],{i,0,r-1}])],{n,5},{r,1,n-1}]

Table[Together[(Dst[-1,2r,2r+2]/Dst[-1,2r,2r])/((3-mu)*Pochhammer[mu+2r-2,2r]*Pochhammer[mu+2r-1,2r]/Pochhammer[2r+1,2r]/Pochhammer[2r+1,2r+1])],{r,5}]

RD[r_,n_]:=Which[
n>r,-Pochhammer[mu+2n-3,2r+1]*Pochhammer[mu+2n-1,2r]*Pochhammer[mu+2n+4r,n-r]^2*Pochhammer[mu/2+2n+r+1/2,n-r-1]^2/Pochhammer[n-r,n-r]^2/Pochhammer[2n+2,2r+1]/Pochhammer[2n+1,2r]/Pochhammer[mu/2+n+2r+1/2,n-r-1]^2 (* this is quoD *),
n==r,(3-mu)*Pochhammer[mu+2r-2,2r]*Pochhammer[mu+2r-1,2r]/Pochhammer[2r+1,2r]/Pochhammer[2r+1,2r+1],
n<r,Pochhammer[mu+2n-2,2r]*Pochhammer[mu+2n-1,2r]/Pochhammer[2n+1,2r]/Pochhammer[2n+2,2r]];
Table[Together[Dst[-1,2r,2n]/Product[RD[r,i],{i,0,n-1}]],{n,5},{r,0,10}]

Corollary 22  (Family E):  (Subscript[D, 2r-1,1])(2n)/(Subscript[D, 2r-1,1])(2n-1)

quoE=Pochhammer[mu/2+2n+r-1/2,n-r]* Pochhammer[mu+2n+4r-4,n-r+1]/Pochhammer[n-r+1,n-r+1] /Pochhammer[mu/2+n+2r-3/2,n-r];
Table[Together[(Dst[2r-1,1,2n]/Dst[2r-1,1,2n-1])/quoE],{n,5},{r,n}]

Table[Together[(Dst[1,2r-1,2n]/Dst[1,2r-1,2n-1])/(Pochhammer[mu+2 n-1,2 r-2]/Pochhammer[2n+1,2r-2]*quoE)],{n,5},{r,n}]

TableForm[Table[Factor[Dst[2r-1,1,n+1]/Dst[2r-1,1,n]],{r,4},{n,2r}]]

Corollary 23  (Family F):  (Subscript[D, -1,2r-1])(2n+1)/(Subscript[D, -1,2r-1])(2n)

quoF=2Pochhammer[-2+2 n+mu,2r]*Pochhammer[-2+4 r+2 n+mu,n-r-1]*Pochhammer[mu/2+2n+r-1/2,n-r]/Pochhammer[2n,2r]/Pochhammer[n-r+1,n-r]/Pochhammer[mu/2+n+2r-1/2,n-r-1];
Table[Factor[(Dst[-1,2r-1,2n+1]/Dst[-1,2r-1,2n])/quoF],{n,4},{r,n}]

Encore: Guessing another closed form for (Subscript[D, 1,2r-1])(2n)/(Subscript[D, 1,2r-1])(2n-1)

dets=Table[Dst[1,2r-1,n,100003],{n,40},{r,20}];
data=Table[dets[[n+1,r]]/dets[[n,r]],{n,1,39,2},{r,20}];
guess=GuessMultRE[data,{f[m,r],f[m+1,r],f[m,r+1]},{m,r},6,StartPoint->{1,1},Constraints->m>r];
gb=Factor[OreGroebnerBasis[NormalizeCoefficients/@ToOrePolynomial[guess,f[m,r]]]]

Reconst[x_,v_]:=With[{a=Round[v/x]},(mu+x*a-v)/a];
Factor[gb=NormalizeCoefficients/@Expand[gb/.x_Integer/;Abs[x]>10000:>Reconst[x,100003]]]

sol=RSolve[ApplyOreOperator[gb[[2]],f[m]]==0,f[m],m][[1,1,2]];
Factor[gb=DFiniteTimes[Annihilator[1/sol,{S[m],S[r]}],gb]]

sol=Simplify[sol*RSolve[ApplyOreOperator[gb[[1]],f[r]]==0,f[r],r][[1,1,2]]]

Union[Flatten[Table[(sol/.mu->100003)/C[1]^2/data[[m,r]],{r,6},{m,r,6}]]]

sol=sol/C[1]^2/(1/2);

sol=FullSimplify[sol,Element[{r,m},Integers]&&r≥1&&m≥1&&m≥r]

quo=sol//.Gamma[a_]/Gamma[b_]/;IntegerQ[Expand[a-b-m]]:>Pochhammer[b,Expand[a-b]]/.
Pochhammer[(mu-1)/6,m] Pochhammer[(1+mu)/6,m] Pochhammer[(3+mu)/6,m]->Pochhammer[mu/2-1/2,3m]/3^(3m)/.
Gamma[a_]:>Gamma[Expand[a]]/.
Gamma[a_+Rational[b_,2]]:>With[{z=a+b/2-1/2},2^Expand[1-2z]*Sqrt[Pi]*Gamma[Expand[2z]]/Gamma[z]]/.
Gamma[m-r+1]/Gamma[2m-2r+2]->1/Pochhammer[m-r+1,m-r+1]/.
Gamma[mu/2-1+r+2 m]/Gamma[mu/2-2+2 r+m]->Pochhammer[mu/2+2r+m-2,m-r+1]/.
Gamma[mu/2+2m+r-2]/Gamma[mu/2-1]->Pochhammer[mu/2-1,2m+r-1]/.
Gamma[-3+3 r+3 m+mu]/Gamma[-2+2 r+4 m+mu]->1/Pochhammer[-3+3 r+3 m+mu,1-r+m]/.
Gamma[mu-2]/Gamma[mu+2m-1]->1/Pochhammer[mu-2,2m+1]/.
Gamma[mu+2m+2r-3]/Gamma[mu+4m+2r-4]->1/Pochhammer[mu+2m+2r-3,2m-1]/.
(2m)!/(2(m+r-1))!->1/Pochhammer[2m+1,2r-2]/.
Pochhammer[mu/2-1/2,3m]->Pochhammer[mu/2-1/2,2m+r-1]*Pochhammer[mu/2-1/2+2m+r-1,m-r+1]/.
Pochhammer[mu/2-1,2m+r-1]*Pochhammer[mu/2-1/2,2m+r-1]->Pochhammer[mu-2,4m+2r-2]/2^(4m+2r-2)/.
Pochhammer[mu-2,4m+2r-2]/Pochhammer[mu-2,2m+1]->Pochhammer[mu+2m-1,4m+2r-2m-3]/.
Pochhammer[mu+2m-1,2m+2r-3]->Pochhammer[mu+2m-1,2r-2]*Pochhammer[mu+2m+2r-3,2m-1]

Table[Det[DstMat[1,2r-1,2m,37]]/Det[DstMat[1,2r-1,2m-1,37]]/(quo/.{mu->37}),{r,10},{m,r,10}]

TraditionalForm[quo/.{mu->μ}]