######################################################################
# CLD: Save this file as CLD. To use it, stay in the                 #
# same directory, get into Maple (by typing: maple <Enter> )         #
# and then type:  read CLD : <Enter>                                 #
# Then follow the instructions given there                           #
#                                                                    #
# Written by Doron Zeilberger, Rutgers University ,                  #
# zeilberg@math.rutgers.edu.                                         # 
######################################################################

#Created: Oct. 11, 2002
#This version: Oct. 11, 2002
#CLD: A Maple package to automatically conjecture and
#rigorously prove explicit evaluations for Hankel and
#Topelitz determinants that evaluate to hyperhypergeometric
#expressions, using Dodgson's rule
#It accompanies the article:
# Liebe Opa Paul: Ich Bin Auch Ein Experimental Scientist
#Please report bugs to zeilberg@math.rutgers.edu



print(`Created: Oct. 11, 2002.`):
print(`This version: Oct. 11, 2002`):
lprint(``):
print(`CLD: A Maple package to automatically conjecture and`):
print(`rigorously prove explicit evaluations for Hankel and`):
print(`Topelitz determinants that evaluate to super-nice expressions`):
print(`using Dodgson's rule`):
print(``):
print(`Written by Doron Zeilberger, zeilberg@math.rutgers.edu`):
print(``):
print(`It accompanies the article:`):
print(` "Liebe Opa Paul: Ich Bin Auch Ein Experimental Scientist"`):
print(` available from Zeilberger's website`):
lprint(``):
print(`Please report bugs to zeilberg@math.rutgers.edu`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/`):
 print(`For a list of the IMPORTANT procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);`):
 print(`For a list of ALL  procedures, type, ezra1(); `):
 print(``):
ezra1:=proc()
if args=NULL then
 print(`Contains the following procedures: CheckHankel, `):
 print(` FindfH, FindgH, Findf1H, Findf2H, Findg1H, findg2H, `):
 print(` EvalH, EvalHp, EvalHpaper`):
 print(` GenPol2, GenRat1, GenRat2, , Guessf1H, Guessf2H, Guessg2H`):
 print(` GuessRat1,GuessRat2, Hankel, IsYaya`):
 print(``):
 print(`For help with a specific procdure type: ezra(proc_name); `):
 else
   ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then
 print(`The main procedures are:`):
 print(` EvalH, EvalHpaper, EvalT, EvalTpaper`):
 print(`For a list of all procedures, type: ezra1(); `):
fi:

if nops([args])=1 and op(1,[args])=CheckHankel then
print(`CheckHankel(PP,n,r): checks whether the  yaya PP`):
print(`satisfies the Dodgson Recurrence for Hankel determinants`):
print(`For example, try :`):
print(` CheckHankel([[n*(n+r),n+r],[n+r,(n+r)/(r-1)]],n,r);`):
fi:

if nops([args])=1 and op(1,[args])=FindfH then
print(`FindfH(h,r,n0,r0): Hankel(h,r,n0,r0)/Hankel(h,r,n0-1,r0)`):
fi:

if nops([args])=1 and op(1,[args])=FindgH then
print(`FindgH(h,r,n0,r0): Hankel(h,r,n0,r0)/Hankel(h,r,n0,r0-1)`):
fi:

if nops([args])=1 and op(1,[args])=Findf1H then
print(`Findf1H(h,r,n0,r0): FindfH(h,r,n0,r0)/FindfH(h,r,n0-1,r0)`):
fi:

if nops([args])=1 and op(1,[args])=Findf2H then
print(`Findf2H(h,r,n0,r0): FindfH(h,r,n0,r0)/FindfH(h,r,n0,r0-1)`):
fi:
 
if nops([args])=1 and op(1,[args])=Findg1H then
print(`Findg1H(h,r,n0,r0): FindgH(h,r,n0,r0)/FindgH(h,r,n0-1,r0)`):
fi:

if nops([args])=1 and op(1,[args])=Findg2H then
print(`Findg2H(h,r,n0,r0): FindgH(h,r,n0,r0)/FindgH(h,r,n0,r0-1)`):
fi:

if nops([args])=1 and op(1,[args])=EvalH then
print(`EvalH(h,n,r,MaxDeg): Guesses a Yaya expression `):
print(`of degree <=MaxDeg, for`):
print(`the Hankel determinant det(h(r+i+j)) (0<=i,j<=n-1)`):
print(`and proves it at the same time. For example, try:`):
print(` EvalH(r!,n,r,4); `):
fi:

if nops([args])=1 and op(1,[args])=EvalHp then
print(`EvalHp(h,n,r,MaxDeg,p,Degp): Guesses a Yaya expression `):
print(`of degree <=MaxDeg, for`):
print(`the Hankel determinant det(h(r+i+j)) (0<=i,j<=n-1)`):
print(`where h also depends on a parameter p, and the ratios`):
print(`are expected to be of degree<=Degp in p`):
print(`and proves it at the same time. For example, try:`):
print(`EvalHp(1/(m+r)!/(m-r)!,n,r,4,m,2);`):
fi:

if nops([args])=1 and op(1,[args])=EvalHpaper then
print(`EvalHpaper(h,n,r,MaxDeg):Inputs an expression (hypegeometric)`):
print(` h in the variable r, tries to find  `):
print(` the Hankel determinant det(h(r+i+j)) (0<=i,j<=n-1) `):
print(` belonging to the yaya ansatz `):
print(` and proves it at the same time. If succesful, outputs a paper. `):
print(` For example try: EvalHpaper(r!,n,r,4);`):
fi:

if nops([args])=1 and op(1,[args])=EvalT then
print(`EvalT(h,n,r,MaxDeg): Guesses a Yaya expression `):
print(`of degree <=MaxDeg, for`):
print(`the Toeplitz determinant det(h(r+i+j)) (0<=i,j<=n-1)`):
print(`and proves it at the same time. For example, try:`):
print(` EvalT(1/(r+11)!,n,r,4); `):
fi:

if nops([args])=1 and op(1,[args])=EvalTpaper then
print(`EvalTpaper(h,n,r,MaxDeg):Inputs an expression (hypegeometric)`):
print(` h in the variable r, tries to find  `):
print(` the Toeplitz determinant det(h(r+i-j)) (0<=i,j<=n-1) `):
print(` belonging to the yaya ansatz `):
print(` and proves it at the same time. If succesful, outputs a paper. `):
print(` For example try: EvalTpaper(1/(r+11)!,n,r,4);`):
fi:


if nops([args])=1 and op(1,[args])=GenPol2 then
print(`GenPol2(n,r,deg,a): Given an integer deg, and symbols n and r`):
print(`and another symbol a, outputs a generic polynomial in n,r, of`):
print(`total degree deg, with coefficients indexed by a, followed`):
print(`by the set of coefficients. For example, try: GenPol2(n,r,1,a);`):
fi:

if nops([args])=1 and op(1,[args])=GenRat1 then
print(` GenRat1(n,deg,a,b): Given an integer deg, and symbol n`):
print(`and two other symbols a,b outputs a generic rational function in n of`):
print(`degree deg (of both numerator and denominator), `):
print(`with numerator (denom.) coefficients indexed by a (b), followed`):
print(`by the set of coefficients. For example, try: GenRat1(n,1,a,b); `):
fi:

if nops([args])=1 and op(1,[args])=GenRat2 then
print(` GenRat2(n,r,deg,a,b): Given an integer deg, and symbols n and r`):
print(`and another symbol a, outputs a generic rational function in n,r, of`):
print(`total degree deg (of both numerator and denominator), `):
print(`with numerator (denom.) coefficients indexed by a (b), followed`):
print(`by the set of coefficients. For example, try: GenRat2(n,r,1,a,b); `):
fi:

if nops([args])=1 and op(1,[args])=Guessf1H then
print(`Guessf1H(h,n,r,MaxDeg): Given an expression h in the variable r`):
print(`Guesses the (hopefully rational) function`):
print(` (in the variables r and n)  up to MaxDeg`):
print(`fitting Findf1H(h,r,n0,r0): For example, try:`):
print(`Guessf1H(r!,n,r,1);`):
fi:

if nops([args])=1 and op(1,[args])=Guessf2H then
print(`Guessf2H(h,n,r,MaxDeg): Given an expression h in the variable r`):
print(`Guesses the (hopefully rational) function`):
print(` (in the variables r and n)  up to MaxDeg`):
print(`fitting Findf2H(h,r,n0,r0): For example, try:`):
print(`Guessf2H(r!,n,r,1);`):
fi:

if nops([args])=1 and op(1,[args])=Guessg2H then
print(`Guessg2H(h,n,r,MaxDeg): Given an expression h in the variable r`):
print(`Guesses the (hopefully rational) function`):
print(` (in the variables r and n)  up to MaxDeg`):
print(`fitting Findg2H(h,r,n0,r0): For example, try:`):
print(`Guessg2H(r!,n,r,1);`):
fi:

if nops([args])=1 and op(1,[args])=GuessRat1 then
print(`GuessRat1(Data,n,deg): Given a list of`):
print(`data, Data,`):
print(` tries to finds a rational function R(n) (whose numerator`):
print(` and denominator) are of degree<=deg, in n  `):
print(` such that R(i)=Data[i] `):
print(` For example, try: GuessRat1([1,1,1,1,1],n,0);`):
fi:
 
if nops([args])=1 and op(1,[args])=GuessRat2 then
print(`GuessRat2(Data,n,r,deg): Given, Data,  a list of lists of`):
print(`data tries to finds a rational function R(n,r) (whose numerator`):
print(` and denominator) are of degree<=deg, in n and r `):
print(` such that R(i,j)=Data[i][j] `):
print(` For example, try: GuessRat2([[1,1,1],[1,1,1],[1,1,1]],n,r,1);`):
fi:

if nops([args])=1 and op(1,[args])=Hankel then
print(`Hankel(h,r,n0,r0): Given an expression h(r), and integers`):
print(`n0,r0>=0, finds det(h(r0+i+j),0<=i,j<=n0-1)`):
print(` by using Dogson's rule. For example, try `):
print(` Hankel(r!,r,3,2); `):
fi:

if nops([args])=1 and op(1,[args])=IsYaya then
print(`IsYaya(yaya,n,r): Check whether YaYa is a legal yaya in n and r`):
print(`For example, try IsYaya([[(n+r)/n,n+r] , [n+r, (n+r)/r]],n,r);`):
fi:

end:


#IsYaya(yaya,n,r): Check whether YaYa is a legal yaya in n and r
#For example, try IsYaya([[(n+r)/n,n+r] , [n+r, (n+r)/r]],n,r);
IsYaya:=proc(PP,n,r) local f1,f2,f3:
if not type(PP,list) then
print(PP, `should be a list`):
RETURN(false):
fi:

f1:=PP[1]:f2:=PP[2]:f3:=PP[3]:


if normal(f1*subs(n=n-1,f2)-subs(r=r-1,f1)*f2)<>0 then
RETURN(false):
fi:

if normal(f2*subs(n=n-1,f3)-subs(r=r-1,f2)*f3)<>0 then
RETURN(false):
fi:
true:
end:


########General procedures for guessing

#GenPol1(n,deg,a): Given an integer deg, and symbol n
#and another symbol a, outputs a generic polynomial in n, of
#degree deg, with coefficients indexed by a, followed
#by the set of coefficients. For example, try: GenPol1(n,1,a);
GenPol1:=proc(n,deg,a)
local i,gu,mu:
gu:=0:
mu:={}:

for i from 0 to deg do
 gu:=gu+a[i]*n^i:
 mu:=mu union {a[i]}:
od:
gu,mu:
end:

#GenPol2(n,r,deg,a): Given an integer deg, and symbols n and r
#and another symbol a, outputs a generic polynomial in n,r, of
#total degree deg, with coefficients indexed by a, followed
#by the set of coefficients. For example, try: GenPol2(n,r,1,a);
GenPol2:=proc(n,r,deg,a)
local i,j,gu,mu:
gu:=0:
mu:={}:

for i from 0 to deg do
for j from 0 to i do
 gu:=gu+a[j,i-j]*n^j*r^(i-j):
 mu:=mu union {a[j,i-j]}:
od:
od:
gu,mu:
end:


#GenRat1(n,deg,a,b): Given an integer deg, and symbol n 
#and another symbol a, outputs a generic rational function in n of
#degree deg (of both numerator and denominator), 
#with numerator (denom.) coefficients indexed by a (b), followed
#by the set of coefficients. For example, try: GenRat1(n,r,1,a,b); 
GenRat1:=proc(n,deg,a,b) local gu1,gu2:
gu1:=GenPol1(n,deg,a):
gu2:=GenPol1(n,deg,b):
gu1[1]/gu2[1],gu1[2] union gu2[2]:
end:

#GenRat2(n,r,deg,a,b): Given an integer deg, and symbols n and r
#and another symbol a, outputs a generic rational function in n,r, of
#total degree deg (of both numerator and denominator), 
#with numerator (denom.) coefficients indexed by a (b), followed
#by the set of coefficients. For example, try: GenRat2(n,r,1,a,b); 
GenRat2:=proc(n,r,deg,a,b) local gu1,gu2:
gu1:=GenPol2(n,r,deg,a):
gu2:=GenPol2(n,r,deg,b):
gu1[1]/gu2[1],gu1[2] union gu2[2]:
end:

#GuessRat2(Data,n,r,deg): Given a list of lists of
#data tries to finds a rational function R(n,r) (whose numerator
#and denominator) are of degree<=deg, in n and r 
#such that R(i,j)=Data[i][j]
GuessRat2:=proc(Data,n,r,deg)
local che,i,j,eq,var,gu,a,b,R:
che:=convert([seq(nops(Data[i]), i=1..nops(Data))],`+`):
if not (che>(deg+1)*deg+3) then
 ERROR(`Insufficient data`):
fi:

gu:=GenRat2(n,r,deg,a,b):
var:=gu[2]:
R:=gu[1]:
eq:={}:

for i from 1 to nops(Data) do
for j from 1 to nops(Data[i]) do
eq:=eq union {numer(normal(subs({n=i,r=j},R)-Data[i][j]))}:
od:
od:
var:=solve(eq,var):
if subs(var,numer(R))=0 then
RETURN(0):
fi:
factor(normal(subs(var,R))):
end:


#GuessRat1(Data,n,deg): Given a list of 
#data tries to finds a rational function R(i) (whose numerator
#and denominator) are of degree<=deg, in n 
#such that R(i)=Data[i]
GuessRat1:=proc(Data,n,deg)
local che,i,eq,var,gu,a,b,R:
che:=nops(Data):
if not (che>deg+2) then
 ERROR(`Insufficient data`):
fi:

gu:=GenRat1(n,deg,a,b):
var:=gu[2]:
R:=gu[1]:
eq:={}:

for i from 1 to nops(Data) do
eq:=eq union {numer(normal(subs({n=i},R)-Data[i]))}:
od:
var:=solve(eq,var):
if subs(var,numer(R))=0 then
RETURN(0):
fi:
factor(normal(subs(var,R))):
end:


########End general procedures for guessing

##########Hankel-Specific##############
#CheckHankel(PP,n,r): checks whether the  yaya PP
#satisfies the Dodgson Recurrence for Hankel determinants
#For example, try CheckHankel([[n*(n+r),n+r],[n+r,(n+r)/(r-1)]],n,r);
CheckHankel:=proc(PP,n,r)
local f1,f2,f3,gu:
if not IsYaya(PP,n,r) then
print(PP, `is not a yaya`):
RETURN(false):
fi:
f1:=PP[1]:f2:=PP[2]:f3:=PP[3]:


gu:=normal(subs({n=n-1,r=r+2},f2)*subs({n=n-1,r=r+1},f2)/f1-
subs(r=r+2,f2)*subs({n=n-1,r=r+2},f2)/
subs(r=r+2,f3)/subs(r=r+1,f1)):
if gu=1 then RETURN(true) else RETURN(false) fi:
end:


#Hankel(h,r,n0,r0): Given an expression h(r), and integers
#n0,r0>=0, finds det(h(r0+i+j),0<=i,j<=n0-1)
#by using Dogson's rule. For example, try
#Hankel(r!,r,3,2);
Hankel:=proc(h,r,n0,r0)
option remember:

if n0=0 then
 RETURN(1):
fi:

if n0=1 then
 RETURN(subs(r=r0,h)):
fi:

(Hankel(h,r,n0-1,r0)*Hankel(h,r,n0-1,r0+2)-Hankel(h,r,n0-1,r0+1)^2)/
Hankel(h,r,n0-2,r0+2):
end:


FindfH:=proc(h,r,n0,r0):Hankel(h,r,n0,r0)/Hankel(h,r,n0-1,r0):end:
FindgH:=proc(h,r,n0,r0):Hankel(h,r,n0,r0)/Hankel(h,r,n0,r0-1):end:
Findf1H:=proc(h,r,n0,r0):simplify(FindfH(h,r,n0,r0)/FindfH(h,r,n0-1,r0)):end:
Findf2H:=proc(h,r,n0,r0):simplify(FindfH(h,r,n0,r0)/FindfH(h,r,n0,r0-1)):end:
Findg1H:=proc(h,r,n0,r0):simplify(FindgH(h,r,n0,r0)/FindgH(h,r,n0-1,r0)):end:
Findg2H:=proc(h,r,n0,r0):simplify(FindgH(h,r,n0,r0)/FindgH(h,r,n0,r0-1)):end:


#Guessf1H1(h,n,r,Deg): Guesses the (hopefully rational) function
#fittin Findf1H(h,r,n0,r0): For example, try:
#Guessf1H1(r!,r,1);
Guessf1H1:=proc(h,n,r,Deg)
local Data,n0,r0,lu:
Data:=[]:
for n0 from 2 to Deg+3 do
lu:=[]:
for r0 from 2 to Deg+3 do
lu:=[op(lu),Findf1H(h,r,n0,r0)]:
od:
Data:=[op(Data),normal(simplify(lu))]:
od:
subs({r=r-1,n=n-1},GuessRat2(Data,n,r,Deg)):
end:


#Guessf1H(h,n,r,MaxDeg): Guesses the (hopefully rational) function
#fittin Findf1H(h,r,n0,r0): up to maximal degree MaxDeg For example, try:
#Guessf1H(r!,r,1);

Guessf1H:=proc(h,n,r,MaxDeg)
local Deg,gu:
for Deg from  0 to MaxDeg do
gu:=Guessf1H1(h,n,r,Deg):
if gu<>0 then
RETURN(gu):
fi:
od:
print(`There is nothing of degree <= `, MaxDeg):
0:
end:


#Guessf2H1(h,n,r,Deg): Guesses the (hopefully rational) function
#fittin Findf2H(h,r,n0,r0): For example, try:
#Guessf2H1(r!,r,1);
Guessf2H1:=proc(h,n,r,Deg)
local Data,n0,r0,lu:
Data:=[]:
for n0 from 2 to Deg+3 do
lu:=[]:
for r0 from 2 to Deg+3 do
lu:=[op(lu),Findf2H(h,r,n0,r0)]:
od:
Data:=[op(Data),lu]:
od:
subs({r=r-1,n=n-1},GuessRat2(Data,n,r,Deg)):
end:


#Guessf2H(h,n,r,MaxDeg): Guesses the (hopefully rational) function
#fittin Findf2H(h,r,n0,r0): up to maximal degree MaxDeg For example, try:
#Guessf2H(r!,r,1);

Guessf2H:=proc(h,n,r,MaxDeg)
local Deg,gu:
for Deg from  0 to MaxDeg do
gu:=Guessf2H1(h,n,r,Deg):
if gu<>0 then
RETURN(gu):
fi:
od:
print(`There is nothing of degree <= `, MaxDeg):
0:
end:

 
#Guessg2H1(h,n,r,Deg): Guesses the (hopefully rational) function
#fittin Findg2H(h,r,n0,r0): For example, try:
#Guessg2H1(r!,r,1);
Guessg2H1:=proc(h,n,r,Deg)
local Data,n0,r0,lu:
Data:=[]:
for n0 from 2 to Deg+3 do
lu:=[]:
for r0 from 2 to Deg+3 do
lu:=[op(lu),Findg2H(h,r,n0,r0)]:
od:
Data:=[op(Data),lu]:
od:
subs({r=r-1,n=n-1},GuessRat2(Data,n,r,Deg)):
end:


#Guessg2H(h,n,r,MaxDeg): Guesses the (hopefully rational) function
#fittin Findg2H(h,r,n0,r0): up to maximal degree MaxDeg For example, try:
#Guessg2H(r!,r,1);

Guessg2H:=proc(h,n,r,MaxDeg)
local Deg,gu:
for Deg from  0 to MaxDeg do
gu:=Guessg2H1(h,n,r,Deg):
if gu<>0 then
RETURN(gu):
fi:
od:
print(`There is nothing of degree <= `, MaxDeg):
0:
end:


#EvalH(h,n,r,MaxDeg): Guesses a Yaya expression 
#of degree <=MaxDeg, for
#the Hankel determinant det(h(r+i+j)) (0<=i,j<=n-1)
#and proves it at the same time. For example, try:
#EvalH(r!,n,r,4);
EvalH:=proc(h,n,r,MaxDeg)
local gu1,gu2,gu3,Y,mu:

gu1:=Guessf1H(h,n,r,MaxDeg):
if gu1=0 then
 RETURN(0):
fi:
gu2:=Guessf2H(h,n,r,MaxDeg):
if gu2=0 then
 RETURN(0):
fi:

gu3:=Guessg2H(h,n,r,MaxDeg):
if gu3=0 then
 RETURN(0):
fi:

Y:=[gu1,gu2,gu3]:
if gu1=0 or gu2=0 or gu3=0 then
print(`Nothing found`):
RETURN(Y):
fi:

mu:=normal(simplify(h/subs(r=r-1,h))):
mu:=normal(simplify(mu/subs(r=r-1,mu))):
mu:=normal(subs(n=1,gu3)-mu):
if mu<>0 then
print(`Initial conditions do not match`):
RETURN(false):
fi:

if not CheckHankel(Y,n,r) then
print(`The guess`, Y, `is wrong `):
RETURN(0):
else
RETURN(Y):
fi:
end:





#EvalHpaper(h,n,r,MaxDeg):Inputs an expression (hypegeometric)
#h in the variable r, tries to find 
#the Hankel determinant det(h(r+i+j)) (0<=i,j<=n-1)
#belonging to the yaya ansatz
#and proves it at the same time. If succesful, outputs a paper
#For example try: EvalHpaper(r!,n,r,4);
EvalHpaper:=proc(h,n,r,MaxDeg)
local gu,f1:
gu:=EvalH(h,n,r,MaxDeg):
f1:=gu[1]:
if gu=0 then
print(`Sorry, no paper`):
RETURN(0):
fi:

print(`An Explicit Evaluation Of a Certain Hankel Determinant` ):
print(`By S. B. Ekhad (c/o D. Zeilberger, zeilberg@math.rutgers.edu)`):
print(``):
print(`Theorem: Let h(r):=`):
print(h):
print(`Let H(n,r) be the Hankel determinant`):
print(`det (h(r+i+j)), 0<=i,j<=n-1 `):
print(`Then we have the explicit evaluation`):
print(`H(n,r)= `):
print(h^n*Product(Product(subs(n=n[2],f1),n[2]=2..n[1]),n[1]=1..n)):
print(``):
print(`Proof: Let's call the expression on the r.h.s. G(n,r).`):
print(`Obviously  H(n,r)=G(n,r), for n=0,1`):
print(`It is readily checked that`):
print( G(n,r)=(G(n-1,r)*G(n-1,r+2)-G(n-1,r+2)^2)/G(n-2,r+2) ):
print(` Since, by Dodgson's rule: `):
print( H(n,r)=(H(n-1,r)*H(n-1,r+2)-H(n-1,r+2)^2)/H(n-2,r+2) ):
print(` It follows by induction on n that, for all n>=0 `):
print(H(n,r)=G(n,r)):
print(` QED `):
end:

#EvalHp(h,n,r,MaxDeg,p,Degp): Guesses a Yaya expression 
#of degree <=MaxDeg, for
#the Hankel determinant det(h(r+i+j)) (0<=i,j<=n-1)
#where h also depends on a parameter p
#and proves it at the same time. For example, try:
#EvalHp(1/(m+r)!/(m-r)!,n,r,4,m);
EvalHp:=proc(h,n,r,MaxDeg,p,Degp)
local K,lu,h1,i1,lu11,lu12,lu22,f11,f12,f22:
K:=MaxDeg+17:
h1:=subs(p=K,h):
lu:=EvalH(h1,n,r,MaxDeg):
if lu=0 then
RETURN(0):
fi:

lu:=[seq(EvalH(subs(p=i1,h),n,r,MaxDeg),i1=K+1..K+Degp+2)]:
lu11:=[seq(lu[i1][1],i1=1..nops(lu))]:
lu12:=[seq(lu[i1][2],i1=1..nops(lu))]:
lu22:=[seq(lu[i1][3],i1=1..nops(lu))]:

f11:=GuessRat1a(lu11,p,MaxDeg):
if f11=0 then
RETURN(0):
fi:

f12:=GuessRat1a(lu12,p,MaxDeg):
if f12=0 then
RETURN(0):
fi:

f22:=GuessRat1a(lu22,p,MaxDeg):
if f22=0 then
RETURN(0):
fi:

subs(p=p-K,[f11,f12,f22]):


end:

#########End of Hankel-specific

##########Toeplitz-Specific##############
#CheckToeplitz(PP,n,r): checks whether the  yaya PP
#satisfies the Dodgson Recurrence for Toeplitz determinants
#For example, try CheckToeplitz([[n*(n+r),n+r],[n+r,(n+r)/(r-1)]],n,r);
CheckToeplitz:=proc(PP,n,r)
local f1,f2,g2,gu:
if not IsYaya(PP,n,r) then
print(PP, `is not a yaya`):
RETURN(false):
fi:
f1:=PP[1]:f2:=PP[1]:g2:=PP[3]:


gu:=normal(1/f1-
subs({n=n-1,r=r+1},f2)*subs({n=n-2,r=r+1},g2)/
subs({n=n-1},f2)/f1):

if gu=1 then RETURN(true) else RETURN(false) fi:
end:


#Toeplitz(h,r,n0,r0): Given an expression h(r), and integers
#n0,r0>=0, finds det(h(r0+i-j),0<=i,j<=n0-1)
#by using Dogson's rule. For example, try
#Toeplitz((r+6)!,r,3,2);
Toeplitz:=proc(h,r,n0,r0)
option remember:

if n0=0 then
 RETURN(1):
fi:

if n0=1 then
 RETURN(subs(r=r0,h)):
fi:

(Toeplitz(h,r,n0-1,r0)^2-Toeplitz(h,r,n0-1,r0+1)*Toeplitz(h,r,n0-1,r0-1))/
Toeplitz(h,r,n0-2,r0):
end:

FindfT:=proc(h,r,n0,r0):simplify(Toeplitz(h,r,n0,r0)/Toeplitz(h,r,n0-1,r0)):
end:
FindgT:=proc(h,r,n0,r0):simplify(Toeplitz(h,r,n0,r0)/Toeplitz(h,r,n0,r0-1)):
end:
Findf1T:=proc(h,r,n0,r0):simplify(FindfT(h,r,n0,r0)/FindfT(h,r,n0-1,r0)):end:
Findf2T:=proc(h,r,n0,r0):simplify(FindfT(h,r,n0,r0)/FindfT(h,r,n0,r0-1)):end:
Findg1T:=proc(h,r,n0,r0):simplify(FindgT(h,r,n0,r0)/FindgT(h,r,n0-1,r0)):end:
Findg2T:=proc(h,r,n0,r0):simplify(FindgT(h,r,n0,r0)/FindgT(h,r,n0,r0-1)):end:

#Guessf1T1(h,n,r,Deg): Guesses the (hopefully rational) function
#fittin Findf1T(h,r,n0,r0): For example, try:
#Guessf1T1(r!,r,1);
Guessf1T1:=proc(h,n,r,Deg)
local Data,n0,r0,lu:
Data:=[]:
for n0 from 2 to Deg+3 do
lu:=[]:
for r0 from 2 to Deg+3 do
lu:=[op(lu),Findf1T(h,r,n0,r0)]:
od:
Data:=[op(Data),normal(simplify(lu))]:
od:
subs({r=r-1,n=n-1},GuessRat2(Data,n,r,Deg)):
end:


#Guessf1T(h,n,r,MaxDeg): Guesses the (hopefully rational) function
#fittin Findf1H(h,r,n0,r0): up to maximal degree MaxDeg For example, try:
#Guessf1T((r+10)!,r,1);

Guessf1T:=proc(h,n,r,MaxDeg)
local Deg,gu:
for Deg from  0 to MaxDeg do
gu:=Guessf1T1(h,n,r,Deg):
if gu<>0 then
RETURN(gu):
fi:
od:
print(`There is nothing of degree <= `, MaxDeg):
0:
end:


#Guessf2T1(h,n,r,Deg): Guesses the (hopefully rational) function
#fittin Findf2H(h,r,n0,r0): For example, try:
#Guessf2T1((r+5)!,r,1);
Guessf2T1:=proc(h,n,r,Deg)
local Data,n0,r0,lu:
Data:=[]:
for n0 from 2 to Deg+3 do
lu:=[]:
for r0 from 2 to Deg+3 do
lu:=[op(lu),Findf2T(h,r,n0,r0)]:
od:
Data:=[op(Data),lu]:
od:
subs({r=r-1,n=n-1},GuessRat2(Data,n,r,Deg)):
end:


#Guessf2T(h,n,r,MaxDeg): Guesses the (hopefully rational) function
#fittin Findf2T(h,r,n0,r0): up to maximal degree MaxDeg For example, try:
#Guessf2T((r+10)!,r,1);

Guessf2T:=proc(h,n,r,MaxDeg)
local Deg,gu:
for Deg from  0 to MaxDeg do
gu:=Guessf2T1(h,n,r,Deg):
if gu<>0 then
RETURN(gu):
fi:
od:
print(`There is nothing of degree <= `, MaxDeg):
0:
end:

 
#Guessg2T1(h,n,r,Deg): Guesses the (hopefully rational) function
#fittin Findg2T(h,r,n0,r0): For example, try:
#Guessg2T1(r!,r,1);
Guessg2T1:=proc(h,n,r,Deg)
local Data,n0,r0,lu:
Data:=[]:
for n0 from 2 to Deg+3 do
lu:=[]:
for r0 from 2 to Deg+3 do
lu:=[op(lu),Findg2T(h,r,n0,r0)]:
od:
Data:=[op(Data),lu]:
od:
subs({r=r-1,n=n-1},GuessRat2(Data,n,r,Deg)):
end:


#Guessg2T(h,n,r,MaxDeg): Guesses the (hopefully rational) function
#fittin Findg2T(h,r,n0,r0): up to maximal degree MaxDeg For example, try:
#Guessg2T(r!,r,1);

Guessg2T:=proc(h,n,r,MaxDeg)
local Deg,gu:
for Deg from  0 to MaxDeg do
gu:=Guessg2T1(h,n,r,Deg):
if gu<>0 then
RETURN(gu):
fi:
od:
print(`There is nothing of degree <= `, MaxDeg):
0:
end:


#EvalT(h,n,r,MaxDeg): Guesses a Yaya expression 
#of degree <=MaxDeg, for
#the Toeplitz determinant det(h(r+i+j)) (0<=i,j<=n-1)
#and proves it at the same time. For example, try:
#EvalT(r!,n,r,4);
EvalT:=proc(h,n,r,MaxDeg)
local gu1,gu2,gu3,Y,mu:

gu1:=Guessf1T(h,n,r,MaxDeg):
if gu1=0 then
 RETURN(0):
fi:
gu2:=Guessf2T(h,n,r,MaxDeg):
if gu2=0 then
 RETURN(0):
fi:

gu3:=Guessg2T(h,n,r,MaxDeg):
if gu3=0 then
 RETURN(0):
fi:

Y:=[gu1,gu2,gu3]:
if gu1=0 or gu2=0 or gu3=0 then
print(`Nothing found`):
RETURN(Y):
fi:

mu:=normal(simplify(h/subs(r=r-1,h))):
mu:=normal(simplify(mu/subs(r=r-1,mu))):
mu:=normal(subs(n=1,gu3)-mu):
if mu<>0 then
print(`Initial conditions do not match`):
RETURN(false):
fi:

if not CheckToeplitz(Y,n,r) then
print(`The guess`, Y, `is wrong `):
RETURN(0):
else
RETURN(Y):
fi:
end:





#EvalTpaper(h,n,r,MaxDeg):Inputs an expression (hypegeometric)
#h in the variable r, tries to find 
#the Toeplitz determinant det(h(r+i+j)) (0<=i,j<=n-1)
#belonging to the yaya ansatz
#and proves it at the same time. If succesful, outputs a paper
#For example try: EvalTpaper((r+10)!,n,r,4);
EvalTpaper:=proc(h,n,r,MaxDeg)
local gu,f1:
gu:=EvalT(h,n,r,MaxDeg):
f1:=gu[1]:
if gu=0 then
print(`Sorry, no paper`):
RETURN(0):
fi:

print(`An Explicit Evaluation Of a Certain Toeplitz Determinant` ):
print(`By S. B. Ekhad (c/o D. Zeilberger, zeilberg@math.rutgers.edu)`):
print(``):
print(`Theorem: Let h(r):=`):
print(h):
print(`Let H(n,r) be the Toeplitz determinant`):
print(`det (h(r+i-j)), 0<=i,j<=n-1 `):
print(`Then we have the explicit evaluation`):
print(`H(n,r)= `):
print(h^n*Product(Product(subs(n=n[2],f1),n[2]=2..n[1]),n[1]=1..n)):
print(``):
print(`Proof: Let's call the expression on the r.h.s. G(n,r).`):
print(`Obviously  H(n,r)=G(n,r), for n=0,1`):
print(`It is readily checked that`):
print( G(n,r)=(G(n-1,r)^2-G(n-1,r+1)*G(n-1,r-1))/G(n-2,r) ):
print(` Since, by Dodgson's rule: `):
print( H(n,r)=(H(n-1,r)^2-H(n-1,r+1)*H(n-1,r-1))/H(n-2,r) ):
print(` It follows by induction on n that, for all n>=0 `):
print(H(n,r)=G(n,r)):
print(` QED `):
end:

#########End of Toeplitz-specific############



#GuessRat1a(Data,n,deg): Given a list of 
#data tries to finds a rational function R(i) (whose numerator
#and denominator) are of degree<=deg, in n 
#such that R(i)=Data[i]
GuessRat1a:=proc(Data,n,deg)
local lu,d:

for d from 0 to deg do
lu:=GuessRat1(Data,n,d):
if lu<>0 then
RETURN(lu):
fi:
od:
0:
end: