This Mathematica notebook accompanies the article "Zeilberger's Holonomic Ansatz for Pfaffians" by Masao Ishikawa and Christoph Koutschan (Proceedings of the ISSAC 2012, pp. 227-233).
The packages that are required to reproduce the computations can be downloaded freely from the following webpage:
http://www.risc.jku.at/research/combinat/software/

<<Guess.m
<<HolonomicFunctions.m
<<Hyper.m

HolonomicFunctions package by Christoph Koutschan, RISC-Linz, Version 1.6 (12.04.2012)
⟶ Type  ?HolonomicFunctions  for help

Motzkin Number Pfaffian (Theorem 2)

Evaluate the Pfaffian Pf(a(Subscript[(i,j), 1≤i,j≤2n])) where a(i,j)=(j-i)M(i+j-3) with M(n) denoting the Motzkin numbers: M(n)==(∑_(k=0)^n)(1/(k+1))(Binomial[n, 2 k])(Binomial[2k, k]). The evaluation of this Pfaffian was conjectured in the paper "Pfaffian decomposition and a Pfaffian analogue of q-Catalan Hankel determinants'' (Masao Ishikawa, Hiroyuki Tagawa, Jiang Zeng, arXiv:1009.2004), see Formula (6.3) there.

MotzkinNumber[0|1]:=1;MotzkinNumber[n_Integer/;n>1]:=MotzkinNumber[n]=(3*(n-1)*MotzkinNumber[n-2]+(2n+1)*MotzkinNumber[n-1])/(n+2);

Clear[a,aij];
a[i_Integer,j_Integer]:=(j-i)*MotzkinNumber[i+j-3];
aij=DFiniteTimes[Annihilator[j-i,{S[i],S[j]}],First[CreativeTelescoping[Binomial[i+j-3,2k]*CatalanNumber[k],S[k]-1,{S[i],S[j]}]]]

Now we guess an implicit description (linear recurrences) for the auxiliary function (Subscript[c, 2n,i]). We set the option AdditionalEquations to Infinity, which means that all data is used for guessing. Thus we can be sure that the resulting recurrences hold for all values in the given array. In other words, it doesn't make a difference whether we produce values for (Subscript[c, 2n,i]) from its definition (recurrences plus initial values) or from the data array.

Timing[
dim=30;data={};
Do[
matrix=Table[a[i,j],{j,2n-1},{i,2n-1}];
vec=First[LinSolveQ[matrix]];
vec=PadRight[vec/Last[vec],2*dim];
AppendTo[data,vec];
,{n,dim}];

(* The parameters of the following commands are the result of trial and error. *)
guess1=GuessMultRE[data,{f[n,i],f[n,i+1],f[n+1,i],f[n,i+2],f[n+1,i+1]},{n,i},3,StartPoint->{1,1},AdditionalEquations->Infinity];
guess2=GuessMultRE[data,{f[n,i],f[n,i+1],f[n+1,i],f[n,i+2],f[n+2,i]},{n,i},6,StartPoint->{1,1},AdditionalEquations->Infinity];
guess3=GuessMultRE[data,{f[n,i],f[n,i+1],f[n+1,i],f[n,i+2],f[n,i+3]},{n,i},4,StartPoint->{1,1},AdditionalEquations->Infinity];
c2ni=OreGroebnerBasis[ToOrePolynomial[First/@{guess1,guess2,guess3},f[n,i]]];
]

Factor[c2ni]

UnderTheStaircase[c2ni]

Factor[FGLM[c2ni,Lexicographic]]

The recurrences that generate the annihilating ideal c2ni do not have any singularities, and thus we have to specify only 4 initial values, corresponding to the monomials under the stairs:

AnnihilatorSingularities[c2ni,{1,1}]

Boundary conditions for (Subscript[c, 2n,i])

We show that (Subscript[c, 2n,i])=0 for i≤0 and for i≥2n. This will be useful for the summations later ("natural boundaries").

Factor[c2ni[[2]]]

Provided with the appropriate initial conditions ((Subscript[c, 2,0])=(Subscript[c, 4,0])=0), we see that this recurrence produces zeros on the line i=0, since the terms ((Subscript[S, i])^2) and (Subscript[S, i]) vanish. Similarly for i=-1, since the term ((Subscript[S, i])^2) still vanishes (provided that (Subscript[c, 2,-1])=(Subscript[c, 4,-1])=0). Because of these two zero rows, it is clear that everything beyond them (i.e., for i<-1) must be zero as well. The following computation shows that setting the initial conditions to 0 is compatible with the recurrences:

test=ApplyOreOperator[c2ni,c[n,i]];Union[Flatten[Table[test,{n,1,5},{i,-5,5}]/.c[_,i_?NonPositive]->0/.c[a__]:>data[[a]]]]

Now recall the leading coefficient and the support of the first defining recurrence of (Subscript[c, 2n,i]):

{Factor[LeadingCoefficient[c2ni[[1]]]],Support[c2ni[[1]]]}

Since this coefficient does not vanish for any integer point in the area n≥2 and i≥2n, we can use this recurrence to produce the values of (Subscript[c, 2n,i]) in this area. The support of this recurrence indicates that we need only to show that (Subscript[c, 2n,2n])=0.

Factor[LeadingCoefficient[First[DFiniteSubstitute[c2ni,{i->2n}]]]]

The leading coefficient having no relevant singularities and the initial values being zero (by construction) show this. With the same argument it is clear that also (Subscript[c, 2,i])=0.

Factor[LeadingCoefficient[First[DFiniteSubstitute[c2ni,{n->1}]]]]

Identity (1)

We compute an annihilating operator for (Subscript[c, 2n,2n-1]). Its leading coefficient has no nonnegative integer roots, and it has the operator ((Subscript[S, n])-1) as a right factor. Therefore it annihilates any constant sequence.

Factor[diag=DFiniteSubstitute[c2ni,{i->2n-1}]]

OreReduce[diag,Annihilator[1,S[n]]]

The four initial values are 1 by construction and therefore (Subscript[c, 2n,2n-1])=1 for all n.

Identity (2)

Compute a creative telescoping operator for the left-hand side.

alg=OreAlgebra[S[n],S[i],S[j]];
Timing[ByteCount[smnd=DFiniteTimes[ToOrePolynomial[Append[aij,S[n]-1],alg],ToOrePolynomial[Append[c2ni,S[j]-1],alg]]]]

Timing[fct=FindCreativeTelescoping[smnd,S[i]-1];]

Factor[First[fct]]

These two operators annihilate the left-hand-side expression of identity (2). If we consider them as the defining recurrences for a bivariate sequence (in n and j) then we need the following initial values to fill the area 1≤j<2n:

AnnihilatorSingularities[fct[[1]],{1,1},Assumptions->j<2n]

ReleaseHold[Hold[Sum[data[[n,i]]*a[i,j],{i,1,2n-1}]]/.(First/@%)]

This concludes the proof of (2).

Remark: note that the above reasoning is about the maximal possible area. If one tries to extend it further, the first step being j=2n, the second recurrence found above breaks down:

Collect[ApplyOreOperator[fct[[1,2]],f[n,j]]/.j->j-2,f[__],Factor]

AnnihilatorSingularities[fct[[1]],{1,1},Assumptions->j≤2n]

Indeed, the values for j=2n are nonzero as we demonstrate in the next section.

Identity (3)

ai2n=DFiniteSubstitute[aij,{j->2n},Algebra->OreAlgebra[S[n],S[i]]]

Compute recurrences for the summand of identity (3):

Timing[ByteCount[smnd=DFiniteTimes[c2ni,ai2n]]]

Compute a creative telescoping operator in order to deal with the summation in identity (3):

Timing[ByteCount/@(fct=FindCreativeTelescoping[smnd,S[i]-1])]

fct=<<"fct_6.3.m";

Check the correctness of the creative telescoping operator (show that it is a member of the left ideal generated by the recurrences of smnd):

{{principalPart},{{deltaPart}}}=fct;
Timing[OreReduce[principalPart+(S[i]-1)**deltaPart,smnd]]

Hence the ratio  r(n)=((b(2 n))/(b(2 n-2))) satisfies the following recurrence (actually from n=1 on, since the Pfaffian evaluation is also true for n=0):

rec=ApplyOreOperator[Factor[principalPart],r[n]]

Together with the initial values we get a closed form for this quotient.

Table[Sum[data[[n,i]]*a[i,2n],{i,1,2n-1}],{n,1,2}]

RSolve[{rec==0,r[1]==1,r[2]==5},r[n],n]

It follows that b(2 n)==(∏_(k=1)^n)((b(2 k))/(b(2 k-2)))==(∏_(k=1)^n)(4 k-3)==(∏_(k=0)^(n-1))(4 k+1).

Central Delannoy Number Pfaffian (Theorem 3)

Evaluate the Pfaffian Pf(a(Subscript[(i,j), 1≤i,j≤2n])) where a(i,j)=(j-i)D(i+j-3) with D(n) denoting the central Delannoy numbers: D(n)==(∑_(k=0)^n)(Binomial[n, k]) (Binomial[n+k, k]). The evaluation of this Pfaffian was conjectured in the paper "Pfaffian decomposition and a Pfaffian analogue of q-Catalan Hankel determinants'' (Masao Ishikawa, Hiroyuki Tagawa, Jiang Zeng, arXiv:1009.2004), see Formula (6.4) there.

CentralDelannoyNumber[0]:=1;
CentralDelannoyNumber[1]:=3;
CentralDelannoyNumber[n_Integer/;n>1]:=CentralDelannoyNumber[n]=(-(n-1)*CentralDelannoyNumber[n-2]+3*(2n-1)*CentralDelannoyNumber[n-1])/n;

Clear[a,aij];
a[i_Integer,j_Integer]:=(j-i)*CentralDelannoyNumber[i+j-3];
aij=Annihilator[(j-i)*Sum[Binomial[i+j-3,k]*Binomial[i+j-3+k,k],{k,0,i+j-3}],{S[i],S[j]}];
Factor[aij]

Now we guess an implicit description (linear recurrences) for the auxiliary function (Subscript[c, 2n,i]). We set the option AdditionalEquations to Infinity, which means that all data is used for guessing. Thus we can be sure that the resulting recurrences hold for all values in the given array. In other words, it doesn't make a difference whether we produce values for (Subscript[c, 2n,i]) from its definition (recurrences plus initial values) or from the data array.

Timing[
dim=30;data={};
Do[
matrix=Table[a[i,j],{j,2n-1},{i,2n-1}];
vec=First[LinSolveQ[matrix]];
vec=PadRight[vec/Last[vec],2*dim];
AppendTo[data,vec];
,{n,dim}];

(* The parameters of the following commands are the result of trial and error. *)
guess1=GuessMultRE[data,{f[n,i],f[n,i+1],f[n+1,i],f[n,i+2],f[n+1,i+1]},{n,i},4,StartPoint->{1,1},AdditionalEquations->Infinity];
guess2=GuessMultRE[data,{f[n,i],f[n,i+1],f[n+1,i],f[n,i+2],f[n+2,i]},{n,i},8,StartPoint->{1,1},AdditionalEquations->Infinity];
guess3=GuessMultRE[data,{f[n,i],f[n,i+1],f[n+1,i],f[n,i+2],f[n,i+3]},{n,i},5,StartPoint->{1,1},AdditionalEquations->Infinity];
c2ni=OreGroebnerBasis[ToOrePolynomial[First/@{guess1,guess2,guess3},f[n,i]]];
]

Factor[c2ni]

The recurrences become slightly nicer when we transform the Gröbner basis to the lexicographic monomial order (however, we will work with the original c2ni, which corresponds to degree lexicographic order):

Factor[FGLM[c2ni,Lexicographic]]

UnderTheStaircase[c2ni]

We determine the points for which the recurrences in c2ni cannot be applied. We need to include the initial conditions for those points explicitly.

AnnihilatorSingularities[c2ni,{1,1}]

Boundary conditions for (Subscript[c, 2n,i])

We show that (Subscript[c, 2n,i])=0 for i≤0 and for i≥2n. This will be useful for the summations later ("natural boundaries").

Factor[c2ni[[2]]]

Provided with the appropriate initial conditions ((Subscript[c, 2,0])=(Subscript[c, 4,0])=0), we see that this recurrence produces zeros on the line i=0, since the terms ((Subscript[S, i])^2) and (Subscript[S, i]) vanish. Similarly for i=-1, since the term ((Subscript[S, i])^2) still vanishes (provided that (Subscript[c, 2,-1])=(Subscript[c, 4,-1])=0). Because of these two zero rows, it is clear that everything beyond them (i.e., for i<-1) must be zero as well. The following computation shows that setting the initial conditions to 0 is compatible with the recurrences:

test=ApplyOreOperator[c2ni,c[n,i]];Union[Flatten[Table[test,{n,1,5},{i,-5,5}]/.c[_,i_?NonPositive]->0/.c[a__]:>data[[a]]]]

Now recall the leading coefficient and the support of the first defining recurrence of (Subscript[c, 2n,i]):

{Factor[LeadingCoefficient[c2ni[[1]]]],Support[c2ni[[1]]]}

Since this coefficient does not vanish for any integer point in the area n≥2 and i≥2n, we can use this recurrence to produce the values of (Subscript[c, 2n,i]) in this area. The support of this recurrence indicates that we need only to show that (Subscript[c, 2n,2n])=0.

Factor[LeadingCoefficient[First[DFiniteSubstitute[c2ni,{i->2n}]]]]

The leading coefficient having no relevant singularities and the initial values being zero (by construction) show this. With the same argument it is clear that also (Subscript[c, 2,i])=0.

Factor[LeadingCoefficient[First[DFiniteSubstitute[c2ni,{n->1}]]]]

Identity (1)

We compute an annihilating operator for (Subscript[c, 2n,2n-1]). Its leading coefficient has no nonnegative integer roots, and it has the operator ((Subscript[S, n])-1) as a right factor. Therefore it annihilates any constant sequence.

Factor[diag=DFiniteSubstitute[c2ni,{i->2n-1}]]

OreReduce[diag,Annihilator[1,S[n]]]

The four initial values are 1 by construction and therefore (Subscript[c, 2n,2n-1])=1 for all n.

Identity (2)

Compute a creative telescoping operator for the left-hand side.

alg=OreAlgebra[S[n],S[i],S[j]];
Timing[ByteCount[smnd=DFiniteTimes[ToOrePolynomial[Append[aij,S[n]-1],alg],ToOrePolynomial[Append[c2ni,S[j]-1],alg]]]]

Timing[ByteCount/@(fct=FindCreativeTelescoping[smnd,S[i]-1])]

Factor[First[fct]]

These two operators annihilate the left-hand-side expression of identity (2). If we consider them as the defining recurrences for a bivariate sequence (in n and j) then we need the following initial values to fill the area 1≤j<2n:

AnnihilatorSingularities[fct[[1]],{1,1},Assumptions->j<2n]

ReleaseHold[Hold[Sum[data[[n,i]]*a[i,j],{i,1,2n-1}]]/.(First/@%)]

This concludes the proof of (2).

Remark: note that the above reasoning is about the maximal possible area. If one tries to extend it further, the first step being j=2n, the second recurrence found above breaks down:

Collect[ApplyOreOperator[fct[[1,2]],f[n,j]]/.j->j-2,f[__],Factor]

AnnihilatorSingularities[fct[[1]],{1,1},Assumptions->j≤2n]

Indeed, the values for j=2n are nonzero as we demonstrate in the next section.

Identity (3)

ai2n=DFiniteSubstitute[aij,{j->2n},Algebra->OreAlgebra[S[n],S[i]]]

Compute recurrences for the summand of identity (3):

Timing[ByteCount[smnd=DFiniteTimes[c2ni,ai2n]]]

Compute a creative telescoping operator in order to deal with the summation in identity (3):

Timing[ByteCount/@(fct=FindCreativeTelescoping[smnd,S[i]-1])]

fct=<<"fct_6.4.m";

Check the correctness of the creative telescoping operator (show that it is a member of the left ideal generated by the recurrences of smnd):

{{principalPart},{{deltaPart}}}=fct;
Timing[OreReduce[principalPart+(S[i]-1)**deltaPart,smnd]]

Hence the ratio  r(n)=((b(2 n))/(b(2 n-2))) satisfies the following recurrence (for n≥2, since the evaluation is conjectured to be true for n≥1):

rec=ApplyOreOperator[Factor[principalPart],r[n]]

Indeed, the recurrence is valid for all n>1, but it does not hold for n=1:

Table[rec,{n,1,10}]/.r[n_]:>Sum[data[[n,i]]*a[i,2n],{i,1,2n-1}]

Together with the initial values we get a closed form for this quotient. We use Marko Petkovsek's implementation of his algorithm Hyper to find the hypergeometric solutions of this recurrence. There is only one, but it turns out that it is exactly the one that we are looking for (two initial values match!).

{hyp}=Hyper[rec,r[n]]

Table[Sum[data[[n,i]]*a[i,2n],{i,1,2n-1}],{n,1,10}]

RSolve[{r[n+1]/r[n]==hyp,r[2]==72,r[3]==1120/3},r[n],n]

Table[%[[1,1,2]],{n,1,10}]

The solution matches except for the first value. This means that (Subscript[r, 1])=1 and (Subscript[r, n])=(((2^(2 n-1)) (2 n-1) (4 n-5))/(2 n-3)) for n≥2.

It follows that b(2 n)==(∏_(k=1)^n)((b(2 k))/(b(2 k-2)))==(∏_(k=1)^n)(Subscript[r, k])=(∏_(k=1)^n)(((2^(2 k-1)) (2 k-1) (4 k-5))/(2 k-3))==(2^(n^2))(∏_(k=1)^n)(( (2 k-1) (4 k-5))/(2 k-3))=-(2^(n^2))(2n-1)(∏_(k=1)^n) (4 k-5)=(2^(n^2))(2n-1)(∏_(k=1)^(n-1)) (4 k-1).

Narayana Polynomial Pfaffian (Theorem 4)

Evaluate the Pfaffian Pf(a(Subscript[(i,j), 1≤i,j≤2n])) where a(i,j)=(j-i)N(i+j-2,x) with N(n,x) denoting the n-th Narayana polynomial: N(n,x)==(∑_(k=0)^n)(1/n)(Binomial[n, k]) (Binomial[n, k-1])(x^k). The evaluation of this Pfaffian was conjectured in the paper "Pfaffian decomposition and a Pfaffian analogue of q-Catalan Hankel determinants'' (Masao Ishikawa, Hiroyuki Tagawa, Jiang Zeng, arXiv:1009.2004), see Formula (6.6) there.

NarayanaN[1,x_]:=x;
NarayanaN[2,x_]:=x*(x+1);
NarayanaN[n_Integer/;n>2,x_]:=NarayanaN[n,x]=
Expand[(-(x-1)^2*(n-2)*NarayanaN[n-2,x]+(x+1)*(2n-1)*NarayanaN[n-1,x])/(n+1)];

Clear[a,aij];
a[i_Integer,j_Integer]:=(j-i)*NarayanaN[i+j-2,x];
aij=Annihilator[(j-i)*Sum[1/(i+j-2)*Binomial[i+j-2,k]*Binomial[i+j-2,k-1]*x^k,{k,0,i+j-2}],{S[i],S[j]}];
Factor[aij]

Now we guess an implicit description (linear recurrences) for the auxiliary function (Subscript[c, 2n,i]).

Timing[
dim=26;data={};
Do[
matrix=Table[a[i,j],{j,2n-1},{i,2n-1}];
vec=First[LinSolveUniv[matrix,x]];
vec=PadRight[vec/Last[vec],2*dim];
AppendTo[data,vec];
,{n,dim}];

(* The parameters of the following commands are the result of trial and error. *)
guess1=GuessMultRE[data,{f[n,i],f[n,i+1],f[n+1,i],f[n,i+2],f[n+1,i+1]},{n,i},3,StartPoint->{1,1},Constraints->2i≤n,AdditionalEquations->Infinity];
guess2=GuessMultRE[data,{f[n,i],f[n,i+1],f[n+1,i],f[n,i+2],f[n+2,i]},{n,i},6,StartPoint->{1,1},Constraints->2i≤n,AdditionalEquations->Infinity];
guess3=GuessMultRE[data,{f[n,i],f[n,i+1],f[n+1,i],f[n,i+2],f[n,i+3]},{n,i},4,StartPoint->{1,1},Constraints->2i≤n,AdditionalEquations->Infinity];
c2ni=OreGroebnerBasis[ToOrePolynomial[First/@{guess1,guess2,guess3},f[n,i]]];
]

c2ni=<<"c2ni_6.6.m";

Factor[c2ni]

We determine the (finitely many) points for which the recurrences in c2ni cannot be applied. There are for such points in general, and some more in the special case x=-1. So for the moment, let's assume that x≠-1.

AnnihilatorSingularities[c2ni,{1,1}]

Boundary conditions for (Subscript[c, 2n,i])

We show that (Subscript[c, 2n,i])=0 for i≤0 and for i≥2n. This will be useful for the summations later ("natural boundaries").

Factor[c2ni[[2]]]

Provided with the appropriate initial conditions ((Subscript[c, 2,0])=(Subscript[c, 4,0])=0), we see that this recurrence produces zeros on the line i=0, since the terms ((Subscript[S, i])^2) and (Subscript[S, i]) vanish. Similarly for i=-1, since the term ((Subscript[S, i])^2) still vanishes (provided that (Subscript[c, 2,-1])=(Subscript[c, 4,-1])=0). Because of these two zero rows, it is clear that everything beyond them (i.e., for i<-1) must be zero as well. The following computation shows that setting the initial conditions to 0 is compatible with the recurrences:

test=ApplyOreOperator[c2ni,c[n,i]];Union[Flatten[Expand[Table[test,{n,1,5},{i,-5,5}]/.c[_,i_?NonPositive]->0/.c[a__]:>data[[a]]]]]

Now recall the leading coefficient and the support of the first defining recurrence of (Subscript[c, 2n,i]):

{Factor[LeadingCoefficient[c2ni[[1]]]],Support[c2ni[[1]]]}

Since this coefficient does not vanish for any integer point in the area n≥2 and i≥2n (again assuming x≠-1), we can use this recurrence to produce the values of (Subscript[c, 2n,i]) in this area. The support of this recurrence indicates that we need only to show that (Subscript[c, 2n,2n])=0.

Factor[LeadingCoefficient[First[DFiniteSubstitute[c2ni,{i->2n}]]]]

CylindricalDecomposition[Implies[n>1&&x<-1,%>0],{n,x}]

The leading coefficient having no relevant singularities (assuming x<-1) and the initial values being zero (by construction) show this. With the same argument it is clear that also (Subscript[c, 2,i])=0.

Factor[LeadingCoefficient[First[DFiniteSubstitute[c2ni,{n->1}]]]]

Identity (1)

We compute an annihilating operator for (Subscript[c, 2n,2n-1]). Its leading coefficient has no nonnegative integer roots (for x<-1), and it has the operator ((Subscript[S, n])-1) as a right factor. Therefore it annihilates any constant sequence.

Timing[ByteCount[diag=First[DFiniteSubstitute[c2ni,{i->2n-1}]]]]

OreReduce[diag,Annihilator[1,S[n]]]

Factor[LeadingCoefficient[diag]]

CylindricalDecomposition[Implies[n>1&&x<-1,%>0],{n,x}]

The four initial values are 1 by construction and therefore (Subscript[c, 2n,2n-1])=1 for all n.

Table[data[[n,2n-1]],{n,4}]

Identity (2)

Compute a creative telescoping operator for the left-hand side.

alg=OreAlgebra[S[n],S[i],S[j]];
Timing[ByteCount[smnd=DFiniteTimes[ToOrePolynomial[Append[aij,S[n]-1],alg],ToOrePolynomial[Append[c2ni,S[j]-1],alg]]]]

Timing[ByteCount/@(fct=FindCreativeTelescoping[smnd,S[i]-1])]

Factor[First[fct]]

These two operators annihilate the left-hand-side expression of identity (2). If we consider them as the defining recurrences for a bivariate sequence (in n and j) then we need the following initial values to fill the area 1≤j<2n:

AnnihilatorSingularities[fct[[1]],{1,1},Assumptions->j<2n]

Expand[ReleaseHold[Hold[Sum[data[[n,i]]*a[i,j],{i,1,2n-1}]]/.(First/@%)]]

This concludes the proof of (2).

Remark: note that the above reasoning is about the maximal possible area. If one tries to extend it further, the first step being j=2n, the second recurrence found above breaks down:

Collect[ApplyOreOperator[fct[[1,2]],f[n,j]]/.j->j-2,f[__],Factor]

AnnihilatorSingularities[fct[[1]],{1,1},Assumptions->j≤2n]

Indeed, the values for j=2n are nonzero as we demonstrate in the next section.

Identity (3)

ai2n=DFiniteSubstitute[aij,{j->2n},Algebra->OreAlgebra[S[n],S[i]]]

Compute recurrences for the summand of identity (3):

Timing[ByteCount[smnd=DFiniteTimes[c2ni,ai2n]]]

Compute a creative telescoping operator in order to deal with the summation in identity (3). 
With evaluation/interpolation it took us 1955s to compute such an operator.

fct=<<"fct_6.6.m";

Check the correctness of the creative telescoping operator (show that it is a member of the left ideal generated by the recurrences of smnd):

{{principalPart},{{deltaPart}}}=fct;
Timing[OreReduce[principalPart+(S[i]-1)**deltaPart,smnd]]

Hence the ratio  r(n)=((b(2 n))/(b(2 n-2))) satisfies the following recurrence (for n≥1, since the evaluation of the Pfaffian holds for n≥0):

rec=ApplyOreOperator[Factor[principalPart],r[n]]

Together with the initial values we get a closed form for this quotient. We use Marko Petkovsek's implementation of his algorithm Hyper to find the hypergeometric solutions of this recurrence. There is only one, but it turns out that it is exactly the one that we are looking for (two initial values match!).

{hyp}=Hyper[rec,r[n]]

Table[Together[Sum[data[[n,i]]*a[i,2n],{i,1,2n-1}]],{n,1,10}]

RSolve[{r[n+1]/r[n]==hyp,r[1]==x,r[2]==5x^3},r[n],n]

Mathematica can also solve this recurrence directly, but this takes much longer:

Timing[RSolve[{rec==0,r[1]==x,r[2]==5x^3},r[n],n]]

Table[%[[2,1,1,2]],{n,1,10}]

It follows that b(2 n)==(∏_(k=1)^n)((b(2 k))/(b(2 k-2)))==(∏_(k=1)^n)(Subscript[r, k])=(∏_(k=1)^n)((4k-3) (x^(2k-1)))==(x^(n^2))(∏_(k=1)^n)(4 k-3)==(x^(n^2))(∏_(k=0)^(n-1))(4 k+1).

Final remark: in some steps we have assumed that x is a real number smaller than -1. Hence for now the evaluation is proven only for x<-1. But for specific n, the Pfaffian is a polynomial in x (of a certain degree), as well as the right-hand side evaluation. Thus their difference is a polynomial in x which is zero for all x<-1. By the fundamental theorem of algebra it follows that this polynomial is identically zero, and therefore the evaluation of the Pfaffian is true for all complex numbers x.