#VERSION OF Aug. 2, 2002
#Previous VERSION: 20, 2000 
#DODGSON: A Maple PACKAGE FOR Conjecturing and proving
#Determinant Identities by Dodgson's Condensation Method
#Version of June 4, 1997
#Written by Tewodros Amdeberhan (deVry Inst.) and
#Doron Zeilberger, Rutgers University
 
#Accompanies the paper "Through the AUTOMATED Looking Glass" 
#published in Adv. Appl. Math. v. 27 (2001), pp. 225-230 
 
print(` DODGSON: A Maple package for conjecturing and proving `):
print(` Determinant Identities by Dodgson's Condensation Method `):
print(`Accompanies the paper "Through the AUTOMATED Looking Glass" `  ):
print(` published in Adv. Appl. Math. v. 27 (2001), pp. 225-230 `):
 print(` Revised Version of : Aug. 2 2002 `):
 print(`Former Version: Oct. 20 2000`):
 print(`Tested on Maple7 and below`):
 print(`written by Tewodros Amdeberhan`):
 print(`and Doron Zeilberger(zeilberg@math.rutgers.edu).`):
 lprint(``):
 print(`The most current version is always available from`):
 print(`http://www.math.rutgers.edu/~zeilberg`):
 print(`For a list of the procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
ezra:=proc()
if args=NULL then

 print(`Contains the following procedures:`):
 print(` Hankelf, Hankelg, HankelProofTerse, HankelProof `):
 print(` Toeplitzf, Toeplitzg , ToeplitzProof`):
fi:
 
 
if nops([args])=1 and op(1,[args])=`Hankelf` then
print(` Hankelf(h,n,r,N,R) or Hankelf(h,n,r,N,R,MAX) `):
print(` Given a procedure h, (whose argument is of type integer) `):
print(` and variable names `):
print(` n,r, and a name for the shift operator in n,  N `):
print(`and a name for the shift operator in r, R `):
print(` and (optionaly) an integer MAX  that bounds the order+degrees `):
print(`of the recurrences looked for, (default=10). `):
print(``):
print(` Define A(n,r):=det(h(r+i+j),0<=i,j<=n-1)) and `):
print(` f(n,r):=A(n+1,r)/A(n,r)  looks empirically for`):
print(` linear recurrence operators,ope1(n,r,N) and ope2(n,r,R) `):
print(` that annihilate f(n,r) ` ):
print(` The output is : ope1,ope2 `):
print(``):
print(`Warning: During any Maple session you can't redefine`):
print(`the procedure h `);
print(` Once again: Hankelf(h,n,r,N,R) or Hankelf(h,n,r,N,R,MAX) `):
print(` For example: Hankelf(n->n!,n,r,N,R); `):
fi:
 
 
 
if nops([args])=1 and op(1,[args])=`Hankelg` then
print(` Hankelg(h,n,r,N,R) or Hankelg(h,n,r,N,R,MAX) `):
print(` Given a procedure h, (whose argument is of type integer) `):
print(` and variable names `):
print(` n,r, and a name for the shift operator in n,  N `):
print(`and shift operator in r, R `):
print(` and (optionaly) an integer MAX  that bounds the order+degrees `):
print(`of the recurrences looked for, (default=10). `):
print(``):
print(` Define A(n,r):=det(h(r+i+j),0<=i,j<=n-1)) and `):
print(` g(n,r):=A(n,r+1)/A(n,r)  looks empirically for`):
print(` linear recurrence operators,ope1(n,r,N) and ope2(n,r,R) `):
print(` that annihilate g(n,r) ` ):
print(` The output is : ope1,ope2 `):
print(``):
print(`Warning: During any Maple session you can't redefine`):
print(`the procedure h `);
print(` Once again: Hankelg(h,n,r,N,R) or Hankelg(h,n,r,N,R,MAX) `):
fi:
 
if nops([args])=1 and op(1,[args])=`HankelProofTerse` then
print(`HankelProofTerse(h,n,r)  `):
print(`Given a procedure h (from integers to integers) tries`):
print(`to guess empirically, and then proves it rigorously,`):
print(`an explicit formula for the Hankel determinant`):
print(`A(n,r):=det(h(r+i+j)) (0<=i,j<=n-1)`):
print(` For example, try:HankelProofTerse(n->n!,n,r)`):
fi:

if nops([args])=1 and op(1,[args])=`HankelProof` then
print(`HankelProof(h,n,r): `):
print(` Verbose version of HankelProofTerse(h,n,r) (q.v.) e.g. type: `):
print(` HankelProof(n->n!,n,r): `):
fi:
 
if nops([args])=1 and op(1,[args])=`Toeplitzf` then
print(` Toeplitzf(h,n,r,N,R) or Toeplitzf(h,n,r,N,R,MAX) `):
print(` Given a procedure h, (whose argument is of type integer) `):
print(` and variable names `):
print(` n,r, and a name for the shift operator in n,  N `):
print(` and (optionaly) an integer MAX  that bounds the order+degrees `):
print(`of the recurrences looked for, (default=10). `):
print(``):
print(` Define A(n,r):=det(h(r+i-j),0<=i,j<=n-1)) and `):
print(` f(n,r):=A(n+1,r)/A(n,r)  looks empirically for`):
print(` linear recurrence operators,ope1(n,r,N) and ope2(n,r,R) `):
print(` that annihilate f(n,r) ` ):
print(` The output is : ope1,ope2 `):
print(`Warning: During any Maple session you can't redefine`):
print(`the procedure h `);
print(``):
print(` Once again: Toeplitzf(h,n,r,N,R) or Toeplitzf(h,n,r,N,R,MAX) `):
fi:
 
 
 
if nops([args])=1 and op(1,[args])=`Toeplitzg` then
print(` Toeplitzg(h,n,r,N,R) or Toeplitzg(h,n,r,N,R,MAX) `):
print(` Given an expression in t, h(t), variable names `):
print(` n,r, and a names for the shift operator in n,  N `):
print(`and shift operator in r, R `):
print(` and (optionaly) an integer MAX  that bounds the order+degrees `):
print(`of the recurrences looked for, (default=10). `):
print(``):
print(` Define A(n,r):=det(h(r+i-j),0<=i,j<=n-1)) and `):
print(` g(n,r):=A(n,r+1)/A(n,r)  looks empirically for`):
print(` linear recurrence operators,ope1(n,r,N) and ope2(n,r,R) `):
print(` that annihilate g(n,r) ` ):
print(` The output is : ope1,ope2 `):
print(`Warning: During any Maple session you can't redefine`):
print(`the procedure h `);
print(``):
print(` Once again: Toeplitzg(h,n,r,N,R) or Toeplitzg(h,n,r,N,R,MAX) `):
fi:
 
if nops([args])=1 and op(1,[args])=`ToeplitzProof` then
print(`ToeplitzProof(h,n,r)  `):
print(`Given a procedure h (from integers to integers) tries`):
print(`to guess empirically, and then proves it rigorously,`):
print(`an explicit formula for the Toeplitz determinant`):
print(`A(n,r):=det(h(r+i-j)) (0<=i,j<=n-1)`):
print(` For example, try:ToeplitzProof(n->1/(a+n)!,n,r)`):
fi:
 
end:
 
#########Start Hankel
AHankel:=proc(h,i,j)
option remember:
if i=0 then
 RETURN(1):
elif i=1 then
 RETURN(h(j)):
else
RETURN((AHankel(h,i-1,j+2)*AHankel(h,i-1,j)-
AHankel(h,i-1,j+1)^2)/AHankel(h,i-2,j+2)):
fi:
end:
 
 
fHankel:=proc(h,i,j):
normal(expand(AHankel(h,i+1,j)/AHankel(h,i,j))):
end:
 
 
gHankel:=proc(h,i,j):
normal(expand(AHankel(h,i,j+1)/AHankel(h,i,j))):
end:
 
 
HankelfA11:=proc(h,ORD,DEG1,DEG2,n1,n2,N1)
local ope1,i,j,c,eq,var,n1a,n2a,ope2,lu,c1,eq1,i1,i2,Godel:
 
 
if not type(ORD,integer) or ORD<0 then
ERROR(`Second argument, is the ORDER, must be a non-negative integer`):
fi:
 
if not type(DEG1,integer) or DEG1<0 then
print(`Third argument is the maximal degree of the coefficients`):
print(` in`,n1):
ERROR(`must be a non-negative integer`):
fi:
 
 
if not type(DEG2,integer) or DEG2<0 then
print(`Fourth argument is the maximal degree of the coefficients`):
print(` in`,n2):
ERROR(`must be a non-negative integer`):
fi:
 
 
ope1:=0:
 
var:={}:
 
for i1 from 0 to DEG1 do
for i2 from 0 to DEG2 do
for j from 0 to ORD do
ope1:=ope1+c[i1,i2,j]*n1^i1*n2^i2*N1^j:
var:=var union {c[i1,i2,j]}:
od:
od:
od:
 
eq:={}:
 
 
Godel:=trunc(evalf(sqrt((DEG1+1)*(DEG2+1)*(ORD+1))))+3:
for n1a from 0 to  Godel do
for n2a from 0 to  Godel do
lu:=0:
 
for i from 0 to ORD do
lu:=lu+subs({n1=n1a,n2=n2a},coeff(ope1,N1,i))*fHankel(h,n1a+i,n2a):
od:
eq:=eq union {lu}:
 
od:od:
 
var:=solve(eq,var):
 
if var=NULL then
RETURN(0):
fi:
 
 
c1:={}:
 
for i from 1 to nops(var) do
eq1:=op(i,var):
 
 if op(1,eq1)=op(2,eq1) then
    c1:=c1 union {op(1,eq1)}:
 fi:
od:
 
if nops(c1)>1 then
#print(`Order or Degrees too high, lower them`):
  RETURN(0):
fi:
 
if nops(c1)=0 then
RETURN(0):
fi:
 
c1:=op(1,c1):
 
ope1:=subs(var,ope1):
 
ope2:=0:
 
ope1:=normal(ope1/c1):
 
for i from 0 to ORD do
ope2:=ope2+factor(coeff(ope1,N1,i))*N1^i:
od:
 
ope2:
end:
 
 
HankelfA1:=proc(h,n,r,N,MAX)
local mu,ORD,DEG1,DEG2,gu:
 
 
 
for mu from 0 to MAX do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
 
gu:=HankelfA11(h,1,DEG1,DEG2,n,r,N):
 
if gu<>0 then
RETURN(gu):
fi:
od:
od:
 
 
for ORD from 2 to MAX do
for mu from 0 to MAX-ORD do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
gu:=HankelfA11(h,ORD,DEG1,DEG2,n,r,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
od:
od:
 
print(`No recurrence found whose ORDER+Sum of Degrees <=`,MAX):
 
0:
end:
 
 
 
 
HankelfA:=proc():
 
if nargs=5 then
HankelfA1(args):
elif nargs=4 then
HankelfA1(args,10):
else
ERROR(`HankelfA1(h,n,r,N,MAX) or HankelfA1(h,n,r,N)`):
fi:
end:
 
 
 
 
 
HankelfB11:=proc(h,ORD,DEG1,DEG2,n1,n2,N2)
local ope1,i,j,c,eq,var,n1a,n2a,ope2,lu,c1,eq1,i1,i2,Godel:
 
 
if not type(ORD,integer) or ORD<0 then
ERROR(`Second argument, is the ORDER, must be a non-negative integer`):
fi:
 
if not type(DEG1,integer) or DEG1<0 then
print(`Third argument is the maximal degree of the coefficients`):
print(` in`,n1):
ERROR(`must be a non-negative integer`):
fi:
 
 
if not type(DEG2,integer) or DEG2<0 then
print(`Fourth argument is the maximal degree of the coefficients`):
print(` in`,n2):
ERROR(`must be a non-negative integer`):
fi:
 
 
ope1:=0:
 
var:={}:
 
for i1 from 0 to DEG1 do
for i2 from 0 to DEG2 do
for j from 0 to ORD do
ope1:=ope1+c[i1,i2,j]*n1^i1*n2^i2*N2^j:
var:=var union {c[i1,i2,j]}:
od:
od:
od:
 
eq:={}:
 
 
Godel:=trunc(evalf(sqrt((DEG1+1)*(DEG2+1)*(ORD+1))))+3:
for n1a from 0 to  Godel do
for n2a from 0 to  Godel do
lu:=0:
 
for i from 0 to ORD do
lu:=lu+subs({n1=n1a,n2=n2a},coeff(ope1,N2,i))*fHankel(h,n1a,n2a+i):
od:
eq:=eq union {lu}:
 
od:od:
 
var:=solve(eq,var):
 
if var=NULL then
RETURN(0):
fi:
 
 
c1:={}:
 
for i from 1 to nops(var) do
eq1:=op(i,var):
 
 if op(1,eq1)=op(2,eq1) then
    c1:=c1 union {op(1,eq1)}:
 fi:
od:
 
if nops(c1)>1 then
#print(`Order or Degrees too high, lower them`):
  RETURN(0):
fi:
 
if nops(c1)=0 then
RETURN(0):
fi:
 
c1:=op(1,c1):
 
ope1:=subs(var,ope1):
 
ope2:=0:
 
ope1:=normal(ope1/c1):
 
for i from 0 to ORD do
ope2:=ope2+factor(coeff(ope1,N2,i))*N2^i:
od:
 
ope2:
end:
 
 
 
 
 
HankelfB1:=proc(h,n,r,R,MAX)
local mu,ORD,DEG1,DEG2,gu:
 
 
 
for mu from 0 to MAX do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
 
gu:=HankelfB11(h,1,DEG1,DEG2,n,r,R):
if gu<>0  then
RETURN(gu):
fi:
od:
od:
 
 
for ORD from 2 to MAX do
for mu from 0 to MAX-ORD do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
gu:=HankelfB11(h,ORD,DEG1,DEG2,n,r,R):
if gu<>0  then
 RETURN(gu):
fi:
od:
od:
od:
 
print(`No recurrence found whose ORDER+Sum of Degrees <=`,MAX):
 
0:
end:
 
 
 
 
HankelfB:=proc():
 
if nargs=5 then
HankelfB1(args):
elif nargs=4 then
HankelfB1(args,10):
else
ERROR(`HankelfB(h,n,r,R,MAX) or HankelfB(h,n,r,R)`):
fi:
end:
 
 
 
 
 
Hankelf:=proc()
local h,n,r,N,R,MAX:
h:=args[1]:  n:=args[2]: r:=args[3]: N:=args[4]:
R:=args[5]:
 
if nargs=6 then
MAX:=args[6]:
fi:
 
if nargs=5 then
RETURN(HankelfA(h,n,r,N),HankelfB(h,n,r,R)):
elif nargs=6 then
RETURN(HankelfA(h,n,r,N,MAX),HankelfB(h,n,r,R,MAX)):
else
ERROR(`Hankelf(h,n,r,N,R) or Hankelf(h,n,r,N,R,MAX)`):
fi:
 
end:
 
 
 
 
HankelgA11:=proc(h,ORD,DEG1,DEG2,n1,n2,N1)
local ope1,i,j,c,eq,var,n1a,n2a,ope2,lu,c1,eq1,i1,i2,Godel:
 
 
if not type(ORD,integer) or ORD<0 then
ERROR(`Second argument, is the ORDER, must be a non-negative integer`):
fi:
 
if not type(DEG1,integer) or DEG1<0 then
print(`Third argument is the maximal degree of the coefficients`):
print(` in`,n1):
ERROR(`must be a non-negative integer`):
fi:
 
 
if not type(DEG2,integer) or DEG2<0 then
print(`Fourth argument is the maximal degree of the coefficients`):
print(` in`,n2):
ERROR(`must be a non-negative integer`):
fi:
 
 
 
ope1:=0:
 
var:={}:
 
for i1 from 0 to DEG1 do
for i2 from 0 to DEG2 do
for j from 0 to ORD do
ope1:=ope1+c[i1,i2,j]*n1^i1*n2^i2*N1^j:
var:=var union {c[i1,i2,j]}:
od:
od:
od:
 
eq:={}:
 
 
Godel:=trunc(evalf(sqrt((DEG1+1)*(DEG2+1)*(ORD+1))))+3:
for n1a from 0 to  Godel do
for n2a from 0 to  Godel do
lu:=0:
 
for i from 0 to ORD do
lu:=lu+subs({n1=n1a,n2=n2a},coeff(ope1,N1,i))*gHankel(h,n1a+i,n2a):
od:
eq:=eq union {lu}:
 
od:od:
 
var:=solve(eq,var):
 
if var=NULL then
RETURN(0):
fi:
 
 
c1:={}:
 
for i from 1 to nops(var) do
eq1:=op(i,var):
 
 if op(1,eq1)=op(2,eq1) then
    c1:=c1 union {op(1,eq1)}:
 fi:
od:
 
if nops(c1)>1 then
#print(`Order or Degrees too high, lower them`):
  RETURN(0):
fi:
 
if nops(c1)=0 then
RETURN(0):
fi:
 
c1:=op(1,c1):
 
ope1:=subs(var,ope1):
 
ope2:=0:
 
ope1:=normal(ope1/c1):
 
for i from 0 to ORD do
ope2:=ope2+factor(coeff(ope1,N1,i))*N1^i:
od:
 
ope2:
end:
 
 
HankelgA1:=proc(h,n,r,N,MAX)
local mu,ORD,DEG1,DEG2,gu:
 
 
 
for mu from 0 to MAX do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
 
gu:=HankelgA11(h,1,DEG1,DEG2,n,r,N):
 
if gu<>0 then
RETURN(gu):
fi:
od:
od:
 
 
for ORD from 2 to MAX do
for mu from 0 to MAX-ORD do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
gu:=HankelgA11(h,ORD,DEG1,DEG2,n,r,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
od:
od:
 
print(`No recurrence found whose ORDER+Sum of Degrees <=`,MAX):
 
0:
end:
 
 
 
 
HankelgA:=proc():
 
if nargs=5 then
HankelgA1(args):
elif nargs=4 then
HankelgA1(args,10):
else
ERROR(`HankelgA1(h,n,r,N,MAX) or HankelgA1(h,n,r,N)`):
fi:
end:
 
 
 
 
 
HankelgB11:=proc(h,ORD,DEG1,DEG2,n1,n2,N2)
local ope1,i,j,c,eq,var,n1a,n2a,ope2,lu,c1,eq1,i1,i2,Godel:
 
 
if not type(ORD,integer) or ORD<0 then
ERROR(`Second argument, is the ORDER, must be a non-negative integer`):
fi:
 
if not type(DEG1,integer) or DEG1<0 then
print(`Third argument is the maximal degree of the coefficients`):
print(` in`,n1):
ERROR(`must be a non-negative integer`):
fi:
 
 
if not type(DEG2,integer) or DEG2<0 then
print(`Fourth argument is the maximal degree of the coefficients`):
print(` in`,n2):
ERROR(`must be a non-negative integer`):
fi:
 
 
ope1:=0:
 
var:={}:
 
for i1 from 0 to DEG1 do
for i2 from 0 to DEG2 do
for j from 0 to ORD do
ope1:=ope1+c[i1,i2,j]*n1^i1*n2^i2*N2^j:
var:=var union {c[i1,i2,j]}:
od:
od:
od:
 
eq:={}:
 
 
Godel:=trunc(evalf(sqrt((DEG1+1)*(DEG2+1)*(ORD+1))))+3:
for n1a from 0 to  Godel do
for n2a from 0 to  Godel do
lu:=0:
 
for i from 0 to ORD do
lu:=lu+subs({n1=n1a,n2=n2a},coeff(ope1,N2,i))*gHankel(h,n1a,n2a+i):
od:
eq:=eq union {lu}:
 
od:od:
 
var:=solve(eq,var):
 
if var=NULL then
RETURN(0):
fi:
 
 
c1:={}:
 
for i from 1 to nops(var) do
eq1:=op(i,var):
 
 if op(1,eq1)=op(2,eq1) then
    c1:=c1 union {op(1,eq1)}:
 fi:
od:
 
if nops(c1)>1 then
#print(`Order or Degrees too high, lower them`):
  RETURN(0):
fi:
 
if nops(c1)=0 then
RETURN(0):
fi:
 
c1:=op(1,c1):
 
ope1:=subs(var,ope1):
 
ope2:=0:
 
ope1:=normal(ope1/c1):
 
for i from 0 to ORD do
ope2:=ope2+factor(coeff(ope1,N2,i))*N2^i:
od:
 
ope2:
end:
 
 
HankelgB1:=proc(h,n,r,R,MAX)
local mu,ORD,DEG1,DEG2,gu:
 
 
 
for mu from 0 to MAX do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
 
gu:=HankelgB11(h,1,DEG1,DEG2,n,r,R):
if gu<>0  then
RETURN(gu):
fi:
od:
od:
 
 
for ORD from 2 to MAX do
for mu from 0 to MAX-ORD do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
gu:=HankelgB11(h,ORD,DEG1,DEG2,n,r,R):
if gu<>0  then
 RETURN(gu):
fi:
od:
od:
od:
 
print(`No recurrence found whose ORDER+Sum of Degrees <=`,MAX):
 
0:
end:
 
 
 
 
HankelgB:=proc():
 
if nargs=5 then
HankelgB1(args):
elif nargs=4 then
HankelgB1(args,10):
else
ERROR(`HankelgB1(h,n,r,R,MAX) or HankelgB1(h,n,r,R)`):
fi:
end:
 
 
 
 
 
Hankelg:=proc()
local h,n,r,N,R,MAX:
h:=args[1]: n:=args[2]: r:=args[3]: N:=args[4]:
R:=args[5]:
 
if nargs=6 then
MAX:=args[6]:
fi:
 
if nargs=5 then
RETURN(HankelgA(h,n,r,N),HankelgB(h,n,r,R)):
elif nargs=6 then
RETURN(HankelgA(h,n,r,N,MAX),HankelgB(h,n,r,R,MAX)):
else
ERROR(`Hankelg(h,n,r,N,R) or Hankelg(h,n,r,N,R,MAX)`):
fi:
 
end:
 
 
###########Start Toeplitz
 
AToeplitz:=proc(h,i,j)
option remember:
if i=0 then
 RETURN(1):
elif i=1 then
 RETURN(h(j)):
else
RETURN(
normal(expand((AToeplitz(h,i-1,j)^2-
AToeplitz(h,i-1,j+1)*AToeplitz(h,i-1,j-1))/AToeplitz(h,i-2,j)))):
fi:
end:
 
 
fToeplitz:=proc(h,i,j):
normal(AToeplitz(h,i+1,j)/AToeplitz(h,i,j)):
end:
 
 
gToeplitz:=proc(h,i,j):
normal(AToeplitz(h,i,j+1)/AToeplitz(h,i,j)):
end:
 
 
 
 
ToeplitzfA11:=proc(h,ORD,DEG1,DEG2,n1,n2,N1)
local ope1,i,j,c,eq,var,n1a,n2a,ope2,lu,c1,eq1,i1,i2,Godel:
 
 
if not type(ORD,integer) or ORD<0 then
ERROR(`Second argument, is the ORDER, must be a non-negative integer`):
fi:
 
if not type(DEG1,integer) or DEG1<0 then
print(`Third argument is the maximal degree of the coefficients`):
print(` in`,n1):
ERROR(`must be a non-negative integer`):
fi:
 
 
if not type(DEG2,integer) or DEG2<0 then
print(`Fourth argument is the maximal degree of the coefficients`):
print(` in`,n2):
ERROR(`must be a non-negative integer`):
fi:
 
 
 
ope1:=0:
 
var:={}:
 
for i1 from 0 to DEG1 do
for i2 from 0 to DEG2 do
for j from 0 to ORD do
ope1:=ope1+c[i1,i2,j]*n1^i1*n2^i2*N1^j:
var:=var union {c[i1,i2,j]}:
od:
od:
od:
 
eq:={}:
 
 
Godel:=trunc(evalf(sqrt((DEG1+1)*(DEG2+1)*(ORD+1))))+3:
for n1a from 0 to  Godel do
for n2a from 0 to  Godel do
lu:=0:
 
for i from 0 to ORD do
lu:=lu+subs({n1=n1a,n2=n2a},coeff(ope1,N1,i))*fToeplitz(h,n1a+i,n2a):
od:
eq:=eq union {lu}:
 
od:od:
 
var:=solve(eq,var):
 
if var=NULL then
RETURN(0):
fi:
 
 
c1:={}:
 
for i from 1 to nops(var) do
eq1:=op(i,var):
 
 if op(1,eq1)=op(2,eq1) then
    c1:=c1 union {op(1,eq1)}:
 fi:
od:
 
if nops(c1)>1 then
#print(`Order or Degrees too high, lower them`):
  RETURN(0):
fi:
 
if nops(c1)=0 then
RETURN(0):
fi:
 
c1:=op(1,c1):
 
ope1:=subs(var,ope1):
 
ope2:=0:
 
ope1:=normal(ope1/c1):
 
for i from 0 to ORD do
ope2:=ope2+factor(coeff(ope1,N1,i))*N1^i:
od:
 
ope2:
end:
 
 
ToeplitzfA1:=proc(h,n,r,N,MAX)
local mu,ORD,DEG1,DEG2,gu:
 
 
 
for mu from 0 to MAX do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
 
gu:=ToeplitzfA11(h,1,DEG1,DEG2,n,r,N):
 
if gu<>0 then
RETURN(gu):
fi:
od:
od:
 
 
for ORD from 2 to MAX do
for mu from 0 to MAX-ORD do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
gu:=ToeplitzfA11(h,ORD,DEG1,DEG2,n,r,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
od:
od:
 
print(`No recurrence found whose ORDER+Sum of Degrees <=`,MAX):
 
0:
end:
 
 
 
 
ToeplitzfA:=proc():
 
if nargs=5 then
ToeplitzfA1(args):
elif nargs=4 then
ToeplitzfA1(args,10):
else
ERROR(`ToeplitzfA1(h,n,r,N,MAX) or ToeplitzfA1(h,n,r,N)`):
fi:
end:
 
 
 
 
 
ToeplitzfB11:=proc(h,ORD,DEG1,DEG2,n1,n2,N2)
local ope1,i,j,c,eq,var,n1a,n2a,ope2,lu,c1,eq1,i1,i2,Godel:
 
 
if not type(ORD,integer) or ORD<0 then
ERROR(`Second argument, is the ORDER, must be a non-negative integer`):
fi:
 
if not type(DEG1,integer) or DEG1<0 then
print(`Third argument is the maximal degree of the coefficients`):
print(` in`,n1):
ERROR(`must be a non-negative integer`):
fi:
 
 
if not type(DEG2,integer) or DEG2<0 then
print(`Fourth argument is the maximal degree of the coefficients`):
print(` in`,n2):
ERROR(`must be a non-negative integer`):
fi:
 
 
ope1:=0:
 
var:={}:
 
for i1 from 0 to DEG1 do
for i2 from 0 to DEG2 do
for j from 0 to ORD do
ope1:=ope1+c[i1,i2,j]*n1^i1*n2^i2*N2^j:
var:=var union {c[i1,i2,j]}:
od:
od:
od:
 
eq:={}:
 
 
Godel:=trunc(evalf(sqrt((DEG1+1)*(DEG2+1)*(ORD+1))))+3:
for n1a from 0 to  Godel do
for n2a from 0 to  Godel do
lu:=0:
 
for i from 0 to ORD do
lu:=lu+subs({n1=n1a,n2=n2a},coeff(ope1,N2,i))*fToeplitz(h,n1a,n2a+i):
od:
eq:=eq union {lu}:
 
od:od:
 
var:=solve(eq,var):
 
if var=NULL then
RETURN(0):
fi:
 
 
c1:={}:
 
for i from 1 to nops(var) do
eq1:=op(i,var):
 
 if op(1,eq1)=op(2,eq1) then
    c1:=c1 union {op(1,eq1)}:
 fi:
od:
 
if nops(c1)>1 then
#print(`Order or Degrees too high, lower them`):
  RETURN(0):
fi:
 
if nops(c1)=0 then
RETURN(0):
fi:
 
c1:=op(1,c1):
 
ope1:=subs(var,ope1):
 
ope2:=0:
 
ope1:=normal(ope1/c1):
 
for i from 0 to ORD do
ope2:=ope2+factor(coeff(ope1,N2,i))*N2^i:
od:
 
ope2:
end:
 
 
 
 
 
ToeplitzfB1:=proc(h,n,r,R,MAX)
local mu,ORD,DEG1,DEG2,gu:
 
 
 
for mu from 0 to MAX do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
 
gu:=ToeplitzfB11(h,1,DEG1,DEG2,n,r,R):
if gu<>0  then
RETURN(gu):
fi:
od:
od:
 
 
for ORD from 2 to MAX do
for mu from 0 to MAX-ORD do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
gu:=ToeplitzfB11(h,ORD,DEG1,DEG2,n,r,R):
if gu<>0  then
 RETURN(gu):
fi:
od:
od:
od:
 
print(`No recurrence found whose ORDER+Sum of Degrees <=`,MAX):
 
0:
end:
 
 
 
 
ToeplitzfB:=proc():
 
if nargs=5 then
ToeplitzfB1(args):
elif nargs=4 then
ToeplitzfB1(args,10):
else
ERROR(`ToeplitzfB(h,n,r,R,MAX) or ToeplitzfB(h,n,r,R)`):
fi:
end:
 
 
 
 
 
Toeplitzf:=proc()
local h,n,r,N,R,MAX:
h:=args[1]:  n:=args[2]: r:=args[3]: N:=args[4]:
R:=args[5]:
 
if nargs=6 then
MAX:=args[6]:
fi:
 
if nargs=5 then
RETURN(ToeplitzfA(h,n,r,N),ToeplitzfB(h,n,r,R)):
elif nargs=6 then
RETURN(ToeplitzfA(h,n,r,N,MAX),ToeplitzfB(h,n,r,R,MAX)):
else
ERROR(`Toeplitzf(h,n,r,N,R) or Toeplitzf(h,n,r,N,R,MAX)`):
fi:
 
end:
 
 
 
 
ToeplitzgA11:=proc(h,ORD,DEG1,DEG2,n1,n2,N1)
local ope1,i,j,c,eq,var,n1a,n2a,ope2,lu,c1,eq1,i1,i2,Godel:
 
 
if not type(ORD,integer) or ORD<0 then
ERROR(`Second argument, is the ORDER, must be a non-negative integer`):
fi:
 
if not type(DEG1,integer) or DEG1<0 then
print(`Third argument is the maximal degree of the coefficients`):
print(` in`,n1):
ERROR(`must be a non-negative integer`):
fi:
 
 
if not type(DEG2,integer) or DEG2<0 then
print(`Fourth argument is the maximal degree of the coefficients`):
print(` in`,n2):
ERROR(`must be a non-negative integer`):
fi:
 
 
ope1:=0:
 
var:={}:
 
for i1 from 0 to DEG1 do
for i2 from 0 to DEG2 do
for j from 0 to ORD do
ope1:=ope1+c[i1,i2,j]*n1^i1*n2^i2*N1^j:
var:=var union {c[i1,i2,j]}:
od:
od:
od:
 
eq:={}:
 
 
Godel:=trunc(evalf(sqrt((DEG1+1)*(DEG2+1)*(ORD+1))))+3:
for n1a from 0 to  Godel do
for n2a from 0 to  Godel do
lu:=0:
 
for i from 0 to ORD do
lu:=lu+subs({n1=n1a,n2=n2a},coeff(ope1,N1,i))*gToeplitz(h,n1a+i,n2a):
od:
eq:=eq union {lu}:
 
od:od:
 
var:=solve(eq,var):
 
if var=NULL then
RETURN(0):
fi:
 
 
c1:={}:
 
for i from 1 to nops(var) do
eq1:=op(i,var):
 
 if op(1,eq1)=op(2,eq1) then
    c1:=c1 union {op(1,eq1)}:
 fi:
od:
 
if nops(c1)>1 then
#print(`Order or Degrees too high, lower them`):
  RETURN(0):
fi:
 
if nops(c1)=0 then
RETURN(0):
fi:
 
c1:=op(1,c1):
 
ope1:=subs(var,ope1):
 
ope2:=0:
 
ope1:=normal(ope1/c1):
 
for i from 0 to ORD do
ope2:=ope2+factor(coeff(ope1,N1,i))*N1^i:
od:
 
ope2:
end:
 
 
ToeplitzgA1:=proc(h,n,r,N,MAX)
local mu,ORD,DEG1,DEG2,gu:
 
 
 
for mu from 0 to MAX do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
 
gu:=ToeplitzgA11(h,1,DEG1,DEG2,n,r,N):
 
if gu<>0 then
RETURN(gu):
fi:
od:
od:
 
 
for ORD from 2 to MAX do
for mu from 0 to MAX-ORD do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
gu:=ToeplitzgA11(h,ORD,DEG1,DEG2,n,r,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
od:
od:
 
print(`No recurrence found whose ORDER+Sum of Degrees <=`,MAX):
 
0:
end:
 
 
 
 
ToeplitzgA:=proc():
 
if nargs=5 then
ToeplitzgA1(args):
elif nargs=4 then
ToeplitzgA1(args,10):
else
ERROR(`ToeplitzgA1(h,n,r,N,MAX) or ToeplitzgA1(h,n,r,N)`):
fi:
end:
 
 
 
 
 
ToeplitzgB11:=proc(h,ORD,DEG1,DEG2,n1,n2,N2)
local ope1,i,j,c,eq,var,n1a,n2a,ope2,lu,c1,eq1,i1,i2,Godel:
 
 
if not type(ORD,integer) or ORD<0 then
ERROR(`Second argument, is the ORDER, must be a non-negative integer`):
fi:
 
if not type(DEG1,integer) or DEG1<0 then
print(`Third argument is the maximal degree of the coefficients`):
print(` in`,n1):
ERROR(`must be a non-negative integer`):
fi:
 
 
if not type(DEG2,integer) or DEG2<0 then
print(`Fourth argument is the maximal degree of the coefficients`):
print(` in`,n2):
ERROR(`must be a non-negative integer`):
fi:
 
 
ope1:=0:
 
var:={}:
 
for i1 from 0 to DEG1 do
for i2 from 0 to DEG2 do
for j from 0 to ORD do
ope1:=ope1+c[i1,i2,j]*n1^i1*n2^i2*N2^j:
var:=var union {c[i1,i2,j]}:
od:
od:
od:
 
eq:={}:
 
 
Godel:=trunc(evalf(sqrt((DEG1+1)*(DEG2+1)*(ORD+1))))+3:
for n1a from 0 to  Godel do
for n2a from 0 to  Godel do
lu:=0:
 
for i from 0 to ORD do
lu:=lu+subs({n1=n1a,n2=n2a},coeff(ope1,N2,i))*gToeplitz(h,n1a,n2a+i):
od:
eq:=eq union {lu}:
 
od:od:
 
var:=solve(eq,var):
 
if var=NULL then
RETURN(0):
fi:
 
 
c1:={}:
 
for i from 1 to nops(var) do
eq1:=op(i,var):
 
 if op(1,eq1)=op(2,eq1) then
    c1:=c1 union {op(1,eq1)}:
 fi:
od:
 
if nops(c1)>1 then
#print(`Order or Degrees too high, lower them`):
  RETURN(0):
fi:
 
if nops(c1)=0 then
RETURN(0):
fi:
 
c1:=op(1,c1):
 
ope1:=subs(var,ope1):
 
ope2:=0:
 
ope1:=normal(ope1/c1):
 
for i from 0 to ORD do
ope2:=ope2+factor(coeff(ope1,N2,i))*N2^i:
od:
 
ope2:
end:
 
 
ToeplitzgB1:=proc(h,n,r,R,MAX)
local mu,ORD,DEG1,DEG2,gu:
 
 
 
for mu from 0 to MAX do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
 
gu:=ToeplitzgB11(h,1,DEG1,DEG2,n,r,R):
if gu<>0  then
RETURN(gu):
fi:
od:
od:
 
 
for ORD from 2 to MAX do
for mu from 0 to MAX-ORD do
for DEG1 from 0 to mu do
 DEG2:=mu-DEG1:
gu:=ToeplitzgB11(h,ORD,DEG1,DEG2,n,r,R):
if gu<>0  then
 RETURN(gu):
fi:
od:
od:
od:
 
print(`No recurrence found whose ORDER+Sum of Degrees <=`,MAX):
 
0:
end:
 
 
 
 
ToeplitzgB:=proc():
 
if nargs=5 then
ToeplitzgB1(args):
elif nargs=4 then
ToeplitzgB1(args,10):
else
ERROR(`ToeplitzgB1(h,n,r,R,MAX) or ToeplitzgB1(h,n,r,R)`):
fi:
end:
 
 
 
 
 
Toeplitzg:=proc()
local h,n,r,N,R,MAX:
h:=args[1]: n:=args[2]: r:=args[3]: N:=args[4]:
R:=args[5]:
 
if nargs=6 then
MAX:=args[6]:
fi:
 
if nargs=5 then
RETURN(ToeplitzgA(h,n,r,N),ToeplitzgB(h,n,r,R)):
elif nargs=6 then
RETURN(ToeplitzgA(h,n,r,N,MAX),ToeplitzgB(h,n,r,R,MAX)):
else
ERROR(`Toeplitzg(h,n,r,N,R) or Toeplitzg(h,n,r,N,R,MAX)`):
fi:
 
end:
 
 
 

Ptor:=proc(ope1,ope2,n,r,N,R)
local gu1,gu2,f,g:
gu1:=rsolve(coeff(ope1,N,0)*f(n)+coeff(ope1,N,1)*f(n+1),f(n)):
gu1:=convert(expand(gu1),factorial):
gu2:=rsolve(coeff(ope2,R,0)*g(r)+coeff(ope2,R,1)*g(r+1),g(r)):
gu2:=convert(expand(gu2),factorial):
normal(gu1/f(0)*subs(n=0,gu2/g(0))):
end:


#HankelProofTesre(h,n,r):
#Given a procedure h (from integers to integers) tries
#to guess empirically, and then proves it rigorously,
#an explicit formula for the Hankel determinant
#A(n,r):=det(h(r+i+j)) (0<=i,j<=n-1)
# For example, try:HankelProofTerse(n->n!,n,r)
HankelProofTerse:=proc(h,n,r)
local N,R,gu,lu,m,i,j:
lu:=Hankelf(h,n,r,N,R):

if lu=0 then
 print(`Nothing found`):
 RETURN(0):
fi:

if (degree(lu[1],N)<>1 or degree(lu[2],R)<>1) then
 print(`The recurrences are of higher order then one`):
 RETURN(lu):
fi:
gu:=Ptor(lu,n,r,N,R):

print(`Proof of a Determinant Evaluation`):
print(`By Shalosh B. EKhad, c/o D. Zeilberger, Rutgers Univ. `):
print(``):
print(`Theorem: Let A(m,r) be the m by m Hankel determinant of the matrix`):
print(`whose (i,j)-entry is:`, h(r+i+j),  `( 0<=i,j<=m-1). Then `): 
print(``):
print(A(m,r)=Product(gu,n=0..m-1)):
print(``):
print(`Proof: Let's call the product on the right side B(m,r). `):
print(`The identity is trivial for m=0  and m=1 .`):
print(` By Dodgson's Rule: `):
print(A(m,r)= (A(m-1,r)*A(m-1,r+2)-A(m-1,r+1)^2)/A(m-2,r+2) ):
print(`It is completely routine to verify that `):
print(B(m,r)= (B(m-1,r)*B(m-1,r+2)-B(m-1,r+1)^2)/B(m-2,r+2) ):
print(`Hence the theorem follows by induction on m , QED `):
end:


#HankelProof(h,n,r):
#Verbose version of HankelProofTerse(h,n,r) (q.v.) e.g. type:
#HankelProof(n->n!,n,r):
HankelProof:=proc(h,n,r)
local N,R,lu1,lu2,m,i,j,f,g,vu,F,G,ku:
lu1:=Hankelf(h,n,r,N,R):
lu2:=Hankelg(h,n,r,N,R):

if lu1=0 or lu2=0 then
 print(`Nothing found`):
 RETURN(0):
fi:

if (degree(lu1[1],N)<>1 or degree(lu1[2],R)<>1) 
or  (degree(lu2[1],N)<>1 or degree(lu2[2],R)<>1) 
then
 print(`The recurrences are of higher order then one`):
 RETURN(lu1,lu2):
fi:

f:=subs(n=n-1,Ptor(lu1,n,r,N,R)):
g:=subs(r=r-1,Ptor(lu2,n,r,N,R)):


vu:=normal(expand(f*subs(n=n-1,g)-subs(r=r-1,f)*g)):
if vu<>0 then
ERROR(`Something is wrong`):
fi:

print(`Proof of a Determinant Evaluation`):
print(`By Shalosh B. EKhad, c/o D. Zeilberger, Rutgers Univ. `):
print(``):
print(`Theorem: Let A(m,r) be the m by m Hankel determinant of the matrix`):
print(`whose (i,j)-entry is:`, h(r+i+j),  `( 0<=i,j<=m-1) . Then `): 
print(``):
print(A(m,r)=Product(f,n=1..m)):
print(``):
print(`Proof: Let's call the product on the right side B(m,r). `):
print(`The identity is trivial for m=0  and m=1 .`):
print(` By Dodgson's Rule: `):
print(A(m,r)= (A(m-1,r)*A(m-1,r+2)-A(m-1,r+1)^2)/A(m-2,r+2) ):
print(`The proof would follow by induction (on m) once we prove that `):
print(B(m,r)= (B(m-1,r)*B(m-1,r+2)-B(m-1,r+1)^2)/B(m-2,r+2) ):
print(`Or equivalently: `):
print(1= (B(m-1,r)/B(m,r))*(B(m-1,r+2)/B(m-2,r+2))-
(B(m-1,r+1)/B(m,r))*(B(m-1,r+1)/B(m-2,r+2)) ):
print(`Or equivalently: `):
print(`1= (B(m-1,r)/B(m,r))*(B(m-1,r+2)/B(m-2,r+2))-`):
print(`(B(m-1,r+1)/B(m-1,r))*(B(m-1,r)/B(m,r)) `):
print(`*(B(m-1,r+1)/B(m-2,r+1))*(B(m-2,r+1)/B(m-2,r+2) ) `):
print(``):
print(`Let F(m,r):=B(m,r)/B(m-1,r), G(m,r):=B(m,r)/B(m,r-1)`):
print(`Hence, in terms of F(m,r) and G(m,r) we have to verify that `):
print(1=F(m-1,r+2)/F(m,r)-G(m-1,r+1)/F(m,r)*F(m-1,r+1)/G(m-2,r+2)):
print(``):
print(`It is readily seen that F(m,r) equals `):
print(subs(n=m,f)):
print(` and that G(m,r) equals `):
print(subs(n=m,g)):
print(`and the above identity is `):
ku:=subs({n=n-1,r=r+2},f)/f-subs({n=n-1,r=r+1},g)/f*subs({n=n-1,r=r+1},f)/
subs({n=n-2,r=r+2},g):
ku:=normal(expand(ku)):
print(evalb(ku=1)):
print(`Hence the theorem follows by induction on m , QED `):
end:


#ToeplitzProof(h,n,r):
#Toeplitz-Analog of HankelProof (q.v.) try
#ToeplitzProof(n->1/(a+n)!,n,r):
ToeplitzProof:=proc(h,n,r)
local N,R,lu1,lu2,m,i,j,f,g,vu,F,G,ku:
print(`Needs more work`):
lu1:=Toeplitzf(h,n,r,N,R):
lu2:=Toeplitzg(h,n,r,N,R):

if lu1=0 or lu2=0 then
 print(`Nothing found`):
 RETURN(0):
fi:

if (degree(lu1[1],N)<>1 or degree(lu1[2],R)<>1) 
or  (degree(lu2[1],N)<>1 or degree(lu2[2],R)<>1) 
then
 print(`The recurrences are of higher order then one`):
 RETURN(lu1,lu2):
fi:

f:=subs(n=n-1,Ptor(lu1,n,r,N,R)):
g:=subs(r=r-1,Ptor(lu2,n,r,N,R)):


vu:=normal(expand(f*subs(n=n-1,g)-subs(r=r-1,f)*g)):
if vu<>0 then
ERROR(`Something is wrong`):
fi:
RETURN(f,g,vu):

print(`Proof of a Determinant Evaluation`):
print(`By Shalosh B. EKhad, c/o D. Zeilberger, Rutgers Univ. `):
print(``):
print(`Theorem: Let A(m,r) be the m by m Hankel determinant of the matrix`):
print(`whose (i,j)-entry is:`, h(r+i+j),  `( 0<=i,j<=m-1) . Then `): 
print(``):
print(A(m,r)=Product(f,n=1..m)):
print(``):
print(`Proof: Let's call the product on the right side B(m,r). `):
print(`The identity is trivial for m=0  and m=1 .`):
print(` By Dodgson's Rule: `):
print(A(m,r)= (A(m-1,r)^2-A(m-1,r-1)*A(m-1,r+1))/A(m-2,r) ):
print(`The proof would follow by induction (on m) once we prove that `):
print(B(m,r)= (B(m-1,r)^2-B(m-1,r-1)*B(m-1,r+1))/B(m-2,r) ):
print(`Or equivalently: `):
print(1= (B(m-1,r)/B(m,r))*(B(m-1,r)/B(m-2,r))-
(B(m-1,r+1)/B(m,r))*(B(m-1,r+1)/B(m-2,r)) ):
print(`Or equivalently: `):
print(`1= (B(m-1,r)/B(m,r))*(B(m-1,r)/B(m-2,r))-`):
print(`(B(m-1,r+1)/B(m-1,r))*(B(m-1,r)/B(m,r)) `):
print(`*(B(m-1,r+1)/B(m-1,r))*(B(m-1,r)/B(m-2,r) ) `):
print(``):
print(`Let F(m,r):=B(m,r)/B(m-1,r), G(m,r):=B(m,r)/B(m,r-1)`):
print(`Hence, in terms of F(m,r) and G(m,r) we have to verify that `):
#print(1=F(m-1,r)/F(m,r)-G(m-1,r+1)/F(m,r)*F(m-1,r+1)*F(m-1,r)):
print(1=F(m-1,r)/F(m,r)*(1-G(m-1,r+1)*F(m-1,r+1))):
print(``):
print(`It is readily seen that F(m,r) equals `):
print(subs(n=m,f)):
print(` and that G(m,r) equals `):
print(subs(n=m,g)):
print(`and the above identity is `):
ku:=subs({n=n-1},f)/f*(1-subs({n=n-1,r=r+1},g)*subs({n=n-1,r=r+1},f)):
ku:=normal(expand(ku)):
RETURN(ku):
print(evalb(ku=1)):
print(`Hence the theorem follows by induction on m , QED `):
print(`Needs more work`):
end: