Some remarks on the usage of this notebook: 
- The package HolonomicFunctions can be downloaded from the website
    http://www.risc.uni-linz.ac.at/research/combinat/software/HolonomicFunctions/
- Throughout this notebook, the variable qn is identified with the expression q^n (and analogously for qi and qj).
- The q-shift which sends qn to q*qn and q^n to q^(n+1) is denoted by QS[qn,q^n].

The objective of this notebook is to give a fully rigorous computer proof of the long-standing q-TSPP conjecture, which has been reduced to a determinant evaluation by Soichi Okada, which can been further reduced by the holonomic ansatz to proving three q-holonomic identities:

Normalization:	c(n,n)=1			(n≥1)
Soichi:		(∑_(j=1)^n)a(i,j) c(n,j)=0		(1≤i<n)
Okada:		(∑_(j=1)^n)a(n,j) c(n,j)==((b(n))/(b(n-1)))	(n≥1)

where  a(i,j)=((q^i)+1) (KroneckerDeltaSeq[i,j])-(KroneckerDeltaSeq[i,j+1])+(q^(i+j-1)) ((QBinomial[i+j-2, i-1, q])+q (QBinomial[i+j-1, i, q])) ,
and  b(n)==(∏_(i=1)^n)(∏_(j=i)^n)(∏_(k=j)^n)(((((1-(q^(i+j+k-1)))/(1-(q^(i+j+k-2))))))^2),
and c(n,j) is another discrete function to be defined below.

<<"HolonomicFunctions.m"

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

(* define the entries of the matrix *)
Clear[a];
a[i_Integer,j_Integer]:=a[i,j]=
Expand[FunctionExpand[q^(i+j-1)*(QBinomial[i+j-2,i-1,q]+q*QBinomial[i+j-1,i,q])+(1+q^i)*KroneckerDelta[i,j]-KroneckerDelta[i,j+1]]];

For convenience, we define a simple function that gives all non-negative integer roots of a "univariate" q-polynomial. By that we mean a polynomial in q and some other variable that represents a power of q, e.g., q^n. The output then is a list of values n ≥ 0 such that the given polynomial is zero.

MyQRoots[poly_,var_]:=
With[{factors=First/@FactorList[poly]},
If[MemberQ[factors,0],
"Polynomial is identically zero.",
Sort[Exponent[#,q]&/@Select[factors,MatchQ[#,(var-q^_.)|(q^_.-var)|(1-var)|(var-1)]&]]]]

Our definition of c(n,j)

We pull out of the hat the following annihilating ideal, and use it to define a bivariate sequence uniquely.

annc=<<"ann-qTSPP-deg.m";

The question which initial values have to be given in order to define the sequence uniquely, is considered in the following. In particular we do a detailed analysis of the singularities. For this purpose we extract the leading coefficients of the Gröbner basis annc and shift them according to the leading power products (such that if one of these polynomials is zero for some specific ((Subscript[n, 0]),(Subscript[j, 0])) then the value of c((Subscript[n, 0]),(Subscript[j, 0])) cannot be computed by means of the corresponding recurrence).

lcfs=(LeadingCoefficient[#]/.Thread[{qn,qj}->{qn,qj}/q^LeadingExponent[#]])&/@annc;

We want to define c(n,j) for n ≥ 1 and j ≥ 0. The latter is just for convenience: when we are going to prove (Soichi) it will simplify things if we know that our definition implies c(n,0)=0  for all n. (to be shown later).
With these restrictions it is clear that in the area (4,5)+(ℕ^2)  (i.e.,  all points with j ≥ 4 and n ≥ 5) all recurrences may be used to compute values of c(n,j) without involving values from outside our definition range. The remaining points, two strips of width 4 along the coordinate axes, are investigated separately.

The strip 1 ≤ n ≤ 4

LeadingPowerProduct/@annc

It is clear that in the case n=1, we can only use the first recurrence, for n=2 either the first or the second one, and so on.

n = 1:

MyQRoots[#/.qn->q,qj]&/@Take[lcfs,1]

singn={1,#}&/@First[%]

n = 2:

MyQRoots[#/.qn->q^2,qj]&/@Take[lcfs,2]

singn=Join[singn,{{2,0},{2,1},{2,2},{2,3},{2,4}}];

n = 3:

MyQRoots[#/.qn->q^3,qj]&/@Take[lcfs,3]

singn=Join[singn,{{3,0},{3,1},{3,2},{3,3}}];

n = 4:

MyQRoots[#/.qn->q^4,qj]&/@Take[lcfs,4]

singn=Join[singn,{{4,0},{4,1},{4,2}}];

The strip 0 ≤ j ≤ 3

j = 0:

MyQRoots[#/.qj->1,qn]&/@Take[lcfs,-1]

singj={#,0}&/@Rest[First[%]]

j = 1:

MyQRoots[#/.qj->q,qn]&/@Take[lcfs,-2]

singj=Join[singj,{{1,1},{2,1},{3,1},{4,1},{5,1}}];

j = 2:

MyQRoots[#/.qj->q^2,qn]&/@Take[lcfs,-3]

singj=Join[singj,{{1,2},{2,2},{3,2},{4,2}}];

j = 3:

MyQRoots[#/.qj->q^3,qn]&/@Take[lcfs,-4]

singj=Join[singj,{{1,3},{2,3},{3,3}}];

the area (4,5)+(ℕ^2)

We show that the leading coefficient of the third recurrence does not have any integer zeros in the area under consideration. A fortiori there cannot be any common zeros.

lcf=lcfs[[3]];

A numerical test before we start to be rigorous:

test=lcf/.q->2;
Position[Table[test/.qn->2^n/.qj->2^j,{j,4,20},{n,5,20}],0,{2}]

First we are extracting all the exponents that appear in the expanded leading coefficient. For each choice of n and j, we pick only the one with the highest power in q.

exps=DeleteCases[Flatten[MapIndexed[(Exponent[#1,q]+{n,j}.(#2-1))&,CoefficientList[fac,{qn,qj}],{2}]],-Infinity];

In this picture we see all exponents of q^n (horizontally) and q^j (vertically).

ListPlot[Coefficient[#,{n,j}]&/@exps]

Now we sort out all exponents that are always smaller than some other (again, inside the area under consideration). What we get is the upper right border of the above figure (this was somehow to be expected).

i1=1;
While[i1≤Length[exps],
i2=1;
While[i2≤Length[exps]&&Simplify[exps[[i1]]<exps[[i2]],n≥5&&j≥4]=!=True,i2++];
If[i2>Length[exps],i1++,exps=Drop[exps,{i1}]];
]

exps

Now we compute conditions on n and j, under which each of the above exponents becomes the maximal one.

Table[{k," ",Simplify[Reduce[Simplify[n≥5&&j≥4&&And@@(exps[[k]]>#&/@Drop[exps,{k}])]],Element[{n,j},Integers]]},{k,Length[exps]}]//TableForm

The entry at position 7 simplifies to n ≥ 5 && n-1 ≤ j ≤ n.  
So we see that all points are covered, except the lines j = 2n+1, j = n+1, j = n-2, and j = n-3, as it can also be seen in the following picture:

ArrayPlot[Table[With[{pos=Position[exps,Max[exps]]},If[Length[pos]===1,pos[[1,1]],0]],{n,100,5,-1},{j,4,100}],ColorRules->{0->White,1->Yellow,3->Red,5->Green,7->Blue,9->Black}]

For these exceptional lines, we test separately whether the leading coefficient becomes zero  (and find that none of the zeros lies in the area under consideration).

MyQRoots[lcf/.qj->qn^2*q,qn]

MyQRoots[lcf/.qj->qn*q,qn]

MyQRoots[lcf/.qn->qj*q^2,qj]

MyQRoots[lcf/.qn->qj*q^3,qj]

Conclusion

We have seen that there are no singularities for n ≥ 5 and j ≥ 4. 
Therefore the only singularities are the ones that we have found in the strips 1 ≤ n ≤ 4 and 0 ≤ j ≤ 3.

sing=Union[singn,singj]

We can hence write down a procedure that computes the value of c(n,j) for given integers n and j, using our annihilating ideal plus the above initial values (although this is not the most efficient way).

Clear[c];
localAnnc=OrePolynomialSubstitute[annc,{qn->q^n,qj->q^j}]/.{n->ln,j->lj};
c[nn_Integer/;nn≥1,jj_Integer/;jj≥0]:=c[nn,jj]=
Which[
(* special initial value (j=n=1) *)
jj===nn===1,1,
(* special initial values (j=0) *)
jj===0&&1≤nn≤4,0,
(* special initial values (j>n) *)
MemberQ[{{1,2},{1,3},{1,4},{1,5},{1,6},{2,3},{2,4}},{nn,jj}],0,
(* initial values that need to be computed *)
MemberQ[{{2,1},{2,2},{3,1},{3,2},{3,3},{4,1},{4,2},{5,1}},{nn,jj}],
Together[(-1)^(nn+jj)*Det[Table[Drop[Table[a[i1,j1],{j1,nn}],{jj}],{i1,nn-1}]]/Det[Table[a[i1,j1],{i1,nn-1},{j1,nn-1}]]],
(* use recurrences of annc *)
True,
Together[Solve[
Together[(0==ApplyOreOperator[#,c1[ln,lj]]/.Thread[{ln,lj}->({nn,jj}-LeadingExponent[#])])&/@
Select[localAnnc,Function[op,With[{lex=LeadingExponent[op]},FreeQ[(LeadingExponent[#]-lex+{nn,jj})&/@Support[op],{n1_,j1_}/;(n1<1||j1<0)]]]]/.c1[nn,jj]->x/.c1->c],x][[1,1,2]]]
];

table=Table[c[n,j],{n,9},{j,0,9}];

Map[If[#===0,0,"*"]&,table,{2}]//TableForm

Proof of (Normalization)

In order to get a recurrence that is valid on the diagonal j=n, we find an operator whose support consists of power products with QS[qj,q^j] and QS[qn,q^n] having the same power.

(* requires 3GB of memory *)
Timing[recDiag=First[FindRelation[annc,OrePolynomialSubstitute->{qj->qn},Support->Table[QS[qj,q^j]^i*QS[qn,q^n]^i,{i,0,7}]]];]

recDiag=NormalizeCoefficients[recDiag];

Omitting the shifts in j now gives a univariate recurrence of order 7 that is valid on the diagonal j=n.

recDiag=OrePolynomialSubstitute[recDiag,{QS[qj,q^j]->1},Algebra->OreAlgebra[QS[qn,q^n]]];

Exponent[OrePolynomialListCoefficients[recDiag],qn]

Exponent[OrePolynomialListCoefficients[recDiag],q]

ByteCount[recDiag]

We find that the diagonal recurrence is a multiple of the recurrence of the constant sequence 1.

OreReduce[recDiag,Annihilator[1,QS[qn,q^n]]]

We check whether c(n,n) = 1 for 1 ≤ n ≤ 7 using our definition of c(n,j):

Table[c[n,n],{n,7}]

Select[First/@FactorList[LeadingCoefficient[First[recDiag]]],Exponent[#,qn]===1&]

Together with the fact that the leading coefficient does not have any integer zeros, we have proven that c(n,n)=1 for all n ≥ 1.

Show that we have natural boundaries

We show that our definition of c(n,j) implies c(n,0)=0 and c(n,j)=0 for j > n. These facts will be necessary for the subsequent reasoning, in particular for treating the sums.

The initial values that are required for the following arguments have already been computed (and shown to be 0) when the procedure c[n,j] was defined.

Consider the case j > n. Having a look at the support of the first element in annc

Support[annc[[1]]]

reveals that we need to prove c(n,j)=0 for the 6 diagonals j = n+1, ..., j = n+6. It then follows thate c(n,j)=0 for all j > n by applying this recurrence.

The diagonals are treated in the same way as in the proof of (Normalization).

(* a little procedure to plot the monomials in the support of an operator. *)
PlotSupport[p_OrePolynomial]:=
Module[{points,m1,m2},
points=Reverse[LeadingExponent[#]]&/@Support[p];
{m1,m2}=(Max/@Transpose[points])+1;
Return[ListPlot[points,
ImageSize->Round[m1*50],
AspectRatio->m2/m1,
AxesLabel->ToOrePolynomial[Reverse[OreAlgebraGenerators[OreAlgebra[p]]]],
PlotMarkers->Style["×",Bold,15],
PlotRange->{{-1/2,m1},{-1/2,m2}},
Ticks->{Range[m1]-1,Range[m2]-1}]];
];

The special case c(n, 0) = 0

We first verify that applying the appropriate recurrence from our annihilating ideal implies c(n,0)=0.

LeadingPowerProduct/@annc

rec=OrePolynomialSubstitute[annc[[5]],{qj->1}];

Support[rec]

MyQRoots[LeadingCoefficient[rec],qn]

The diagonal j = n+1

Timing[rec=NormalizeCoefficients[First[FindRelation[annc,
OrePolynomialSubstitute->{qj->qn*q},
Support->Table[QS[qj,q^j]^k*QS[qn,q^n]^k,{k,0,7}],
Modulus->0]]];]

PlotSupport[rec]

MyQRoots[LeadingCoefficient[rec],qn]

The diagonal j = n+2

For the diagonal j = n+2, we can already make use of the zeros on j = n+1. This means that we do not necessarily need an operator whose support lies on one diagonal, but on two diagonals. Therefore the total degree becomes smaller and the computation much easier.

Timing[rec=NormalizeCoefficients[First[FindRelation[annc,
OrePolynomialSubstitute->{qj->qn*q},Support->Flatten[Table[{QS[qj,q^j]^jj*QS[qn,q^n]^nn},{jj,nn,nn+1},{nn,0,3}]],
Modulus->0]]];]

PlotSupport[rec]

MyQRoots[LeadingCoefficient[rec],qn]

The diagonal j = n+3

For j = n+3, the support of the operator may lie on 3 diagonals, analogously as before.

Timing[rec=NormalizeCoefficients[First[FindRelation[annc,
OrePolynomialSubstitute->{qj->qn*q},
Support->(QS[qj,q^j]^#[[1]]*QS[qn,q^n]^#[[2]]&/@
{{0,0},{1,0},{2,0},{3,2},{4,2}}),
Modulus->0]]];]

PlotSupport[rec]

MyQRoots[LeadingCoefficient[rec],qn]

The diagonal j = n+4

Timing[rec=NormalizeCoefficients[First[FindRelation[annc,
OrePolynomialSubstitute->{qj->qn*q},Support->Flatten[Table[{QS[qj,q^j]^jj*QS[qn,q^n]^nn},{jj,nn,nn+3},{nn,0,1}]],
Modulus->0]]];]

PlotSupport[rec]

MyQRoots[LeadingCoefficient[rec],qn]

The diagonal j = n+5

Timing[rec=NormalizeCoefficients[First[FindRelation[annc,
OrePolynomialSubstitute->{qj->qn*q},
Support->Table[QS[qj,q^j]^k,{k,0,4}],
Modulus->0]]];]

PlotSupport[rec]

MyQRoots[LeadingCoefficient[rec],qn]

The diagonal j = n+6

Timing[rec=NormalizeCoefficients[Last[FindRelation[annc,
OrePolynomialSubstitute->{qj->qn*q},
Support->Table[QS[qj,q^j]^k,{k,0,5}],
Modulus->0]]];]

PlotSupport[rec]

MyQRoots[LeadingCoefficient[rec],qn]

Proof of (Soichi)

We have computed two creative telescoping relations for
    c'(n,i,j)=(q^(i+j-1)) ((QBinomial[i+j-2, i-1, q])+q (QBinomial[i+j-1, i, q]))c(n,j)  
whose sizes are 1.5GB and 7GB, respectively.

We can verify that these creative telescoping relations lie in the annihilating ideal of the above function (i.e., the summand in the (Soichi)-sum). So let's first compute this ideal:

(* these computations take about 90 minutes. *)
ann=ChangeOreAlgebra[Join[annc,Annihilator[1,QS[qi,q^i]]],OreAlgebra[QS[qi,q^i],QS[qj,q^j],QS[qn,q^n]]];
ann=DFiniteTimes[ann,Annihilator[(q^(i+j-1)*(-1-q+q^i+q^(i+j))*QBinomial[i+j-2,i-1,q])/(q^i-1),{QS[qi,q^i],QS[qj,q^j],QS[qn,q^n]}]];

Singularity analysis of the annihilating ideal for the summand

We show that there are only finitely many singular points, which are additionally located close to the origin.

OreAlgebra[ann]

lcfs=Numerator[Together[(LeadingCoefficient[#]/.Thread[{qi,qj,qn}->{qi,qj,qn}/q^LeadingExponent[#]])&/@ann]];

LeadingPowerProduct/@ann

Whenever i>1, we can apply the first recurrence since their only singular points are i=0 and i=j=1:

Factor[lcfs[[1]]]

From now on we consider only the case i=1.

On the lines  n = 1, 2, 3, 4 we can apply the last recurrence (except finitely many points near the origin):

MyQRoots[lcfs[[6]]/.qi->q/.qn->q,qj]

MyQRoots[lcfs[[6]]/.qi->q/.qn->q^2,qj]

MyQRoots[lcfs[[6]]/.qi->q/.qn->q^3,qj]

MyQRoots[lcfs[[6]]/.qi->q/.qn->q^4,qj]

On the lines  j = 1, 2, 3, 4 we can apply the second recurrence (except finitely many points near the origin):

MyQRoots[lcfs[[2]]/.qi->q/.qj->q,qn]

MyQRoots[lcfs[[2]]/.qi->q/.qj->q^2,qn]

MyQRoots[lcfs[[2]]/.qi->q/.qj->q^3,qn]

MyQRoots[lcfs[[2]]/.qi->q/.qj->q^4,qn]

In the area i = 1 and j > 4 and n > 4, we can safely apply recurrence no. 4 (no singularities here):

exps=DeleteCases[Flatten[MapIndexed[(Exponent[#1,q]+{n,j}.(#2-1))&,CoefficientList[lcfs[[4]]/.qi->q,{qn,qj}],{2}]],-Infinity];

i1=1;
While[i1≤Length[exps],
i2=1;
While[i2≤Length[exps]&&Simplify[exps[[i1]]<exps[[i2]],n≥5&&j≥5]=!=True,i2++];
If[i2>Length[exps],i1++,exps=Drop[exps,{i1}]];
]

exps

Table[{k," ",Simplify[Reduce[Simplify[n≥5&&j≥5&&And@@(exps[[k]]>#&/@Drop[exps,{k}])]],Element[{n,j},Integers]]},{k,Length[exps]}]//TableForm

The cases j=n+1, j+2=n, j+3=n are not contained in the above regions and have to be checked separately:

MyQRoots[lcfs[[4]]/.qi->q/.qj->q*qn,qn]

MyQRoots[lcfs[[4]]/.qi->q/.qn->q^2*qj,qj]

MyQRoots[lcfs[[4]]/.qi->q/.qn->q^3*qj,qj]

Prove correctness of the certificates

Now we load and assemble the first creative telescoping relation.

ansatz=<<"ansatz1.m";
denom=<<"denom1.m";
solution=<<"solution1.m";
ctrel=ansatz/.myDen[a1__][a2_]:>(denom[[a1]]/.qj->a2)/.solution;

ByteCount[ctrel]

reduce=OreReduce[ctrel,ann];

We find that reduce is zero, proving that the first creative telescoping relation is correct. The reduction takes 12h and needs 14GB of memory.

A strong indication for the correctness is obtained in much less time by performing the reduction in a modular way. We give some examples in the following

Timing[OreReduce[ctrel,ann,OrePolynomialSubstitute->{q->37,qi->98,qj->-113,qn->-5}]]

Timing[OreReduce[ctrel,ann,OrePolynomialSubstitute->{q->-58,qi->20,qj->31}]]

Timing[OreReduce[ctrel,ann,OrePolynomialSubstitute->{qi->1089,qj->81,qn->6}]]

Timing[OreReduce[ctrel,ann,OrePolynomialSubstitute->{q->67,qn->133}]]

Timing[OreReduce[ctrel,ann,OrePolynomialSubstitute->{q->13}]]

Timing[OreReduce[ctrel,ann,OrePolynomialSubstitute->{qn->27}]]

We now assemble and reduce the second creative telescoping relation.

ansatz=<<"ansatz2.m";
denom=<<"denom2.m";
solution=<<"solution2.m";
ctrel=ansatz/.myDen[a1__][a2_]:>(denom[[a1]]/.qj->a2)/.solution;

ByteCount[ctrel]

reduce=OreReduce[ctrel,ann];

The reduction takes 65h, needs 37GB of memory, and finally yields reduce = 0.

Again, some modular computations that can be carried out more easily and that are nearly equally convincing.

Timing[OreReduce[ctrel,ann,OrePolynomialSubstitute->{q->3,qi->13,qj->-18,qn->47}]]

Timing[OreReduce[ctrel,ann,OrePolynomialSubstitute->{q->-17,qi->229,qn->42}]]

Problem statement

After simple rewriting, (Soichi) looks as follows: For 2 <= i < n

TraditionalForm[HoldForm[Sum[q^(i+j-1)*(QBinomial[i+j-2,i-1,q]+q*QBinomial[i+j-1,i,q])*c[n,j],{j,1,n}]==c[n,i-1]-(1+q^i)c[n,i]]]

and for i = 1

TraditionalForm[HoldForm[Sum[q^(i+j-1)*(QBinomial[i+j-2,i-1,q]+q*QBinomial[i+j-1,i,q])*c[n,j],{j,1,n}]/.i->1==-(1+q)c[n,1]]]

Obviously, we can view the case i=1 as a special case of the first formula since c(n,0)=0 holds (this has been proven before). And this is the reason why we have extended our definition of c(n,j) to j=0.

Next, observe that the sum runs over natural boundaries: we have that c'(n,i,j)=0 for j ≤ 0 (because of the binomial coefficients) and that c'(n,i,j}=0 for j>n (because c(n,j)=0 for j>n). This information is sufficient to deduce that when we sum over the creative telescoping relation, the delta part vanishes completely. Hence the above sum is annihilated by either of the two principal parts. We extract them now from the creative telescoping relation.

ansatz=<<"ansatz1.m";solution=<<"solution1.m";
pp1=OrePolynomialSubstitute[ansatz,{q[__]->0}]/.solution;

ByteCount[pp1]

Support[pp1]

ansatz=<<"ansatz2.m";solution=<<"solution2.m";
pp2=OrePolynomialSubstitute[ansatz,{q[__]->0}]/.solution;

ByteCount[pp2]

Support[pp2]

Determine singularities

In this section we want to show that the two recurrences pp1 and pp2 uniquely define a bivariate sequence, provided that some finite set of initial values is given. In particular we show that they do not have relevant singularities.

lcfs=LeadingCoefficient/@{pp1,pp2};

LeadingPowerProduct/@{pp1,pp2}

lcfs[[1]]=Numerator[Factor[lcfs[[1]]/.qi->qi/q^3]];

lcfs[[2]]=Numerator[Factor[lcfs[[2]]/.qn->qn/q^4]];

The leading coefficients consist of some small factors and one big factor each. We have a look at the small factors and treat the big ones separately.

Table[ByteCount/@(List@@lcfs[[k]]),{k,2}]

Most[lcfs[[1]]]

→ is zero for i=0,1,2,  and i=n, and  i=n+1.

Most[lcfs[[2]]]

→ is zero for n=1,2,3, and i=n, and i=n-1, and i=n+1, and (i,n)=(1,0),(0,1),(0,0).

f1=Last[lcfs[[1]]];f2=Last[lcfs[[2]]];

f1 has no zeros

We first test numerically that there are no zeros for 1 <= i < n <= 20:

test=f1/.q->2;
Position[Table[test/.qi->2^i/.qn->2^n,{n,20},{i,n-1}],0,{2}]

Exponent[f1,{qi,qn}]

In the following we will argue by case distinction that for any choice of integers i and n (1 <= i < n), there is one single term that yields the highest power of q which therefore cannot cancel away. This proves that f1 cannot be zero for those choices of i and n.

These are the exponents of q^i and q^n that occur in f1:

ListPlot[Position[CoefficientList[f1,{qi,qn}],Except[0],{2},Heads->False]-1]

Since we are interested in the highest q-powers only, we take for each power  product of q^i and q^n only the term with the highest power of q. In the following list we store the corresponding exponents. If we can prove that for any integers i and n there is a single maximal element in this list, we are done.

exps=DeleteCases[Flatten[
MapIndexed[(#1+Total[{i,n}*(#2-1)])&,
Map[Exponent[#,q]&,CoefficientList[f1,{qi,qn}],{2}],{2}]],-Infinity]

We sort out all elements that are obviously smaller than some other for all i>0 and n>0.

extremal=exps;
Function[exp1,c1=Flatten[PadRight[CoefficientList[exp1,{i,n}],{2,2}]];extremal=DeleteCases[extremal,exp2_/;With[{c2=Flatten[PadRight[CoefficientList[exp2,{i,n}],{2,2}]]},(And@@Thread[c2≤c1])&&Total[c2]<Total[c1]]];]/@exps;

extremal

ListPlot[Coefficient[#,{i,n}]&/@extremal]

We compare all terms with the same n-coefficient and figure out the minimal i0 such that the right-most term (i.e., the one with the largest i-coefficient) always takes the biggest value among these terms for all i > i0.

Do[
With[{row=Select[extremal,Coefficient[#,n]===nn&]/.n->0},Print["nn = ",nn,", i > ",Max[Solve[#==Last[row],i][[1,1,2]]&/@Most[row]]];]
,{nn,14,22}]

Provided that i>4 we can restrict our investigation to the following set of exponents:

extremal=Table[Last[Select[extremal,Coefficient[#,n]===nn&]],{nn,14,22}]

TableForm[Table[Position[extremal,Max[extremal]],{i,0,15},{n,0,15}]]

It is now easy to show that for n big enough the term q^31*(q^i)^9*(q^n)^22 will be the highest power of q in f1, for 4 < i < n.

FullSimplify[And@@((Last[extremal]>#)&/@Most[extremal]),4<i<n&&n≥9]

For n < 9 there is nothing to do, since these (finitely many; note that n < 9 implies i < 8) cases have been checked (numerically) at the beginning of this section. 
We check the cases i=1,2,3,4 separately. After plugging in the corresponding q-power for q^i, we are left with a bivariate polynomial in q and q^n.

MyQRoots[f1/.qi->q,qn]

MyQRoots[f1/.qi->q^2,qn]

MyQRoots[f1/.qi->q^3,qn]

MyQRoots[f1/.qi->q^4,qn]

f2 has zeros at (i,n)=(1,2),(1,3),(1,4),(2,4),(3,4),(3,5)

We do a modular check and find all candidates for zeros in the area 1 <= i < n <= 20:

test=f2/.q->2;
Reverse/@Position[Table[test/.qi->2^i/.qn->2^n,{n,20},{i,n-1}],0,{2}]

f2/.Thread[{qi,qn}->q^#]&/@%

Exponent[f2,{qi,qn}]

In the following we will argue by case distinction that for any choice of integers i and n (1 <= i < n), there is one single term that yields the highest power of q which therefore cannot cancel away. This proves that f2 cannot be zero for those choices of i and n.

These are the exponents of q^i and q^n that occur in f2:

ListPlot[Position[CoefficientList[f2,{qi,qn}],Except[0],{2},Heads->False]-1]

Since we are interested in the highest q-powers only, we take for each power  product of q^i and q^n only the term with the highest power of q. In the following list we store the corresponding exponents. If we can prove that for any integers i and n there is a single maximal element in this list, we are done.

exps=DeleteCases[Flatten[
MapIndexed[(#1+Total[{i,n}*(#2-1)])&,
Map[Exponent[#,q]&,CoefficientList[f2,{qi,qn}],{2}],{2}]],-Infinity]

We sort out all elements that are obviously smaller than some other for all i>0 and n>0.

extremal=exps;
Function[exp1,c1=Flatten[PadRight[CoefficientList[exp1,{i,n}],{2,2}]];extremal=DeleteCases[extremal,exp2_/;With[{c2=Flatten[PadRight[CoefficientList[exp2,{i,n}],{2,2}]]},(And@@Thread[c2≤c1])&&Total[c2]<Total[c1]]];]/@exps;

extremal

ListPlot[Coefficient[#,{i,n}]&/@extremal]

We compare all terms with the same n-coefficient and figure out the minimal i0 such that the right-most term (i.e., the one with the largest i-coefficient) always takes the biggest value among these terms for all i > i0.

Do[
With[{row=Select[extremal,Coefficient[#,n]===nn&]/.n->0},Print["nn = ",nn,", i > ",Max[Solve[#==Last[row],i][[1,1,2]]&/@Most[row]]];]
,{nn,0,22}]

Provided that i>3 we can restrict our investigation to the following set of exponents:

extremal=Table[Last[Select[extremal,Coefficient[#,n]===nn&]],{nn,0,22}]

TableForm[Table[Position[extremal,Max[extremal]],{i,0,15},{n,0,15}]]

It is now easy to show that for 3 < i <= n-7  the term q^9*(q^i)^9*(q^n)^22 will be the highest power of q in f2.

FullSimplify[And@@((Last[extremal]>#)&/@Most[extremal]),3<i≤n-7]

We check the cases i=1, 2, 3, and i=n-6, n-5,..., n-1 separately. After plugging in the corresponding q-power for q^i, we are left with a polynomial in q and q^n only.

MyQRoots[f2/.qi->q,qn]

→ gives zeros at (i,n)=(1,2),(1,3),(1,4)

MyQRoots[f2/.qi->q^2,qn]

→ gives zeros at (i,n)=(2,4)

MyQRoots[f2/.qi->q^3,qn]

→ gives zeros at (i,n)=(3,4),(3,5)

MyQRoots[f2/.qi->qn/q^6,qn]

MyQRoots[f2/.qi->qn/q^5,qn]

MyQRoots[f2/.qi->qn/q^4,qn]

MyQRoots[f2/.qi->qn/q^3,qn]

→ gives a zero at (i,n)=(1,4)

MyQRoots[f2/.qi->qn/q^2,qn]

→ gives zeros at (i,n)=(1,3),(2,4),(3,5)

MyQRoots[f2/.qi->qn/q,qn]

→ gives zeros at (i,n)=(1,2),(3,4)

Summarizing, the zeros of f2 are (1,2),(1,3),(1,4),(2,4),(3,4),(3,5), and these are exactly the ones found at the beginning.

Conclusion

Support[pp1]

Support[pp2]

Together with the constraints 1 <= i < n, these supports imply that we have to give initial values for (i=1 and 2 <= n <= 7),  (i=2 and 3 <= n <= 8),  (i=3 and 4 <= n <= 9).

For i <= 3, the first recurrence cannot be used. The leading coefficient of pp2 does not have zeros in this area that are not already part of the initial conditions.

On the lines i=n-1, i=n-2, ..., i=n-6, the second recurrence cannot be used. The leading coefficient of pp1 does not have zeros here.

There are no common zeros of the two leading coefficients that could be relevant for the remaining area.

The situation is depicted in the picture below. Red are the points where the first recurrence can be applied, the area where the second can be used is shown in green. The black points are the ones that have to be given as initial conditions.

Show[Graphics[Join[
{Line[{{0,0},{13,13}}],Text["i = n",{25/2,23/2}]},
{Opacity[1/3],Red,Polygon[{{4,5},{12,13},{4,13}}],Green,Polygon[{{1,8},{6,13},{1,13}}]},{Opacity[1],Red},
Flatten[Table[Disk[{x,y},1/10],{x,4,12},{y,x+1,13}]],
{Green},
Flatten[Table[Disk[{x,y},1/10],{x,3},{y,x+7,13}]],
{Brown},
Flatten[Table[Disk[{x,y},1/10],{x,4,6},{y,x+7,13}]],
{Black},
Flatten[Table[Disk[{x,y},1/10],{x,3},{y,x+1,x+6}]]
],Frame->True,FrameLabel->{i,n}]]

The right-hand side

We compute an annihilating ideal for the right-hand side and show that pp1 and pp2 are contained in it.

(* Timing for DFiniteSubstitute is 1h, for DFiniteOreAction it is 8h. *)
nm=$NormalizeMethod;$NormalizeMethod="gcd";  (* causes a speed-up in this instance *)
Timing[rhs=DFiniteOreAction[DFiniteSubstitute[annc,{j->i-1},Algebra->OreAlgebra[QS[qn,q^n],QS[qi,q^i]]],1-(1+qi)*QS[qi,q^i]];]

(* takes 63s. red1 = 0. *)
Timing[red1=OreReduce[pp1,rhs];]

(* takes 16h and 4GB of memory. red2 = 0.*)
Timing[red2=OreReduce[pp2,rhs];]
$NormalizeMethod=nm;

We compare initial values for (i=1 and 2 <= n <= 7),  (i=2 and 3 <= n <= 8),  (i=3 and 4 <= n <= 9), see the subsection on singularities above. Again, we do not treat the case i=1 separately, instead we take the degenerate right-hand side into account via c(n,0)=0.

Table[
Sum[FunctionExpand[q^(i+j-1)*(QBinomial[i+j-2,i-1,q]+q*QBinomial[i+j-1,i,q])]*c[n,j],{j,n}]==c[n,i-1]-(1+q^i)c[n,i],
{n,i+1,i+6},{i,3}]//Simplify

This completes the proof of (Soichi). Quod erat demonstrandum.

Proof of (Okada)

We rewrite the sum in (Okada) slightly:

TraditionalForm[HoldForm[Sum[q^(n+j-1)*(q^(n+j)+q^n-q-1)/(q^n-1)*QBinomial[n+j-2,n-1,q]*c[n,j],{j,1,n}]+(q^n+1)-c[n,n-1]==b[n]/b[n-1]]]

We start by computing an annihilating ideal for the summand:

anj=Annihilator[q^(n+j-1)*(q^(n+j)+q^n-q-1)/(q^n-1)*QBinomial[n+j-2,n-1,q],{QS[qn,q^n],QS[qj,q^j]}]

(* takes about 3 hours. *)
ann=DFiniteTimes[annc,anj];

Singularity analysis of the annihilating ideal for the summand

We show that there are only finitely many singular points, which are additionally located close to the origin.

lcfs=Numerator[Together[(LeadingCoefficient[#]/.Thread[{qn,qj}->{qn,qj}/q^LeadingExponent[#]])&/@ann]];

LeadingPowerProduct/@ann

On the lines  n = 1, 2, 3, 4 we can apply the first recurrence (except finitely many points near the origin):

MyQRoots[lcfs[[1]]/.qn->q,qj]

MyQRoots[lcfs[[1]]/.qn->q^2,qj]

MyQRoots[lcfs[[1]]/.qn->q^3,qj]

MyQRoots[lcfs[[1]]/.qn->q^4,qj]

On the lines  j = 1, 2, 3, 4 we can apply the last recurrence (except finitely many points near the origin):

MyQRoots[lcfs[[5]]/.qj->q,qn]

MyQRoots[lcfs[[5]]/.qj->q^2,qn]

MyQRoots[lcfs[[5]]/.qj->q^3,qn]

MyQRoots[lcfs[[5]]/.qj->q^4,qn]

In the area j > 4 and n > 4, we can safely apply recurrence 3 (no singularities here):

exps=DeleteCases[Flatten[MapIndexed[(Exponent[#1,q]+{n,j}.(#2-1))&,CoefficientList[lcfs[[3]],{qn,qj}],{2}]],-Infinity];

i1=1;
While[i1≤Length[exps],
i2=1;
While[i2≤Length[exps]&&Simplify[exps[[i1]]<exps[[i2]],n≥5&&j≥5]=!=True,i2++];
If[i2>Length[exps],i1++,exps=Drop[exps,{i1}]];
]

exps

Table[{k," ",Simplify[Reduce[Simplify[n≥5&&j≥5&&And@@(exps[[k]]>#&/@Drop[exps,{k}])]],Element[{n,j},Integers]]},{k,Length[exps]}]//TableForm

The cases j=n+1, j+2=n, j+3=n are not contained in the above regions and have to be checked separately:

MyQRoots[lcfs[[3]]/.qj->q*qn,qn]

MyQRoots[lcfs[[3]]/.qn->q^2*qj,qj]

MyQRoots[lcfs[[3]]/.qn->q^3*qj,qj]

Prove correctness of the certificate

We have computed a creative telescoping relation for  (q^(n+j-1)) ((QBinomial[n+j-2, n-1, q])+q (QBinomial[n+j-1, n, q]))c(n,j) . We convince ourselves that this operator indeed is in the corresponding annihilating ideal. We first load our solution and assemble the whole operator.

ansatz=<<"ansatz3.m";
denom=<<"denom3.m";
solution=<<"solution3.m";
ctrel=ansatz/.myDen[a1__][a2_]:>(denom[[a1]]/.qj->a2)/.solution;

This time, the reduction is more involved (more reduction steps have to be done). Our first attempt was interrupted by a computer crash, after one week, when it consumed about 33GB and was halfway through the reduction.

Therefore we thought of a different strategy that finally was successful. The idea is to avoid any summation that occurs during the reduction, by considering each term in a polynomial separately. Assume we start with the polynomial p; we perform one reduction step on LT(p), the leading term of p (and later also for each of the remaing terms). The result is in general again a polynomial p' consisting of several terms; we proceed as before by treating them all separately. Note that no addition of rational functions has to be done, only multiplications. Of course, this recursive branching yields an exponential number of terms in the end, and therefore is not a good idea in general, but in our example it is feasible: we obtain about 30 rational functions for each monomial under the stairs, which finally have to be added together. Since there are 10 monomials under the stairs, we have divided the problem into 10 smaller subproblems.
This strategy is implemented in the procedure MyReduce; it stores the intermediate results in the given directory.

MyReduce[poly_OrePolynomial,ann1_List,dir_String]:=Module[{ann,lann,lpps,t,algebra,order,uts,fcount,todo,pos,red},uts=UnderTheStaircase[LeadingExponent/@ann1];
ann=SortBy[ann1,ByteCount];
lann=Length[ann];
lpps=LeadingExponent/@ann;
algebra=Together[OreAlgebra[ann]];order=MonomialOrder[ann[[1]]];
t=OrePolynomial[{#},algebra,order]&/@poly[[1]];
fcount=1;todo={};

(* store all terms of the input in order to treat them later. *)
While[t=!={},Put[First[t],dir<>"myred"<>ToString[fcount]<>".m"];
If[Not[MemberQ[uts,LeadingExponent[First[t]]]],AppendTo[todo,fcount]];
fcount++;
t=Rest[t];
];

While[todo=!={},
t=Get[dir<>"myred"<>ToString[First[todo]]<>".m"];
todo=Rest[todo];
pos=1;
While[pos≤lann&&Not[And@@Thread[lpps[[pos]]≤LeadingExponent[t]]],pos++];
red=OreTimes[OrePolynomial[{{1,LeadingExponent[t]-lpps[[pos]]}},algebra,order],ann[[pos]]];
t=Together[-LeadingCoefficient[t]/LeadingCoefficient[red]];
t=OrePolynomial[{{Together[t*#[[1]]],#[[2]]}},algebra,order]&/@Rest[First[red]];
While[t=!={},Put[First[t],dir<>"myred"<>ToString[fcount]<>".m"];
If[Not[MemberQ[uts,LeadingExponent[First[t]]]],AppendTo[todo,fcount]];
fcount++;
t=Rest[t];
];
];
];

The following command produces 360 files (altogether about 5GB).

MyReduce[ctrel,ann,""]

Since the polynomials under consideration are partially factored, a tremendous speed-up can be achieved by using special purpose procedures for certain polynomial operations.

MyExponent[poly_,var_]:=
If[Head[poly]=!=Times,Exponent[poly,var],Total[Exponent[#,var]&/@(List@@poly)]];

(* returns {gcd(p11,p21), p11/gcd, p21/gcd}. *)
(* The strategy is to compute the gcd pairwise for all factors. *)
MyPolynomialGCD[p11_,p21_]:=
Module[{a1,a2,p1=p11,p2=p21,g,i,j,ex1,ex2,gcd,todo,vars},
vars=Variables[{p1,p2}];
p1=If[Head[p1]===Times,List@@p1,{p1}];
p2=If[Head[p2]===Times,List@@p2,{p2}];
gcd=Intersection[p1,p2];
p1=Complement[p1,gcd];
p2=Complement[p2,gcd];
gcd=Times@@gcd;
ex1=Total[MyExponent[#,vars]]&/@p1;
ex2=Total[MyExponent[#,vars]]&/@p2;
todo=SortBy[Flatten[Table[{i,j,ex1[[i]]+ex2[[j]]},{i,Length[p1]},{j,Length[p2]}],1],Last];
While[todo=!={},
{i,j}=Take[First[todo],2];
todo=Rest[todo];
g=PolynomialGCD[p1[[i]],p2[[j]]];
If[g=!=1,gcd*=g;
p1[[i]]=Together[p1[[i]]/g];p2[[j]]=Together[p2[[j]]/g];
If[p1[[i]]===1,todo=DeleteCases[todo,{i,_,_}]];
If[p2[[j]]===1,todo=DeleteCases[todo,{_,j,_}]];
ex1=Total[MyExponent[g,vars]];
todo=SortBy[todo/.{i,a1_,a2_}:>{i,a1,a2-ex1}/.{a1_,j,a2_}:>{a1,j,a2-ex1},Last];
];
];
Return[{gcd,Times@@p1,Times@@p2}];
];

(* takes a list of rational functions, and computes their sum in togethered form. *)
(* In the option Replace a list of substitutions can be given (for modular computations). *)
Options[TogetherTotal]={Replace->{{}},NormalizeVector->True};
TogetherTotal[pl1_List,opts:((_Rule|_RuleDelayed)...)]:=Module[{i,pl=pl1,replace,normalize,repl,pl0,polymod,pli,vars,subs,plmod,pairs,gcd,p1,p2,max,t,prime,result},{replace,normalize}={Replace,NormalizeVector}/.{opts}/.Options[TogetherTotal];

(* if not otherwise stated, we multiply by the common denominator, and remove content. *)
If[normalize===True,
gcd=PolynomialLCM@@Denominator[pl];
If[gcd=!=1,pl=Together[gcd*pl]];
gcd=Fold[PolynomialGCD,First[#],Rest[#]]&[SortBy[pl,ByteCount]];
If[gcd=!=1,pl=Together[pl/gcd]];
];

result={};
pl0=pl;
Do[
pl=pl0/.replace[[repl]];

(* By modular computations we try to find out which pairs of polynomials have the largest gcd (those will be treated first). *)
vars=Variables[pl];
subs=Table[NextPrime[19,3*i],{i,Length[vars]-1}];
prime=NextPrime[Developer`$MaxMachineInteger,-1];
polymod=Function[poly,If[Head[poly]===Times,PolynomialMod[#,prime]&/@poly,PolynomialMod[poly,prime]]];
plmod=polymod[pl/.Thread[Complement[vars,{#}]->subs]]&/@vars;

(* List of indices which entries in pl are valid. *)
pli=Range[Length[pl]];

pairs=Function[pair,t=Sum[Exponent[PolynomialGCD@@Append[plmod[[i,pair]],Modulus->prime],vars[[i]]],{i,Length[vars]}];
{pair,t}]/@Subsets[pli,{2}];

(* Now we add the polynomials pairwise: first remove the gcd, add the rest, and multiply again be the gcd. *)
While[pairs=!={},
max=Max[#[[2]]&/@pairs];
{p1,p2}=First[pairs[[Position[pairs,{_,max}][[1,1]]]]];
gcd=MyPolynomialGCD[pl[[p1]],pl[[p2]]];
pl[[p1]]=First[gcd]*Expand[Total[Rest[gcd]]];

(* Update modular data and list of pairs and indices. *)
Do[plmod[[i,p1]]=polymod[pl[[p1]]/.Thread[Complement[vars,{vars[[i]]}]->subs]];,{i,Length[vars]}];
pairs=DeleteCases[pairs,{{p1,_}|{_,p1}|{p2,_}|{_,p2},_}];
pli=Complement[pli,{p2}];

If[pl[[p1]]===0,
pli=Complement[pli,{p1}];
,pairs=Join[pairs,Function[pair,t=Sum[Exponent[PolynomialGCD@@Append[plmod[[i,pair]],Modulus->prime],vars[[i]]],{i,Length[vars]}];
{pair,t}]/@({p1,#}&/@Complement[pli,{p1}])];
];
];
AppendTo[result,Total[pl[[pli]]]];

,{repl,Length[replace]}]; (* end of Do *)

Clear[pl,pl0,subs];
Return[result];
];

Now we are equipped with everything that is needed. We demonstrate for one monomial (viz. QS[qn,q^n]^3) how to proceed. The remaining ones can be treated in the same fashion, although they take longer. For the bigger cases the modular approach is highly preferable.

rats={};
Do[
t=Get["myred"<>ToString[nn]<>".m"];
If[LeadingExponent[t]==={3,0},AppendTo[rats,LeadingCoefficient[t]]];
,{nn,360}];

Length[rats]

(* takes about 6 hours. *)
TogetherTotal[rats]

Alternatively, we can do it for several values of one of the variables (according to the degree w.r.t. this variable). 
We first clear denominators and remove content manually.

Timing[rats=Together[rats*(PolynomialLCM@@Denominator[rats])];]

Timing[rats=Together[rats/(Fold[PolynomialGCD,First[#],Rest[#]]&[SortBy[rats,ByteCount]])];]

Max[MyExponent[#,qj]&/@rats]

(* takes 3-4 min per case *)
TogetherTotal[rats,NormalizeVector->False,Replace->Table[{qj->k},{k,-33,34}]]

Conclude the proof

We now assemble only the principal part (recurrence of order 7) of this creative telescoping relation. For the same reasons as in (Soichi), the sum has natural boundaries and hence is annihilated by the principal part.

ansatz=<<"ansatz3.m";solution=<<"solution3.m";
pp=OrePolynomialSubstitute[ansatz,{q[__]->0}]/.myDen[a1__][a2_]:>(denom[[a1]]/.qj->a2)/.solution;

Support[pp]

ByteCount[pp]

In order to obtain a recurrence for the (Okada)-sum, we have to "add" (closure property) the recurrence for  (1+(q^n)) c(n,n)-c(n,n-1)==(1+(q^n))-c(n,n-1)

We first need a recurrence for c(n,n-1), i.e., a recurrence that is valid on the diagonal j=n-1. We find an operator whose support consists of power products with QS[qj,q^j] and QS[qn,q^n] having the same power, and in whose coefficients qj has been replaced by qn/q.

(* takes 10h and needs 3GB of memory. *)
recDiag1=FindRelation[annc,OrePolynomialSubstitute->{qj->qn/q},Support->Table[QS[qj,q^j]^i*QS[qn,q^n]^i,{i,0,7}]];
recDiag1=NormalizeCoefficients/@recDiag1;

Omitting the shifts in j now gives a univariate recurrence that is valid on the diagonal j=n-1.

recDiag1=OrePolynomialSubstitute[recDiag1,{QS[qj,q^j]->1},Algebra->OreAlgebra[QS[qn,q^n]]];

Now we perform the closure property "plus" and obtain a recurrence for the whole (Okada)-sum

TraditionalForm[HoldForm[Sum[a[n,j]c[n,j],{j,1,n}]==Sum[q^(n+j-1)*(QBinomial[n+j-2,n-1,q]+q*QBinomial[n+j-1,n,q])*c[n,j],{j,1,n}]+(1+q^n)-c[n,n-1]]]

(* This computation takes 108h and needs about 3GB of memory. *)
nm=$NormalizeMethod;$NormalizeMethod="gcd";
recPlus=DFinitePlus[{Expand[pp]},Annihilator[1+q^n,QS[qn,q^n]],recDiag1];
$NormalizeMethod=nm;

ByteCount[recPlus]

Support[recPlus]

Now, after we have computed a recurrence for the left-hand side, we care about the right-hand side.

b[n_Integer]:=Product[((1-q^(i+j+k-1))/(1-q^(i+j+k-2)))^2,{i,n},{j,i,n},{k,j,n}]

By simple hand calculations we found:

b[n]/b[n-1]==QPochhammer[q^(2n),q,n]^2/QPochhammer[q^n,q^2,n]^2

FunctionExpand[Table[%,{n,10}]]

Hence we get the following second-order recurrence for the right-hand side:

annb=Together[Factor[Annihilator[QPochhammer[q^(2n),q,n]^2/QPochhammer[q^n,q^2,n]^2,QS[qn,q^n]]]]

The following computation proves that recPlus is a left multiple of the recurrence for b[n]/b[n-1]:

Timing[reduce=OreReduce[First[recPlus],annb];]

reduce

We investigate the leading coefficient of the recurrence recPlus:

lcoeff=LeadingCoefficient[recPlus];

Select[First/@FactorList[lcoeff],Exponent[#,qn]===1&]

→ The recurrence has no singularities for positive n. Alternatively, this fact can be shown as follows (avoiding the costly factorization). Considering lcoeff as a polynomial in q^n, we find that the q-degree of the leading coefficient is maximal among the q-degrees of all coefficients. This proves that for positive n, the polynomial cannot become zero.

Coefficient[lcoeff,qn^Exponent[lcoeff,qn]]

Max[Exponent[#,q]&/@CoefficientList[lcoeff,qn]]

Hence we have to compare 12 initial values for n=1, ..., 12:

Table[Together[Sum[a[n,j]*c[n,j],{j,n}]-b[n]/b[n-1]],{n,12}]

This completes the proof of (Okada). Quod erat demonstrandum.