######################################################################
##LehmerDet.txt: Save this file as  LehmerDet.txt                    #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read `LehmerDet.txt` <Enter>                       #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Aug. 2018
 
print(`Created:  Aug. 2018`):
print(` This is LehmerDet.txt `):
print(`It is one of the packages that accompany the article `):
print(` D.H. Lehmer's Tridiagonal determinant`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):

print(`and also available from Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):
with(linalg):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: GP, GuessPol1, Lseq, LseqF, Rseq, RseqF, Qrn `):
 print(``):

else
ezra(args):
fi:

end:

ezraS:=proc()

if args=NULL then
 print(` The STORY procedures are: GuessLehmer1v `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: CXn, d,  GuessLehmer1, GuessPol, M, QXn`):

 print(` `):


elif nops([args])=1 and op(1,[args])=CXn then
print(`CXn(X,n,q): The explicit expression conjectured by Shalosh, and whose derivation and proof is described  in our paper`):
print(`for the Lehmer n by n determinant, d(X,n,q) alias QXn(X,n,q). Try: `):
print(`CXn(X,5,q); `):

elif nops([args])=1 and op(1,[args])=d then
print(`d(X,n,q): The Lehmer n by n determinant done directly.`):
print(`It is given in his paper "Combinatorial and Cyclotomic properties of certain triadiagonal matrices"`):
print(`Congressus Numerantium, v. 10, (1974), 53-74.  Try:`):
print(`d(X,5,q); `):

elif nops([args])=1 and op(1,[args])=GP then
print(`GP(m,n,q): the Gaussian polynomial G(m,n)(q). Try:`):
print(`GP(4,5,q);`):

elif nops([args])=1 and op(1,[args])=GuessLehmer1 then
print(`GuessLehmer1(X,R,q): guesses the first  R coefficients of X in the Lehmer n by n determinant.  Try:`):
print(`GuessLehmer1(X,5,q);`):

elif nops([args])=1 and op(1,[args])=GuessLehmer1v then
print(`GuessLehmer1v(X,R,q): verbose version of GuessLehmer1(X,R,q) (q.v.). Try: `):
print(`GuessLehmer1v(X,5,q);`):


elif nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(L,q,X,st): inputs a list L of polynomials in q, and a variable X, and a positive integer st. `):
print(` and outputs a polynomial in  X, `):
print(`of degree <=nops(L)-3-st, P(q,X) such`):
print(`that the n-th entry of L is P(q,q^n) for n>=st. If none found, it returns FAIL. Try: `):
print(`GuessPol([seq(GP(n1,2,q),n1=1..20)],q,X,1);`):

elif nops([args])=1 and op(1,[args])=GuessPol1 then
print(`GuessPol1(L,q,X,d,st): inputs a list L of polynomials in q, and a variable X, and a positive integer d`):
print(`and outputs a polynomial  P(q,X) in X, of degree d, such`):
print(`that the n-th entry of L is P(q,q^n) for n>=st. Try: `):
print(`GuessPol1([seq(GP(n1,2,q),n1=1..20)],q,X,2,1);`):

elif nops([args])=1 and op(1,[args])=Lseq then
print(`Lseq(N): the first N terms of the  infinite Lehmer determinant, given in  OEIS sequence A039924 . Done directly. Try: `):
print(`Lseq(30);`):


elif nops([args])=1 and op(1,[args])=LseqF then
print(`LseqF(N): the first N terms of the  infinite Lehmer determinant, given in  OEIS sequence A039924 . Done via the generating function. Try: `):
print(`LseqF(30);`):

elif nops([args])=1 and op(1,[args])=M then
print(`M(X,n,q): The Lehmer n by n matrix given in his paper "Combinatorial and Cyclotomic properties of certain triadiagonal matrices"`):
print(`Congressus Numerantium, v. 10, (1974), 53-74. To get it with his symbols, with n=5, Try:`):
print(`matrix(simplify(M(a^2,5,r^2),symbolic));`):


elif nops([args])=1 and op(1,[args])=Qrn then
print(`Qrn(r,n,q): QXn(q^r,n,q). `):
print(`Try: `):
print(`Qrn(1,5,q); `):

elif nops([args])=1 and op(1,[args])=QXn then
print(`QXn(X,n,q): The Lehmer n by n determinant done via the obvious recurrence described in Lehmer's paper`):
print(`by expanding about the last row. Eq. (1) in his paper. Try: `):
print(`QXn(X,5,q); `):

elif nops([args])=1 and op(1,[args])=Rseq then
print(`Rseq(N): the first N terms of the  reciprocal of infinite Lehmer determinant, given in  OEIS sequence A003316 Done directly.. Try: `):
print(`Rseq(20);`):

elif nops([args])=1 and op(1,[args])=RseqF then
print(`Rseq(N): the first N terms of the  reciprocal of infinite Lehmer determinant, given in  OEIS sequence A003316 . `):
print(` done via Lehmer's generating function. Try: `):
print(`RseqF(20);`):


else
print(`There is no ezra for`,args):
fi:
 
end:






#GP(m,n,q): the Gaussian polynomial G(m,n)(q)
GP:=proc(m,n,q) local i:
sort(expand(normal(mul(1-q^(n+i),i=1..m)/mul(1-q^i,i=1..m)))):
end:


#M(X,n,q): The Lehmer n by n matrix given in his paper "Combinatorial and Cyclotomic properties of certain triadiagonal matrices"
#Congressus Numerantium, v. 10, (1974), 53-74. To get it with his symbols, with n=5, Try:
#M(a,5,r);
M:=proc(X,n,q)  local i:

if n=1 then
  RETURN([[1]]):
fi:

[
[1,sqrt(X),0$(n-2)],
seq([0$(i-2),sqrt(X)*q^((i-2)/2),1,sqrt(X)*q^((i-1)/2),0$(n-i-1)],i=2..n-1),
[0$(n-2),sqrt(X)*q^((n-2)/2),1]
]:
end:

#d(X,n,q): the Lehmer determinant
d:=proc(X,n,q): expand(det(M(X,n,q))): end:


#Lseq(N): the first N terms of the  Lehmer determinant  sequence A039924
Lseq:=proc(N) local gu,x,i:
gu:=Qrn(1,N+2,x):
gu:=taylor(gu,x=0,N+1):
[seq(coeff(gu,x,i),i=1..N)]:
end:

#Rseq(N): the first N terms of the  Lehmer determinant sequence A003316
Rseq:=proc(N) local gu,x,i:
gu:=1/Qrn(1,N+2,x):
gu:=taylor(gu,x=0,N+1):
[seq(coeff(gu,x,i),i=1..N)]:
end:

#LseqF(N): Like Pseq(N) but via the  generating function
LseqF:=proc(N) local q,gu,k,i:
gu:=add((-1)^k*q^(k^2)/mul(1-q^i,i=1..k),k=0..trunc(sqrt(N))+1):
gu:=taylor(gu,q=0,N+1):
[seq(coeff(gu,q,i),i=1..N)]:
end:

#RseqF(N): the first N terms of the  Lehmer determinant sequence A003316 done via the generating function
RseqF:=proc(N) local gu,i,k,q:
gu:=1/add((-1)^k*q^(k^2)/mul(1-q^i,i=1..k),k=0..trunc(sqrt(N))+1):
gu:=taylor(gu,q=0,N+1):
[seq(coeff(gu,q,i),i=1..N)]:
end:

#Qrn(r,n,x): 
Qrn:=proc(r,n,q)  local X:
expand(subs(X=q^r,QXn(X,n,q))):

end:

#GuessPol1(L,q,X,d,st): inputs a list L of polynomials in q, and a variable X, and outputs a polynomial in q and X, of degree d, P(q,X) such
#that the n-th entry of L, starting at st, is P(q,q^n). Try:
#GuessPol1([seq(GF(n1,2,q),n1=1..20)],q,X,2,1);
GuessPol1:=proc(L,q,X,d,st) local a,eq,var,P,i,r:
if nops(L)-st+1<d+3 then
 print(`List too short`):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=0..d)}:
P:=add(a[i]*X^i,i=0..d):

eq:=expand({seq(subs(X=q^r,P)-L[r],r=st..d+st+2)}):

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

P:=subs(var,P):

if expand({seq(normal(subs(X=q^r,P)-L[r]),r=d+st+2..nops(L))})<>{0} then
 print(P, `did not work out`):
 RETURN(FAIL):
fi:

factor(P):

end:




#GuessPol(L,q,X,st): inputs a list L of polynomials in q, and a variable X, and outputs a polynomial in q and X, of degree<=nops(L)-4
#(if found) P(q,X), or FAIL, such
#that the n-th entry of L is P(q,q^n), starting at n=st. Try:
#GuessPol([seq(GF(n1,2,q),n1=1..20)],q,X,1);
GuessPol:=proc(L,q,X,st) local d,gu:
for d from 0 to nops(L)-3-st do
gu:=GuessPol1(L,q,X,d,st):
if gu<>FAIL then
  RETURN(gu):
fi:
od:
FAIL:
end:




#CXn(X,n,q): The guessed Lehmer polynomial QXn(X,n,q)
CXn:=proc(X,n,q) local i:
expand(add((-1)^i*X^i*q^(i*(i-1))*GP(i,n-2*i,q),i=0..trunc(n/2))):
end:


#QXn(X,n,q): The Lehmer n by n determinant via recurrence
QXn:=proc(X,n,q) 
option remember:
if n=1 then
1:
elif n=2 then
 1-X:
else
expand(QXn(X,n-1,q)-q^(n-2)*X*QXn(X,n-2,q)):
fi:
end:

#GuessLehmer1(X,R,q): guesses the first  R coefficients of X in the Lehmer n by n determinant.  Try:
#GuessLehmer1(X,5,q);
GuessLehmer1:=proc(X,R,q) local gu,r,n1,i:
gu:=[seq(QXn(X,n1,q),n1=1..4*R+10)]:
[seq(GuessPol([seq(coeff(gu[i],X,r),i=1..nops(gu))],q,X,2*r),r=1..R)]:
end:




#GuessLehmer1v(X,R,q): A verbose form of GuessLehmer1(X,R,q). Try:
#GuessLehmer1v(X,5,q);
GuessLehmer1v:=proc(X,R,q) local gu,r,n1,i,Le,n,Q,mu,i1,N,j,a,lu:
print(`Let Le(n) by the Lehmer n by n triadiagonal  matrix whose entries are zero `):
print(`except at the main diagonal and the diagonals right above and right below where for i from 1 to n [when applicble]`):
print(Le(n)[i,i]=1, Le(n)[i,i+1]=sqrt(X)*q^((i-1)/2), Le(n)[i,i-1]=sqrt(X)*q^((i-2)/2)):
print(``):
print(`For example the  6 by 6 case looks as follows.`):
print(``):
print(matrix(M(X,6,q))):
print(``):
print(`and Let Q(n)(X,q) be its determinant. We will use George Andrews' pioneering reverse engineering approach to`):
print(`discover an explicit expression for these polynomials. `):
print(`By either directly computing the determinant, or better still using the obvious recurrence pointed out by D.H. Lehmer`):
print(Q(n)(X,q)=Q(n-1)(X,q)-q^(n-2)*X*Q(n-2)(X,q)):
print(``):
gu:=[seq(QXn(X,n1,q),n1=1..4*R+10)]:
print(``):
print(`Just for fun, here are the first`,10, `polynomials. `):
print(``):
for i from 1 to 10 do
print(``):
print(Q(i)(X,q)=gu[i]):
print(``):
od:

mu:=[seq(GuessPol([seq(coeff(gu[i],X,r),i=1..nops(gu))],q,N,2*r),r=1..R)]:

print(`Let's try to tackle each coefficient one at a time, and then guess a polynomial in`, q^n, `that we will call`,  N):
print(`The coefficient of X^0, i.e. the constant term, is always`, 1):

for i from 1 to R do
print(`the sequence of coefficients of`, X^i,  `for n from 1 to`, 2*i+4, ` is`):
print(``):
print(seq(coeff(gu[i1],X,i),i1=1..2*i+4)):
print(``):
print(`You can see that it is 0 until n=`, 2*i-1, `but then it fits the polynomial`):
print(``):
print(mu[i]):
print(``):
od:

print(`Let's look at the sequence of polynomials in`, N=q^n, `that we got so far that are the coefficients of`, X^a, `of `, Q(n)(X,q) ):
print(`for a from 1 to`, R):
print(``):
print(mu):
print(``):
print(`You don't have to be a Ramanujan, or George Andrews to conjecture that the numerator of the coefficient of `, X^a, `is always (so far)`):
print(``):
print(Product(N-q^j,j=a..2*a-1)):
print(``):
print(`Also obvious is the factor`):
print(``):
print(q^((a+1)*a/2)*(1-q)^a):
print(``):
print(`in the denominator. `):
print(``):
lu:=[seq(normal(mu[i]/mul(N-q^j,j=i..2*i-1)*q^((i+1)*i/2 )*(1-q^i)),i=1..nops(mu))]:

print(`dividing by the former and multiplying by the latter, gives the following sequence for a from 1 to,`, R):
print(``):
print(lu):
print(``):
print(`We still have to guess an explicit expression for it. `):
print(``):
print(`Let's look at the sequence of consecutive ratios`):
print(``):
lu:=[seq(expand(normal(lu[i]/lu[i+1])),i=1..nops(lu)-1)]:
print(``):
print(lu):
print(``):
print(`These are obviously`, 1-q^a, `hence a conjectured expression for the coefficient of X^a in the Lehner polynomial is`):
print(``):
print(
Product(N-q^j,j=a..2*a-1)/(q^((a+1)*a/2) *(1-q)^a*Product(1-q^j,j=1..a-1))
) :
print(``):
print(`Plugging in N=q^n we conjecture that`):
print(``):

print(Q(n)(X,q)=Sum((-1)^a*X^a*q^(a*(a-1))*qBinomial(n-a,a),a=0..trunc(n/2) )):

mu:
end: