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


with(SolveTools):

print(`First Written: April 15, 2005: tested for Maple 8 and 9 `):
print(`Version of April 15, 2005:  `):
print():
print(`This is GuessHolo3, A Maple package`):
print(`accompanying Doron Zeilberger's  article: `):
print(`The Holonomic Ansatz: I. Foundations`):
print(`It guesses recurrence equations for discrete variables of TWO variables`):
print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg .`):
print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
 
 

ezra1:=proc()
if args=NULL then

print(` GuessHolo3: A Maple package for Handling Holnomic Functions`):
print(`of THREE discrete variables. The supporting procedures are`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(`CheckCanSmnO ,CoeffP,CoeffPs,Comm, CommSeq, DiagOper,`):
print(`EvalOpeF2 , EvalOpeF3,findrec,Findrec, FindrecF, `):
print(`GH2c, GH2dir, GH2Polc, GH2PolDir,GH2Mdir, `):
print(`GH2Mdir1, GH2MPol, GH2MPolDir, GH2MPolPar`):
print(` GR1, GR1Pol, GR1Rat, Mult, SeqFromRec`):
print(` Yafe1 `):
else
 ezra(args):
fi:
end:
 
ezra:=proc()
if args=NULL then

print(` GuessHolo3: A Maple package for empirically guessing partial recurrence`):
print(`equations satisfied by Discrete Functions of TWO Variables`):
print():
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures:  `):
print(` Asy, BaW,DiagonalStory3, GH2, GH2D,GH2DPol,  GH2Pol `):
print(`GH2PolPar,GH3, GH3Pol,GuessDiag, GuessDiagPar,GuessDiagPol,GuesssDiagPolPar,`):
print(` LatticePathStory3, LaW  `):
print(` OrthantWalk, ReciPolStory3 `):




elif nargs=1 and args[1]=Asy then
print(`Asy(ope,n,N,K): the asymptotic expansion of solutions `):
print(`to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator`):
print(`up to the K's term`):

elif nargs=1 and args[1]=BaW then
print(`BaW(Steps,a): The number of ballot-lattice walks in the k-dim lattice with pos. steps`):
print(`in the set of steps Steps ending at the point a`):
print(`For example, try: BaW({[1,0],[0,1]},[1,1]);`):


elif nargs=1 and args[1]=CanSmn then
print(`CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F): Given a holonomic representaion`):
print(`of a discrete function (let's call it F(m,n)), in terms`):
print(`of the two operators OpeM(m,n,M) and OpeN(m,n,N), expresses`):
print(`F(m+i0,n+j0) in terms of F(m+a,n+b), with 0<=a<degree(opeM,M)`):
print(`and 0<=b<=degree(opeN,N). It first goes in the m direction`):
print(`and then in the n direction. For example, try:`):
print(`CanSmn(M-1,N-1,m,n,M,N,1,1,F);`):

elif nargs=1 and args[1]=CanSnm then
print(`CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F): Given a holonomic representaion`):
print(`of a discrete function (let's call it F(m,n)), in terms`):
print(`of the two operators OpeM(m,n,M) and OpeN(m,n,N), expresses`):
print(`F(m+i0,n+j0) in terms of F(m+a,n+b), with 0<=a<degree(opeM,M)`):
print(`and 0<=b<=degree(opeN,N). It first goes in the n direction`):
print(`and then in the m direction. For example, try:`):
print(`CanSnm(M-1,N-1,m,n,M,N,1,1,F);`):


elif nargs=1 and args[1]=CheckCanSmnO then
print(`CheckCanSmnO(P,x,y,i0,j0): Checks CanSmnO on the reciprocals of P`):
print(`For example:`):
print(`CheckCanSmnO(1-x-y-x*y,x,y,3,4);`):

elif nargs=1 and args[1]=CoeffP then
print(`CoeffP(POL,xList,aList): The coeff. of xList^aList in POL`):
print(`(a polynomial in the variables of xList)`):
print(`For example, try: CoeffP(1-x-y,[x,y],[10,20]);`);


elif nargs=1 and args[1]=CoeffPs then
print(`CoeffPs(POL,xList,aList): The coeff. of xList^aList in POL`):
print(`(a SYMMETRIC polynomial in the variables of xList)`):
print(`For example, try: CoeffPs(1-x-y,[x,y],[10,20]);`);


elif nargs=1 and args[1]=Comm then
print(`Comm(Ope1,Ope2,m,n,M,N): The commutator of Ope1, Ope2`):
print(`i.e. Ope1*Ope2-Ope2*Ope1`):


elif nargs=1 and args[1]=CommSeq then
print(`CommSeq(OPE0,Ope1,m,n,M,N): Given a partial recurrence operator`):
print(`with constant coefficients, OPE0, and a partial recurrence`):
print(`operator with polynomial coefficients, finds the`):
print(`sequence of commuators with OPE0: [Ope1,[OPE0,Ope1],[OPE0,%], ...]`):
print(`and checks that the last one is a multiple of OPE0. Otherwise`):
print(`it returns FAIL`):
print(`For example, try:`):
print(`CommSeq(M-1,N-1,m,n,M,N);`):

elif nargs=1 and args[1]=DiagonalStory3 then
print(`DiagonalStory3(Steps, k,M,MaxC): Given a set of steps, and integers k`):
print(`and M`):
print(`outputs a story about the number of walks to the diagonal, the asymptotics`):
print(`and the ballot version`):
print(`For example, try: DiagonalStory({[1,0,0],[0,1,0],[0,0,1]},2,1000,10);`):

elif nargs=1 and args[1]=DiagOper then
print(`DiagOper(OpeM,OpeN,m,n,M,N):The operator in the MN direction`):
print(`deduced from OpeM(M,m,n) and OpeN(N,m,n)`):
print(`For example, try: DiagOper(M-1,N-1,m,n,M,N);`):

elif nargs=1 and args[1]=EvalOpeF2 then
print(`EvalOpeF2(F,ope,m,n,M,N,m0,n0): evaluates ope(m,n,M,N)F(m0,n0)`):
print(`For example try: `):
print(`EvalOpeF2((i,j)->binomial(i+j,j),M*N-M-N,m,n,M,N,3,5);`):


elif nargs=1 and args[1]=EvalOpeF3 then
print(`EvalOpeF3(F,ope,m,n,k,M,N,K,m0,n0,k0): evaluates ope(m,n,k,M,N,K)F(m0,n0,k0)`):
print(`For example try: `):
print(`EvalOpeF3((i,j,k)->(i+j+k)!/i!/j!/k!,M*N*K-M*N-N*K-M*K,m,n,k,M,N,K,3,5,8);`):

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):

elif nargs=1 and args[1]=GH2 then
print(`GH2(F,m,n,M,N,Patience,s1): Given a function of two variables (m,n) gives the holnomic`):
print(`representation,[Ope1(m,n,M),Ope2(m,n,N)], where M and N`):
print(`are the shift-operator in the m- and n- directions resp. `):
print(`For example, try: GH2((i,j)->i!+j!,m,n,M,N,10,15);`):

elif nargs=1 and args[1]=GH2dir then
print(`GH2dir(F,m,n,M,N): Given a function of two variables (m,n) gives the holnomic`):
print(`representation,[Ope1(m,n,M),Ope2(m,n,N)], where M and N`):
print(`are the shift-operator in the m- and n- directions resp. `):
print(`It uses the DIRECT Way`):
print(`For example, try: GH2dir((i,j)->i!+j!,m,n,M,N);`):

elif nargs=1 and args[1]=GH2c then
print(`GH2c(F,m,n,M,N,Patience,s1): Given a function of two variables (m,n) gives the holnomic`):
print(`representation,[Ope1(m,n,M),Ope2(m,n,N)], where M and N`):
print(`are the shift-operator in the m- and n- directions resp. `):
print(`In Canonical Form`):
print(`For example, try: GH2c((i,j)->i!+j!,m,n,M,N,10,15);`):

elif nargs=1 and args[1]=GH2D then
print(`GH2D(F,m,n,M,N): Given a function of two variables (m,n) guesss a recurrence`):
print(`in the diagonal direction (MN), where M and N are the shift-operators in the m-direction `):
print(`and n-direction resp.`):
print(`For example, try: GH2D((i,j)->i!+j!,m,n,M,N);`):


elif nargs=1 and args[1]=GH2M then
print(`GH2(F,m,n,M,Patience,s1):`):
print(` Given a function of two variables (m,n) guesses a recurrence`):
print(`in the m direction, where M is the shift-operator in the m-direction `):
print(`Patience is the patience, and s1 is the starting point`):
print(`For example, try: GH2M((i,j)->i!+j!,m,n,M,10,11);`):

elif nargs=1 and args[1]=GH2Mdir then
print(`GH2Mdir(F,m,n,M): Given a function of two variables (m,n) guesses a recurrence`):
print(`in the m direction, where M is the shift-operator in the m-direction `):
print(`For example, try: GH2Mdir((i,j)->i!+j!,m,n,M);`):

elif nargs=1 and args[1]=GH2Mdir1 then
print(`GH2Mdir1(F,m,n,M,deg,ORD), guesses an operator ope(m,n,M) of degree deg in (m,n)`):
print(`and order ORD in M annihilating the function `):
print(`for example, try`):
print(`GH2Mdir1((i,j)->(i+j)!/i!/j!,m,n,M,1,1);`):

elif nargs=1 and args[1]=GH2MPol then
print(`GH2MPol(POL,x,y,m,n,M,Patience,s1): The recurrence in the m-direction  of`):
print(`the coeff. of x^m*y^n of 1/POL(x,y)`):
print(`For example, try: GH2MPol(1-x-y,x,y,m,n,M,10,11);`):

elif nargs=1 and args[1]=GH2MPolDir then
print(`GH2MPolDir(POL,x,y,m,n,M,Patience): The recurrence in the m-direction  of`):
print(`the coeff. of x^m*y^n of 1/POL(x,y)`):
print(`Doing it directly`):
print(`For example, try: GH2MPolDir(1-x-y,x,y,m,n,M,10);`):

elif nargs=1 and args[1]=GH2MPolPar then
print(`GH2MPolPar(POL,x,y,m,n,M,par,Patience,s1): The recurrence in the m-direction  of`):
print(`the coeff. of x^m*y^n of 1/POL(x,y), where POL depends on par`):
print(`For example, try: GH2MPolPar(1-x-y+a*x*y,x,y,m,n,M,a,10,11);`):


elif nargs=1 and args[1]=GH2Pol then
print(`GH2Pol(POL,x,y,m,n,M,N,Patience,s1): A holonomic representation`):
print(`the coeff. of x^m*y^n of 1/POL(x,y).`):
print(`Here M and N are the shift operators in the m- and n-directions`):
print(`Patience is he pateice level (a positive integer >5)`):
print(` s1 is the starting place of investigation`):
print(`For example, try: GH2Pol(1-x-y,x,y,m,n,M,N,10,5);`):

elif nargs=1 and args[1]=GH2PolDir then
print(`GH2PolDir(POL,x,y,m,n,M,N): A holonomic representation`):
print(`the coeff. of x^m*y^n of 1/POL(x,y).`):
print(`Here M and N are the shift operators in the m- and n-directions`):
print(`Patience is he pateice level (a positive integer >5)`):
print(`For example, try: GH2PolDir(1-x-y,x,y,m,n,M,N);`):

elif nargs=1 and args[1]=GH2PolPar then
print(`GH2PolPar(POL,x,y,m,n,M,N,par,Patience,s1): A holonomic representation`):
print(`the coeff. of x^m*y^n of 1/POL(x,y).`):
print(`where POL depends on a parameter par.`):
print(`Here M and N are the shift operators in the m- and n-directions`):
print(`Patience is the patience, and s1 is the starting point`):
print(`For example, try: GH2PolPar(1-x-a*y,x,y,m,n,M,N,a,20,11);`):

elif nargs=1 and args[1]=GH2Polc then
print(`GH2Polc(POL,x,y,m,n,M,N): A holonomic representation`):
print(`the coeff. of x^m*y^n of 1/POL(x,y).`):
print(`Here M and N are the shift operators in the m- and n-directions`):
print(`In Canonical Form`):
print(`For example, try: GH2Polc(1-x-y,x,y,m,n,M,N);`):

elif nargs=1 and args[1]=GH2DPol then
print(`GH2DPol(POL,x,y,m,n,M,N): The recurrence in the diagonal of`):
print(`the coeff. of x^m*y^n of 1/POL(x,y).`):
print(`Here M and N are the shift operators in the m- and n-directions`):
print(`For example, try: GH2DPol(1-x-y,x,y,m,n,M,N);`):


elif nargs=1 and args[1]=GH3 then
print(`GH3(F,m,n,k,M,N,K,Patience,s1): `):
print(`Given a function of three variables (m,n,k) guesses a recurrence`):
print(`in the m- ,n- and k-  directions, where M,N,K, are the shift-operator in the m- `):
print(` n- and k-directions respectively.`):
print(`s1 is where it starts exploring`):
print(`For example, try: GH3((i,j,k)->(i+j+k)!/i!/j!/k!,m,n,k,M,N,K,10,11);`):

elif nargs=1 and args[1]=GH3M then
print(`GH3M(F,m,n,k,M,Patience,s1): `):
print(`Given a function of three variables (m,n,k) guesses a recurrence`):
print(`in the m direction, where M is the shift-operator in the m-direction `):
print(`s1 is where it starts exploring`):
print(`For example, try: GH3M((i,j,k)->i!+j!+k!,m,n,k,M,10,11);`):

elif nargs=1 and args[1]=GH3Pol then
print(`GH3Pol(POL,x,y,z,m,n,k,M,N,K,Patience,s1): A holonomic representation`):
print(`the coeff. of x^m*y^n*z^k of 1/POL(x,y,z).`):
print(`Here M and N are the shift operators in the m- and n-directions`):
print(`Patience is he pateice level (a positive integer >5)`):
print(` s1 is the starting place of investigation`):
print(`For example, try: GH3Pol(1-x-y-z,x,y,z,m,n,k,M,N,K,10,5);`):

elif nargs=1 and args[1]=GR1 then
print(`GR1(S1,x): guesses a polynomial from the data-set S1`):

elif nargs=1 and args[1]=GR1Pol then
print(`GR1Pol(S1,N,x): Guesses a poly in N rational in x for`):
print(`the data set S1 of the form {[i,POL(N)]}:`):

elif nargs=1 and args[1]=GR1Rat then
print(`GR1Rat(S1,N,x): Guesses a rat. function in N rational in x for`):
print(`the data set S1 of the form {[i,Rat(N)]}:`):


elif nargs=1 and args[1]=GuessDiag then
print(`GuessDiag(F,n,N,a): Given a function F of two discrete variables, guesses`):
print(`a recurrence for  the diag. F(n+a,n), where N is the shift in N`):
print(`For example, try: GuessDiag((i,j)->i!+j!,n,N,0);`):


elif nargs=1 and args[1]=GuessDiagPar then
print(`GuessDiagPar(F,n,N,a,p): Given a function of two variables (m,n) `):
print(`(featuring a paramter p), guesss a recurrence`):
print(`in the diagonal (n+a,n)`):
print(`For example, try: GuessDiagPar((i,j)->p*i!+j!,n,N,1,p);`):

elif nargs=1 and args[1]=GuessDiagPol then
print(`GuessDiagPol(P,x,y,n,N,a): Given a polynomial P of two  variables x and y, `):
print(`guesses a recurrence for  the coeffs. of x^n*y^(n+a)`):
print(`, where N is the shift in N.`):
print(`For example, try: GuessDiagPol(1-x-y+2*x*y,x,y,n,N,0);`):

elif nargs=1 and args[1]=GuessDiagPolPar then
print(`GuessDiagPolPar(P,x,y,n,N,a): Given a polynomial P of two  variables x and y, `):
print(`That also depends (rationally) on  a parameter p,`):
print(`guesses a recurrence for  the coeffs. of x^n*y^(n+a)`):
print(`, where N is the shift in N.`):
print(`For example, try: GuessDiagPolPar(1-x-y+p*x*y,x,y,n,N,0,p);`):

elif nargs=1 and args[1]=LatticePathStory3 then
print(`LatticePathStory3(Steps,Patience, s1): the story about Lattice Paths with Fundamental Steps`):
print(`Steps: For example, try: LatticePathStory3({[0,0,1],[0,1,0],[1,0,1]},10,15);`):

elif nargs=1 and args[1]=LaW then
print(`LaW(Steps,a): The number of lattice walks in the k-dim lattice with pos. steps`):
print(`in the set of steps Steps ending at the point a`):
print(`For example, try: LaW({[1,0],[0,1]},[1,1]);`):

elif nargs=1 and args[1]=Mult then
print(`Mult(Ope1,Ope2,m,n,M,N): Ope1(m,n,M,N)*Ope2(m,n,M,N)`):


elif nargs=1 and args[1]=OrthantWalk then

print(` OrthantWalk(Pt,k,Steps): the number of walks from the origin to Pt using the steps`):
print(`in Steps, staying in the nops(Pt)-dim orthant`):
print(`For example, try: OrthantWalk([0,0],{[1,1],[-1,0],[0,-1]},10);`):

elif nargs=1 and args[1]=ReciPolStory3 then
print(`ReciPolStory3(POL,x,y,z,Patience, s1): the story about the Maclaurin coefficients of the rational`);
print(`function 1/POL, where P is a polynomial in x,y,z`):
print(`For example, try: ReciPolStory3(1-x-y-z,x,y,z,10,15);`):

elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):

elif nargs=1 and args[1]=Yafe1 then
print(`Yafe1(Ope,M,N): Makes Ope nice`):
print(`For example, try: Yafe1((1/m)*M*N-1,M,N);`):

fi:
 
end:



#####One variable rational-function geussing

#Bdok1(DS1,P,x): checks that P in the single variable x  indeed fits the data set DS1
Bdok1:=proc(DaS,P,x) local i:

for i from 1 to nops(DaS) do
if subs(x=DaS[i][1],denom(P))<>0 and
normal(subs(x=DaS[i][1],P)-DaS[i][2])<>0 then

  RETURN(false):
fi:
od:
true:
end:

#GenRat1(x,dxt,dxb,a): a generic (monic top)
#rational function of degree of top dxt 
#and degree of bottom dxb
#with coeffs.parametrized by a[i] followed by the set
#of coeffs.
GenRat1:=proc(x,dxt,dxb,a) local i:

(add(a[i]*x^i,i=0..dxt))/
add(a[i+dxt+1]*x^i,i=0..dxb)
,{seq(a[i],i=0..dxt+dxb+1)}:
end:


#GR1d1d2(S1,x,d1,d2): guesses a rational function in x of degree d1/d2
#that fits the data in DS1
GR1d1d2:=proc(S1,x,d1,d2)
local eq,var,P,a,var1,i:
P:=GenRat1(x,d1,d2,a):
var:=P[2]:
P:=P[1]:

if nops(S1)<=d1+d2+3 then
 print(`Insufficient data`):
 RETURN(FAIL):
fi:

eq:={seq(numer(normal(subs(x=S1[i][1],P)-S1[i][2])),i=1..d1+d2+3)}:

var1:=Linear(eq,var):
if expand(subs(var1,numer(P)))=0 then
  RETURN(FAIL):
fi:
P:=subs(var1,P):
if Bdok1(S1,P,x) then 
   RETURN(factor(normal(P))):
 else
 RETURN(FAIL):
fi:

end:


#GR1(S1,x): guesses a polynomial from the data-set S1
GR1:=proc(S1,x) local d1,d2,gu,L,i1:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:
for L from 0 to nops(S1)-4 do
 for d1 from 0 to L do
   d2:=L-d1:
 gu:=GR1d1d2(S1,x,d1,d2):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
od:

FAIL:
end:

#GR1start(S1,x,D1,D2): guesses a polynomial from the data-set S1
#starting with degree D1/D2
GR1start:=proc(S1,x,D1,D2) local d1,d2,gu,L,i1:
if {seq(S1[i1][2],i1=1..nops(S1))}={0} then
   RETURN(0):
fi:
 for d1 from D1 to nops(S1)-3-D2 do
  for d2 from D2 to nops(S1)-3-d1 do
 gu:=GR1d1d2(S1,x,d1,d2):
  if gu<>FAIL then
    RETURN(gu):
  fi:
 od:
od:

FAIL:
end:





#GR1Pol(S1,N,x): Guesses a poly in N rational in x for
#the data set S1 of the form {[i,POL(N)]}:
GR1Pol:=proc(S1,N,x)
local gu,S1a,Deg,lu,i,i1,d1,d2:
Deg:=max(seq(degree(S1[i1][2],N),i1=1..nops(S1))):
gu:=0:

S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,0)],i1=1..nops(S1))}:
 lu:=GR1(S1a,x):
  if lu=FAIL or lu=0 then
   RETURN(FAIL):
  else
  gu:=gu+factor(normal(lu)):
  fi:


d1:=degree(numer(lu),x):
d2:=degree(denom(lu),x):
for i from 1 to Deg do
S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,i)],i1=1..nops(S1))}:
 lu:=GR1start(S1a,x,d1,d2):
  if lu=FAIL then
   RETURN(FAIL):
  else
  gu:=gu+factor(normal(lu))*N^i:
  fi:
 od:

gu:
end:


#GR1PolRat(S1,N,x,M): Guesses a poly in N rational in x for
#the data set S1 of the form {[i,POL(N,Rat(M))]}:
GR1PolRat:=proc(S1,N,x,M)
local gu,S1a,Deg,lu,i,i1:
option remember:
Deg:=max(seq(degree(S1[i1][2],N),i1=1..nops(S1))):
gu:=0:
for i from 0 to Deg do
S1a:= {seq([S1[i1][1],coeff(S1[i1][2],N,i)],i1=1..nops(S1))}:
 lu:=GR1Rat(S1a,M,x):
  if lu=FAIL then
   RETURN(FAIL):
  else
  gu:=gu+factor(normal(lu))*N^i:
  fi:
 od:
gu:
end:



####End one-variable rational function guessing




###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 1 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:





 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:

#findrecK(f,DEGREE,ORDER): Is there a unique recurrence if degree DEGREE and order ORDER
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrecK([seq(i,i=1..10)],0,2);
findrecK:=proc(f,DEGREE,ORDER)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,n,N,a:
option remember:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(false):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)=1 then
 RETURN(true):
else
 RETURN(false):
fi:

end:
 
 
#end Findrec





####Begin Asymptotics

###EndAsymptotics


#####Coeffs. of reciprocal of a polynomial


#PolToSteps(POL,xList): Inputs a polynomial POL in
#the set of variables xList and returns the 
#Weighted Steps for the scheme followed by the reciprocal
#of the constant term
PolToSteps:=proc(POL,xList) local k,i,POL1,mu,Steps,ku,ku1,i1,pu:
option remember:
POL1:=expand(POL):
k:=nops(xList):
if not type(POL,`+`) then
 print(`POL must have at least two terms`):
 RETURN(FAIL):
fi:

pu:=POL1:
for i from 1 to nops(xList) do
pu:=coeff(pu,xList[i],0):
od:

if pu=0 then
 print(`The polynomial should have a non-zero constant term`):
 RETURN(FAIL):
fi:

POL1:=expand(POL1/pu):

Steps:={}:

for i from 1 to nops(POL1) do
 mu:=op(i,POL1):

 ku:=[seq(degree(mu,xList[i1]),i1=1..nops(xList))]:
 ku1:=normal(mu/convert([seq(xList[i1]^ku[i1],i1=1..nops(xList))],`*`)):
 Steps:=Steps union {[ku,-ku1]}:
od:
Steps:=Steps minus {[[0$k],-1]}:
Steps,1/pu:
end:


#Coeff1(Steps,a): 
#The number of weighter-lattice walks in the k-dim lattice with pos. steps
#in the set of steps Steps ending at the point a
#For example, try: Coeff1({[[1,0],1],[[0,1],1]},[1,1]);
Coeff1:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

Prev:=
{seq([[seq(a[j]-Steps[i][1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}:
expand(add(Prev[i][2]*Coeff1(Steps,Prev[i][1]),i=1..nops(Prev))):

end:



#Coeff1S(Steps,a): 
#The number of weighted-lattice walks in the k-dim lattice with pos. steps
#in the set of symmetric steps Steps ending at the point a
#For example, try: Coeff1({[[1,0],1],[[0,1],1]},[1,1]);
Coeff1S:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

Prev:=
{seq([[seq(a[j]-Steps[i][
1][j],j=1..k)],Steps[i][2]],i=1..nops(Steps))}:
expand(add(Prev[i][2]*Coeff1S(Steps,sort(Prev[i][1])),i=1..nops(Prev))):

end:


#CoeffP(POL,xList,aList): The coeff. of xList^aList in POL
#(a polynomial in the variables of xList)
CoeffP:=proc(POL,xList,aList)
local Steps,coe:
Steps:=PolToSteps(POL,xList):
coe:=Steps[2]:
Steps:=Steps[1]:
coe*Coeff1(Steps,aList): 
end:

#CoeffPs(POL,xList,aList): The coeff. of xList^aList in 1/POL
#(a Symmetric polynomial in the variables of xList)
CoeffPs:=proc(POL,xList,aList)
local Steps,coe:
Steps:=PolToSteps(POL,xList):
coe:=Steps[2]:
Steps:=Steps[1]:
coe*Coeff1S(Steps,aList): 
end:

#####

#LaW(Steps,a): The number of lattice walks in the k-dim lattice with pos. steps
#in the set of steps Steps ending at the point a
#For example, try: LaW({[1,0],[0,1]},[1,1]);
LaW:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}:
add(LaW(Steps,Prev[i]),i=1..nops(Prev)):

end:

#BaW(Steps,a): The number of ballot-lattice walks in the k-dim lattice with pos. steps
#in the set of steps Steps ending at the point a
#For example, try: BaW({[1,0],[0,1]},[1,1]);
BaW:=proc(Steps,a) local i,j,k,Prev:
option remember:
k:=nops(a):
if a=[0$k] then
  RETURN(1):
fi:
if min(op(a))<0 then
  RETURN(0):
fi:

if min(seq(a[i]-a[i+1],i=1..nops(a)-1))<0 then
 RETURN(0):
fi:

Prev:={seq([seq(a[j]-Steps[i][j],j=1..k)],i=1..nops(Steps))}:
add(BaW(Steps,Prev[i]),i=1..nops(Prev)):

end:



#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)]:
od:

gu:

end:





#FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with
#poly coffs.
#e.g. try FindrecF(i->i!,n,N);
FindrecF:=proc(f,n,N)
local DEGREE, ORDER,ope,L,i,t0:
option remeber:
t0:=time():

for L from 0 do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:

 ope:=findrec([seq(f(i),i=1..(2+DEGREE)*(1+ORDER)+4)],DEGREE,ORDER,n,N):
     if ope<>FAIL then
       RETURN(ope):
      fi:
 od:
od:
FAIL:

end:


#GH2MPol(POL,x,y,m,n,M,Patience,s1): The recurrence in the m-direction  of
#the coeff. of x^m*y^n of 1/POL(x,y)
#For example, try: GH2MPol(1-x-y,x,y,m,n,M,10);
GH2MPol:=proc(P,x,y,m,n,M,Patience,s1) local i,j:
if expand(subs({x=y,y=x},P)-P)=0 then
GH2M((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M,Patience,s1):
else
GH2M((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,Patience,s1):
fi:
end:


GH2:=proc(F,m,n,M,N,Patience,s1) local i,j: GH2M(F,m,n,M,Patience,s1),GH2M((i,j)->F(j,i),n,m,N,Patience,s1):end:

GH2dir:=proc(F,m,n,M,N,MaxC) local i,j: GH2Mdir(F,m,n,M),GH2Mdir((i,j)->F(j,i),n,m,N):end:


GH2PolDir:=proc(P,x,y,m,n,M,N) local ope:
if expand(subs({x=y,y=x},P)-P)=0 then
ope:=GH2Mdir((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M):
RETURN(ope,subs({m=n,n=m,M=N},ope)):
else
GH2dir((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,N):
fi:
end:

#GH2c(F,m,n,M,N,Patience,s1): Canonical Form of GH2
GH2c:=proc(F,m,n,M,N,Patience,s1) local lu,lu1,lu2: 
lu:=GH2(F,m,n,M,N,Patience,s1):
lu1:=numer(normal(lu[1])):lu2:=numer(normal(lu[2])):
Yafe1(lu1,M,N),Yafe1(lu2,M,N):end:


GH2Polc:=proc(P,x,y,m,n,M,N)local lu:
lu:=GH2c((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,N):
Yafe1(lu[1],M,N),Yafe1(lu[2],M,N):end:




GH2DPol:=proc(P,x,y,m,n,M,N):GH2D((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,N):end:


#EvalOpeF2(F,ope,m,n,M,N,m0,n0): evaluates ope(m,n,M,N)F(m0,n0)
#For example try: 
#EvalOpeF2((i,j)->binomial(i+j,j),M*N-M-N,m,n,M,N,3,5);
EvalOpeF2:=proc(F,ope,m,n,M,N,m0,n0)
local gu,i,j:
gu:=0:
for i from 0 to degree(ope,M) do
for j from 0 to degree(ope,N) do
 gu:=gu+subs({m=m0,n=n0},coeff(coeff(ope,M,i),N,j))*F(m0+i,n0+j):
od:
od:
gu:
end:



#CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F): Given a holonomic representaion
#of a discrete function (let's call it F(m,n)), in terms
#of the two operators OpeM(m,n,M) and OpeN(m,n,N), expresses
#F(m+i0,n+j0) in terms of F(m+a,n+b), with 0<=a<degree(opeM,M)
#and 0<=b<=degree(opeN,N). It first goes in the m direction
#and then in the n direction. For example, try:
#CanSmn(M-1,N-1,m,n,M,N,1,1,F);
CanSmn:=proc(OpeM,OpeN,m,n,M,N,i0,j0,F)
local gu,Ordm,Ordn,i1,gu1,j1:
option remember:
Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):


if coeff(OpeM,M,Ordm)<>1 then
 print(OpeM, `should be monic in `, M):
 RETURN(FAIL):
fi:

if coeff(OpeN,N,Ordn)<>1 then
 print(OpeN, `should be monic in `, N):
 RETURN(FAIL):
fi:


if i0<Ordm and j0<Ordn then
 RETURN(F(m+i0,n+j0)):
fi:

if j0<Ordn then
 if i0=Ordm then
  RETURN(add(subs(n=n+j0,normal(-coeff(OpeM,M,i1)))*F(m+i1,n+j0),i1=0..Ordm-1)):
 fi:

gu:=subs(m=m+1,CanSmn(OpeM,OpeN,m,n,M,N,i0-1,j0,F)):
gu:=expand(subs(F(m+Ordm,n+j0)=CanSmn(OpeM,OpeN,m,n,M,N,Ordm,j0,F),gu)):
gu:=add(factor(normal(coeff(gu,F(m+i1,n+j0),1)))*F(m+i1,n+j0),i1=0..Ordm-1):
RETURN(gu):
fi:

#Now j0>=Ordn
if i0<Ordm then


 if j0=Ordn then
  RETURN(add(subs(m=m+i0,normal(-coeff(OpeN,N,i1)))*F(m+i0,n+i1),i1=0..Ordn-1)):
 fi:

gu:=subs(n=n+1,CanSmn(OpeM,OpeN,m,n,M,N,i0,j0-1,F)):
gu:=expand(subs(F(m+i0,n+Ordn)=CanSmn(OpeM,OpeN,m,n,M,N,i0,Ordn,F),gu)):
gu:=add(factor(normal(coeff(gu,F(m+i0,n+i1),1)))*F(m+i0,n+i1),i1=0..Ordn-1):
RETURN(gu):
fi:

#Now i0>=Ordm and j0>=Ordn

gu:=subs(m=m+1,CanSmn(OpeM,OpeN,m,n,M,N,i0-1,j0,F)):

for j1 from 0 to Ordn-1 do
gu:=expand(subs(F(m+Ordm,n+j1)=CanSmn(OpeM,OpeN,m,n,M,N,Ordm,j1,F),gu)):
od:
gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*F(m+i1,n+j1):
od:
od:
RETURN(gu1):

end:





#CanSmnO(OpeM,OpeN,m,n,M,N,i0,j0):
#CanSmnO(M-1,N-1,m,n,M,N,1,1);
CanSmnO:=proc(OpeM,OpeN,m,n,M,N,i0,j0)
local F,i1,j1,gu,gu1,Ordm,Ordn:

gu:=CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F):

Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):

gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*M^i1*N^j1:
od:
od:
gu1:=M^i0*N^j0-gu1:
RETURN(gu1):
end:



#CheckCanSmnO(P,x,y,i0,j0): Checks CanSmnO on the reciprocals of P
#For example:
#CheckCanSmnO(1-x-y-x*y,x,y,3,4);
CheckCanSmnO:=proc(P,x,y,i0,j0)
local lu,m,n,M,N,Ope,i,j,m0,n0:

lu:=GH2Pol(P,x,y,m,n,M,N,5):
Ope:=CanSmnO(lu[1],lu[2],m,n,M,N,i0,j0):

evalb(
{seq(seq(EvalOpeF2((i,j)->CoeffP(P,[x,y],[i,j]),Ope,m,n,M,N,m0,n0),m0=0..40),n0=0..40)}
={0})
:
end:




#CanSmnO1(OpeM,OpeN,m,n,M,N,i0,j0):
#CanSmnO1(M-1,N-1,m,n,M,N,1,1);
CanSmnO1:=proc(OpeM,OpeN,m,n,M,N,i0,j0)
local F,i1,j1,gu,gu1,Ordm,Ordn:

gu:=CanSmn(OpeM,OpeN,m,n,M,N,i0,j0,F):

Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):

gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*M^i1*N^j1:
od:
od:
gu1:=gu1:
RETURN(gu1):
end:






#CanSnm(OpeM,OpeN,m,n,M,N,i0,j0,F): Given a holonomic representaion
#of a discrete function (let's call it F(m,n)), in terms
#of the two operators OpeM(m,n,M) and OpeN(m,n,N), expresses
#F(m+i0,n+j0) in terms of F(m+a,n+b), with 0<=a<degree(opeM,M)
#and 0<=b<=degree(opeN,N). It first goes in the n direction
#and then in the m direction. For example, try:
#CanSnn(M-1,N-1,m,n,M,N,1,1,F);
CanSnm:=proc(OpeM,OpeN,m,n,M,N,i0,j0,F)
local gu,Ordm,Ordn,i1,gu1,j1:
option remember:
Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):


if coeff(OpeM,M,Ordm)<>1 then
 print(OpeM, `should be monic in `, M):
 RETURN(FAIL):
fi:

if coeff(OpeN,N,Ordn)<>1 then
 print(OpeN, `should be monic in `, N):
 RETURN(FAIL):
fi:


if i0<Ordm and j0<Ordn then
 RETURN(F(m+i0,n+j0)):
fi:

if j0<Ordn then
 if i0=Ordm then
  RETURN(add(subs(n=n+j0,normal(-coeff(OpeM,M,i1)))*F(m+i1,n+j0),i1=0..Ordm-1)):
 fi:

gu:=subs(m=m+1,CanSnm(OpeM,OpeN,m,n,M,N,i0-1,j0,F)):
gu:=expand(subs(F(m+Ordm,n+j0)=CanSnm(OpeM,OpeN,m,n,M,N,Ordm,j0,F),gu)):
gu:=add(factor(normal(coeff(gu,F(m+i1,n+j0),1)))*F(m+i1,n+j0),i1=0..Ordm-1):
RETURN(gu):
fi:

#Now j0>=Ordn
if i0<Ordm then


 if j0=Ordn then
  RETURN(add(subs(m=m+i0,normal(-coeff(OpeN,N,i1)))*F(m+i0,n+i1),i1=0..Ordn-1)):
 fi:

gu:=subs(n=n+1,CanSnm(OpeM,OpeN,m,n,M,N,i0,j0-1,F)):
gu:=expand(subs(F(m+i0,n+Ordn)=CanSnm(OpeM,OpeN,m,n,M,N,i0,Ordn,F),gu)):
gu:=add(factor(normal(coeff(gu,F(m+i0,n+i1),1)))*F(m+i0,n+i1),i1=0..Ordn-1):
RETURN(gu):
fi:

#Now i0>=Ordm and j0>=Ordn

gu:=subs(n=n+1,CanSmn(OpeM,OpeN,m,n,M,N,i0,j0-1,F)):


for i1 from 0 to Ordm-1 do
gu:=expand(subs(F(m+i1,n+Ordn)=CanSnm(OpeM,OpeN,m,n,M,N,i1,Ordn,F),gu)):
od:


gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*F(m+i1,n+j1):
od:
od:
RETURN(gu1):

end:


#CanSnmO(OpeM,OpeN,m,n,M,N,i0,j0):
#CanSnmO(M-1,N-1,m,n,M,N,1,1);
CanSnmO:=proc(OpeM,OpeN,m,n,M,N,i0,j0)
local F,i1,j1,gu,gu1,Ordm,Ordn:

gu:=CanSnm(OpeM,OpeN,m,n,M,N,i0,j0,F):

Ordm:=degree(OpeM,M):
Ordn:=degree(OpeN,N):

gu1:=0:

for i1 from 0 to Ordm-1 do
for j1 from 0 to Ordn-1 do
gu1:=gu1+factor(coeff(gu,F(m+i1,n+j1),1))*M^i1*N^j1:
od:
od:
gu1:=M^i0*N^j0-gu1:
RETURN(gu1):
end:


#CheckCanSnmO(P,x,y,i0,j0): Checks CanSmnO on the reciprocals of P
#For example:
#CheckCanSnmO(1-x-y-x*y,x,y,3,4);
CheckCanSnmO:=proc(P,x,y,i0,j0)
local lu,m,n,M,N,Ope,i,j,m0,n0:

lu:=GH2Pol(P,x,y,m,n,M,N,5):
Ope:=CanSnmO(lu[1],lu[2],m,n,M,N,i0,j0):

evalb(
{seq(seq(EvalOpeF2((i,j)->CoeffP(P,[x,y],[i,j]),Ope,m,n,M,N,m0,n0),m0=0..40),n0=0..40)}
={0})
:
end:

#Yafe1(Ope,M,N): Makes Ope nice
Yafe1:=proc(Ope,M,N)
local d1,d2,Ope1,i1,i2,gu:
Ope1:=numer(normal(Ope)):
d1:=degree(Ope1,M):
d2:=degree(Ope1,N):
gu:=0:
for i1 from 0 to d1 do
for i2 from 0 to d2 do
 gu:=gu+factor(coeff(coeff(Ope1,M,i1),N,i2))*M^i1*N^i2:
od:
od:

gu:

end:

#Yafe13(Ope,M,N,K): Makes Ope nice
Yafe13:=proc(Ope,M,N,K)
local d1,d2,d3,Ope1,i1,i2,i3,gu:
Ope1:=numer(normal(Ope)):
d1:=degree(Ope1,M):
d2:=degree(Ope1,N):
d3:=degree(Ope1,K):
gu:=0:
for i1 from 0 to d1 do
for i2 from 0 to d2 do
for i3 from 0 to d3 do
 gu:=gu+factor(coeff(coeff(coeff(Ope1,M,i1),N,i2),K,i3))*M^i1*N^i2*K^i3:
od:
od:
od:

gu:

end:


#Mult(Ope1,Ope2,m,n,M,N): Ope1(m,n,M,N)*Ope2(m,n,M,N)
Mult:=proc(Ope1,Ope2,m,n,M,N) local d1,d2,i1,i2,gu:
d1:=degree(Ope1,M):
d2:=degree(Ope1,N):
gu:=0:
for i1 from 0 to d1 do
for i2 from 0 to d2 do
 gu:=expand(gu+factor(coeff(coeff(Ope1,M,i1),N,i2))*M^i1*N^i2
  *subs({m=m+i1,n=n+i2},Ope2) ):
od:
od:
Yafe1(gu,M,N):
end:

#Mult3(Ope1,Ope2,m,n,k,M,N,K): Ope1(m,n,k,M,N,K)*Ope2(m,n,k,M,N,K)
Mult3:=proc(Ope1,Ope2,m,n,k,M,N,K) local d1,d2,d3,i1,i2,i3,gu:
d1:=degree(Ope1,M):
d2:=degree(Ope1,N):
d3:=degree(Ope1,K):
gu:=0:
for i1 from 0 to d1 do
for i2 from 0 to d2 do
for i3 from 0 to d3 do
 gu:=expand(gu+factor(coeff(coeff(coeff(Ope1,M,i1),N,i2),K,i3))*M^i1*N^i2*K^i3
  *subs({m=m+i1,n=n+i2,k=k+i3},Ope2) ):
od:
od:
od:

Yafe1(gu,M,N):
end:

#Comm(Ope1,Ope2,m,n,M,N): The commutator of Ope1, Ope2
#i.e. Ope1*Ope2-Ope2*Ope1
Comm:=proc(Ope1,Ope2,m,n,M,N)
Yafe1(Mult(Ope1,Ope2,m,n,M,N)-Mult(Ope2,Ope1,m,n,M,N),M,N):
end:

#Comm3(Ope1,Ope2,m,n,k,M,N,K): The commutator of Ope1, Ope2
#i.e. Ope1*Ope2-Ope2*Ope1
Comm3:=proc(Ope1,Ope2,m,n,k,M,N,K)
Yafe1(Mult3(Ope1,Ope2,m,n,k,M,N,K)-Mult(Ope2,Ope1,m,n,k,M,N,K),M,N,K):
end:


#CommSeq(OPE0,Ope1,m,n,M,N): Given a partial recurrence operator
#with constant coefficients, OPE0, and a partial recurrence
#operator with polynomial coefficients, finds the
#sequence of commuators with OPE0: [Ope1,[OPE0,Ope1],[OPE0,%], ...]
#and checks that the last one is a multiple of OPE0. Otherwise
#it returns FAIL
CommSeq:=proc(OPE0,Ope1,m,n,M,N)
local gu,lu,gu1:
lu:=Ope1:
gu:=[lu]:
while lu<>0 do
lu:=Comm(OPE0,lu,m,n,M,N):
gu:=[op(gu),lu]:
od:
gu:=[op(1..nops(gu)-1,gu)]:
gu:=[op(1..nops(gu)-1,gu),factor(gu[nops(gu)])]:
gu1:=denom(ormal(gu[nops(gu)]/OPE0));
if gu1<>1 then
 FAIL:
else
 gu:
fi:
end:

#CommSeq3(OPE0,Ope1,m,n,k,M,N,K): Given a partial recurrence operator
#with constant coefficients, OPE0, and a partial recurrence
#operator with polynomial coefficients, finds the
#sequence of commuators with OPE0: [Ope1,[OPE0,Ope1],[OPE0,%], ...]
#and checks that the last one is a multiple of OPE0. Otherwise
#it returns FAIL
CommSeq3:=proc(OPE0,Ope1,m,n,k,M,N,K)
local gu,lu,gu1:
lu:=Ope1:
gu:=[lu]:
while lu<>0 do
lu:=Comm3(OPE0,lu,m,n,k,M,N,K):
gu:=[op(gu),lu]:
od:
gu:=[op(1..nops(gu)-1,gu)]:
gu:=[op(1..nops(gu)-1,gu),factor(gu[nops(gu)])]:
gu1:=denom(ormal(gu[nops(gu)]/OPE0));
if gu1<>1 then
 FAIL:
else
 gu:
fi:
end:



Diag1:=proc(OpeM,OpeN,m,n,M,N,k0)
local gu,a,var,eq,eqTry,i1,j1,OpeD,var1T,i,var1:

OpeD:=M^k0*N^k0:
var:={}:

 gu:=CanSmnO1(OpeM,OpeN,m,n,M,N,k0,k0):
for i from 0 to k0-1 do
 gu:=gu+a[i]*CanSmnO1(OpeM,OpeN,m,n,M,N,i,i):
 var:= var union {a[i]}:
 OpeD:=OpeD+a[i]*M^i*N^i:
od:
gu:=expand(gu):
 eq:={seq(seq(coeff(coeff(gu,M,i1),N,j1),j1=0..degree(gu,N)),
i1=0..degree(gu,M))}:

eqTry:=subs({m=11/5,n=17/19},eq):

var1T:=Linear(eqTry,var):

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

var1:=Linear(eq,var):

if var1=NULL then
 RETURN(FAIL):
else
 RETURN(subs(var1,OpeD)):
fi:

end:

#DiagOper(OpeM,OpeN,m,n,M,N):The operator in the MN direction
#deduced from OpeM(M,m,n) and OpeN(N,m,n)
#For example, try: DiagOper(M-1,N-1,m,n,M,N);
DiagOper:=proc(OpeM,OpeN,m,n,M,N) local k0,gu:
for k0 from 1 do
gu:=Diag1(OpeM,OpeN,m,n,M,N,k0):
 if gu<>FAIL then
   RETURN(gu):
fi:
od:
end:




#GH2MPolPar(POL,x,y,m,n,M,par,Patience,s1): The recurrence in the m-direction  of
#the coeff. of x^m*y^n of 1/POL(x,y)
#where POL also depends (in a rational way)
#on P
#For example, try: GH2MPolPar(1-x-y-a*x*y,x,y,m,n,M,a,10,11);
GH2MPolPar:=proc(P,x,y,m,n,M,par,Patience,s1) local S1,i,Ope,Ope1,i1,j1,m0,n0,N:
S1:={}:
for i from 2 to 5 do
 S1:=S1 union {[i,GH2MPol(subs(par=i,P),x,y,m,n,M,Patience,s1)]}:
od:

for i from 6 do
Ope:= GR1Pol(S1,M,par):
if Ope<>FAIL then

  Ope1:=numer(normal(Ope)):

 if {seq(seq(normal(EvalOpeF2((i1,j1)->CoeffP(P,[x,y],[i1,j1]),
  Ope1,m,n,M,N,m0,n0)),m0=0..10),n0=0..10)}<>{0} then
print( {seq(seq(normal(EvalOpeF2((i1,j1)->CoeffP(P,[x,y],[i1,j1]),
  Ope1,m,n,M,N,m0,n0)),m0=0..10),n0=0..10)}):

  print(Ope, ` failed `):
 else
  RETURN(Ope):
 fi:
fi:
 S1:=S1 union {[i,GH2MPol(subs(par=i,P),x,y,m,n,M,Patience,s1)]}:
od:

end:

#GH2PolPar(P,x,y,m,n,M,N,par,Patience,s1): Like GH2Pol, but P depends
#on a parameter par
GH2PolPar:=proc(P,x,y,m,n,M,N,par,Patience,s1) local ope:
if expand(subs({x=y,y=x},P)-P)=0 then
ope:=GH2MPolPar(P,x,y,m,n,M,par,Patience,s1):
RETURN(ope,subs({m=n,n=m,M=N},ope)):
else
RETURN(GH2MPolPar(P,x,y,m,n,M,par,Patience,s1),
subs({m=n,n=m,M=N},GH2MPolPar(subs({x=y,y=x},P),x,y,m,n,M,par,Patience,s1)))
:
fi:
end:



#GuessDiag(F,n,N,a): Given a function F of two discrete variables, guesses
#a recurrence for  the diag. F(n+a,n), where N is the shift in N
#For example, try: GuessDiag((i,j)->i!+j!,n,N,0);
GuessDiag:=proc(F,n,N,a) local i:
FindrecF(i->F(i+a,i),n,N):
end:

#GuessDiagPol(P,x,y,n,N,a): Given a polynomial of two  variables x and y, 
#guesses a recurrence for  the coeffs. of x^n*y^(n+a)
#, where N is the shift in N
#For example, try: GuessDiagPol(1-x-y+2*x*y,x,y,n,N,0);
GuessDiagPol:=proc(P,x,y,n,N,a) local i,j:
GuessDiag((i,j)->CoeffP(P,[x,y],[i,j]),n,N,a):
end:


#GuessDiagPar(F,n,N,a,p): Given a function of two variables (m,n) 
#(featuring a paramter p), guesss a recurrence
#in the diagonal (n+a,n)
#For example, try: GuessDiagPar((i,j)->p*i!+j!,n,N,1,p);
GuessDiagPar:=proc(F,n,N,a,p)
local i,ope,h,ope1,ope2,ope3,DEG,ORD,gu1,gu2,gu3,S1,
Ope,j:

ope[1]:=FindrecF(h->subs(p=1,F(h+a,h)),n,N):
ope[2]:=FindrecF(h->subs(p=2,F(h+a,h)),n,N):
ope[3]:=FindrecF(h->subs(p=3,F(h+a,h)),n,N):
ope1:=numer(normal(ope[1])):
ope2:=numer(normal(ope[2])):
ope3:=numer(normal(ope[3])):

for i from 4 do

ope1:=numer(normal(ope[i-1])):
ope2:=numer(normal(ope[i-2])):
ope3:=numer(normal(ope[i-3])):
  if nops({degree(ope1,n),degree(ope2,n),degree(ope2,n)})=1 and
   nops({degree(ope1,N),degree(ope2,N),degree(ope2,N)})=1  then
   DEG:=degree(ope3,n):
   ORD:=degree(ope3,N):

gu1:=[seq(subs(p=i-1,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu2:=[seq(subs(p=i-2,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu3:=[seq(subs(p=i-3,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

if findrecK(gu1,DEG,ORD) and findrecK(gu2,DEG,ORD) and
    findrecK(gu3,DEG,ORD) then
      break:

  fi:
fi:

ope[i]:=FindrecF(h->subs(p=i,F(h+a,h)),n,N):

od:

S1:={}:

for j from i-2 to i+6 do

  gu1:=[seq(subs(p=j,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,n,N):
    S1:=S1 union {[j,ope1]}:

   fi:

od:

for j from i+7 do
Ope:= GR1Pol(S1,N,p):
 if Ope<>FAIL then
  RETURN(Ope):
 fi:

  gu1:=[seq(subs(p=j,F(h+a,h)),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,n,N):
    S1:=S1 union {[j,ope1]}:
   fi:

od:
end:




#GuessDiagPolPar(P,x,y,n,N,a,p): Given a polynomial of two  variables x and y, 
#also depending (rationally) on a parameter p
#guesses a recurrence for  the coeffs. of x^n*y^(n+a)
#, where N is the shift in N
#For example, try: GuessDiagPolPar(1-x-y+p*x*y,x,y,n,N,0,p);
GuessDiagPolPar:=proc(P,x,y,n,N,a,p):
GuessDiagPar((i,j)->CoeffP(P,[x,y],[i,j]),n,N,a,p):
end:




#GR1Rat(S1,N,x): Guesses a rational function in N rational in x for
#the data set S1 of the form {[i,RAT(N)]}:
GR1Rat:=proc(S1,N,x)
local Deg1,Deg2,i,i1,S1Top,S1Bot,R1,c,Top,Bot:

Deg1:=max(seq(degree(numer(S1[i1][2]),N),i1=1..nops(S1))):

Deg2:=max(seq(degree(denom(S1[i1][2]),N),i1=1..nops(S1))):

S1Top:={}:
S1Bot:={}:

for i from 1 to nops(S1) do
R1:=S1[i][2]:
 if degree(numer(R1),N)=Deg1 and degree(denom(R1),N)=Deg2  then
 c:=coeff(numer(R1),N,Deg1):
 S1Top:=S1Top union {[S1[i][1],numer(R1)/c]}:
 S1Bot:=S1Bot union {[S1[i][1],denom(R1)/c]}:
 fi:
od:
Top:=GR1Pol(S1Top,N,x):
Bot:=GR1Pol(S1Bot,N,x):

 if Top<>FAIL and Bot<>FAIL then
  RETURN(Top/Bot):
 else
  RETURN(FAIL):
 fi:

end:


#####Direct Way for GH2M: GH2Mdirect

#GenPol2(x,y,d,a): a generic Polynomial
#of dgree d in (x,y) with coeffs. indexed by a
#followed by the coefficient GenPol2(x,y,1,a)
#For example, try
GenPol2:=proc(x,y,d,a) local i,k:
convert([seq(seq(a[i,k-i]*x^i*y^(k-i),i=0..k),k=0..d)],`+`),
{seq(seq(a[i,k-i],i=0..k),k=0..d)}:
end:

#GenOper2(m,n,M,deg,ORD,a): a generic operator of degree deg in (m,n)
#and order ORD in M
GenOper2:=proc(m,n,M,deg,ORD,a) local i,var,gu,lu:
gu:=0:
var:={}:
for i from 0 to ORD do
lu:=GenPol2(m,n,deg,a[i]):
gu:=gu+lu[1]*M^i:
 var:=var union lu[2]:
od:
gu,var:
end:





#GH2Mdir1(F,m,n,M,deg,ORD), guesses an operator ope(m,n,M) of degree deg in (m,n)
#and order ORD in M annihilating the function F
#for example, try
#GH2Mdir1((i,j)->(i+j!,m,n,M,2,2);
GH2Mdir1:=proc(F,m,n,M,deg,ORD) local ope,var,eq,a,m0,n0,k,N,var1,ind1,i:
ope:=GenOper2(m,n,M,deg,ORD,a):
var:=ope[2]:
ope:=ope[1]:

eq:={}:
for k from 5 while nops(eq)<=nops(var)+5 do
for m0 from 0 to k do
 n0:=k-m0:
 eq:=eq union {EvalOpeF2(F,ope,m,n,M,N,m0+2,n0+4)}:
od:
od:

var1:=Linear(eq,var):

ind1:={}:

for i from 1 to nops(var1) do
 if op(1,var1[i])=op(2,var1[i]) then
   ind1:=ind1 union {op(1,var1[i])}:
 fi:
od:

if nops(ind1)=0 then
 RETURN(FAIL):
fi:
if nops(ind1)>1 then
 RETURN(FAIL):
fi:

ope:=subs(ind1[1]=1,subs(var1,ope)):
if {seq(seq(EvalOpeF2(F,ope,m,n,M,N,m0,n0),m0=10..20),n0=10..20)}<>{0} then
 RETURN(FAIL):
else
RETURN(Yafe1(ope,M,N)):
fi:

end:



#GH2Mdir(F,m,n,M), guesses an operator ope(m,n,M) of degree deg in (m,n)
#and order ORD in M annihilating the function F
#for example, try
#GH2Mdir((i,j)->(i+j!,m,n,M);
GH2Mdir:=proc(F,m,n,M) local ope,DEG,ORD,k,LT,i,ope1,ope2,ope3,ope4,j:

ope1:=numer(normal(FindrecF(i->F(i,20),n,N,MaxC))):
ope2:=numer(normal(FindrecF(i->F(i,21),n,N,MaxC))):
ope3:=numer(normal(FindrecF(i->F(i,22),n,N,MaxC))):

for j from 23 do
if nops({degree(ope1,N),degree(ope2,N),degree(ope3,N)})=1 and
  nops({degree(ope1,n),degree(ope2,n),degree(ope3,n)})=1 then
   DEG:=degree(ope3,n):
   ORD:=degree(ope3,N):
 break:
fi:
ope4:=numer(normal(FindrecF(i->F(i,j),n,N,MaxC))):
ope1:=ope2:
ope2:=ope3:
ope3:=ope4:
od:

 ope:=GH2Mdir1(F,m,n,M,DEG,ORD):
  if ope<>FAIL then
 LT:=coeff(ope,M,degree(ope,M)):
 ope:=add(normal(coeff(ope,M,i)/LT)*M^i,i=0..degree(ope,M)):
    RETURN(ope):
  fi:

FAIL:
end:


#GH2MPolDir(POL,x,y,m,n,M,MaxC): The recurrence in the m-direction  of
#the coeff. of x^m*y^n of 1/POL(x,y) the direct way
#For example, try: GH2MPolDir(1-x-y,x,y,m,n,M,10);
GH2MPolDir:=proc(P,x,y,m,n,M,MaxC) local i,j:
if expand(subs({x=y,y=x},P)-P)=0 then
GH2Mdir((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M,MaxC):
else
GH2Mdir((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,MaxC):
fi:
end:



#GH2D(F,m,n,M,N): Given a function of two variables (m,n) guesss a recurrence
#in the diagonal direction (MN), where M and N are the shift-operator in the m-direction 
#and n-directions 
#For example, try: GH2D((i,j)->i!+j!,m,n,M,N);
GH2D:=proc(F,m,n,M,N)
local i,ope,h,ope1,ope2,ope3,DEG,ORD,gu1,gu2,gu3,S1,
Ope,j,K,x,Ope1,m0,n0,i1:

ope[1]:=FindrecF(h->F(h+1,h),n,K):
ope[2]:=FindrecF(h->F(h+2,h),n,K):
ope[3]:=FindrecF(h->F(h+3,h),n,K):
ope1:=numer(normal(ope[1])):
ope2:=numer(normal(ope[2])):
ope3:=numer(normal(ope[3])):

for i from 4 do

ope1:=numer(normal(ope[i-1])):
ope2:=numer(normal(ope[i-2])):
ope3:=numer(normal(ope[i-3])):
  if nops({degree(ope1,n),degree(ope2,n),degree(ope2,n)})=1 and
   nops({degree(ope1,K),degree(ope2,K),degree(ope2,K)})=1  then
   DEG:=degree(ope3,n):
   ORD:=degree(ope3,K):

gu1:=[seq(F(h+i-1,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu2:=[seq(F(h+i-2,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu3:=[seq(F(h+i-3,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

if findrecK(gu1,DEG,ORD) and findrecK(gu2,DEG,ORD) and
    findrecK(gu3,DEG,ORD) then
      break:

  fi:
fi:

ope[i]:=FindrecF(h->F(h+i,h),n,K):

od:

S1:={}:

for j from i-2 to i+6 do

  gu1:=[seq(F(h+j,h),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,n,K):
    S1:=S1 union {[j,ope1]}:

   fi:

od:


for j from i+7 do
Ope:= GR1PolRat(S1,K,x,n):
 if Ope<>FAIL then

 Ope:=subs({x=m-n},Ope):
Ope:=add(factor(normal(coeff(Ope,K,i1)))*K^i1,i1=0..degree(Ope,K)):
 Ope:=subs({K=M*N},Ope):

 Ope1:=numer(normal(Ope)):
 if {seq(seq(normal(EvalOpeF2(F,Ope1,m,n,M,N,m0,n0)),m0=0..40),n0=0..40)}<>{0} then
  print(Ope, ` failed `):
 else
  RETURN(Ope):
 fi:
 fi:

  gu1:=[seq(F(h,j),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,n,K):
    S1:=S1 union {[j,ope1]}:
   fi:

od:
end:


#GH2M(F,m,n,M,Patience,s1): 
#Given a function of two variables (m,n) guesses a recurrence
#in the m direction, where M is the shift-operator in the m-direction 
#It does it the fast way (using GR1Rat instead of GR1)
#s1 is where it starts exploring
#For example, try: GH2M((i,j)->i!+j!,m,n,M,10);
GH2M:=proc(F,m,n,M,Patience,s1)
local i,ope,h,ope1,ope2,ope3,ope4,DEG,ORD,gu1,gu2,gu3,gu4,S1,Ope,j,
m0,n0,N,Ope1:
ope[1]:=FindrecF(h->F(h,s1+1),m,M):
ope[2]:=FindrecF(h->F(h,s1+2),m,M):
ope[3]:=FindrecF(h->F(h,s1+3),m,M):
ope[4]:=FindrecF(h->F(h,s1+4),m,M):
ope1:=numer(normal(ope[1])):
ope2:=numer(normal(ope[2])):
ope3:=numer(normal(ope[3])):
ope4:=numer(normal(ope[4])):



for i from 5 to Patience do

ope1:=numer(normal(ope[i-1])):
ope2:=numer(normal(ope[i-2])):
ope3:=numer(normal(ope[i-3])):
ope4:=numer(normal(ope[i-4])):


  if nops({degree(ope1,m),degree(ope2,m),degree(ope3,m),degree(ope4,m)})=1 and
   nops({degree(ope1,M),degree(ope2,M),degree(ope3,M),degree(ope4,M)})=1  then
   DEG:=degree(ope4,m):
   ORD:=degree(ope4,M):

gu1:=[seq(F(h,s1+i-1),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu2:=[seq(F(h,s1+i-2),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu3:=[seq(F(h,s1+i-3),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
gu4:=[seq(F(h,s1+i-4),h=1..(1+DEG)*(1+ORD)+5+ORD)]:


if findrecK(gu1,DEG,ORD) and findrecK(gu2,DEG,ORD) and
    findrecK(gu3,DEG,ORD)  and  findrecK(gu4,DEG,ORD)  then
      break:
  fi:
fi:
ope[i]:=FindrecF(h->F(h,i+s1),m,M):
od:

if i-1=Patience then
 print(`Couldn't do it with patience`, Patience, `and Starting place`, s1):
 RETURN(FAIL):
fi:

S1:={}:

for j from i-3 to i+6 do
  gu1:=[seq(F(h,j+s1),h=1..(1+DEG)*(1+ORD)+5+ORD)]:
   if findrecK(gu1,DEG,ORD)  then
    ope1:=findrec(gu1,DEG,ORD,m,M):
    S1:=S1 union {[j+s1,ope1]}:
   fi:
od:

if nops( {seq(S1[i][2],i=1..nops(S1))})=1 then
   RETURN(S1[1][2]):
fi:

for j from i+7 do
Ope:= GR1PolRat(S1,M,n,m):
 if Ope<>FAIL then
  Ope1:=numer(normal(Ope)):

 if {seq(seq(normal(EvalOpeF2(F,Ope1,m,n,M,N,m0,n0)),m0=0..40),n0=0..40)}<>{0} then
  print(Ope, ` failed `):
  RETURN(FAIL):
 else

   RETURN(Ope):
 fi:
fi:
  gu1:=[seq(F(h,j+s1),h=1..(1+DEG)*(1+ORD)+5+ORD)]:

   if findrecK(gu1,DEG,ORD) then
    ope1:=findrec(gu1,DEG,ORD,m,M):
    S1:=S1 union {[j+s1,ope1]}:
   fi:

od:
end:


GH2Pol:=proc(P,x,y,m,n,M,N,Patience,s1) local ope:
if expand(subs({x=y,y=x},P)-P)=0 then
ope:=GH2M((i,j)->CoeffPs(P,[x,y],[i,j]),m,n,M,Patience,s1):
RETURN(ope,subs({m=n,n=m,M=N},ope)):
else
GH2M((i,j)->CoeffP(P,[x,y],[i,j]),m,n,M,Patience,s1),
GH2M((i,j)->CoeffP(P,[x,y],[j,i]),m,n,N,Patience,s1)
fi:
end:




###################END TWO VARIABELS


###################Start Three VARIABELS



#EvalOpeF3(F,ope,m,n,k,M,N,K,m0,n0,k0): evaluates ope(m,n,k,M,N,K)F(m0,n0,k0)
#For example try: 
#EvalOpeF3((i,j,k)->(i+j+k)!/i!/j!/k!,M*N*K-M*N-N*K-M*K,m,n,k,M,N,K,3,5,8);
EvalOpeF3:=proc(F,ope,m,n,k,M,N,K,m0,n0,k0)
local gu,i,j,k1:
gu:=0:
for i from 0 to degree(ope,M) do
for j from 0 to degree(ope,N) do
for k1 from 0 to degree(ope,K) do
 gu:=gu+subs({m=m0,n=n0,k=k0},
coeff(coeff( coeff(ope,M,i)   ,N,j),K,k1))*F(m0+i,n0+j,k0+k1):
od:
od:
od:

gu:
end:

#GH3M(F,m,n,k,M,Patience,s1): 
#Given a function of three variables (m,n,k) guesses a recurrence
#in the m direction, where M is the shift-operator in the m-direction 
#s1 is where it starts exploring
#For example, try: GH3M((i,j,k)->i!+j!+k!,m,n,k,M,10);
GH3M:=proc(F,m,n,k,M,Patience,s1)
local i,j,ope,ope1,ope2,ope3,ope4,DEGm,DEGn,ORD,S1,Ope,i1,j1,
m0,n0,k0,Ope1:

ope[1]:=GH2M((i1,j1)->F(i1,j1,s1+1),m,n,M,Patience,s1):
ope[2]:=GH2M((i1,j1)->F(i1,j1,s1+2),m,n,M,Patience,s1):
ope[3]:=GH2M((i1,j1)->F(i1,j1,s1+3),m,n,M,Patience,s1):
ope[4]:=GH2M((i1,j1)->F(i1,j1,s1+4),m,n,M,Patience,s1):
ope1:=numer(normal(ope[1])):
ope2:=numer(normal(ope[2])):
ope3:=numer(normal(ope[3])):
ope4:=numer(normal(ope[4])):



for i from 5 to Patience do

ope1:=numer(normal(ope[i-1])):
ope2:=numer(normal(ope[i-2])):
ope3:=numer(normal(ope[i-3])):
ope4:=numer(normal(ope[i-4])):


  if nops({degree(ope1,m),degree(ope2,m),degree(ope3,m),degree(ope4,m)})=1 and
nops({degree(ope1,n),degree(ope2,n),degree(ope3,n),degree(ope4,n)})=1 and
   nops({degree(ope1,M),degree(ope2,M),degree(ope3,M),degree(ope4,M)})=1  then
   DEGm:=degree(ope4,m):
   DEGn:=degree(ope4,n):
   ORD:=degree(ope4,M):
   break:
fi:
ope[i]:=GH2M((i1,j1)->F(i1,j1,s1+i),m,n,M,Patience,s1):
od:

if i-1=Patience then
 print(`Couldn't do it with patience`, Patience, `and Starting place`, s1):
 RETURN(FAIL):
fi:

S1:={}:

for j from i-3 to i+6 do
ope1:=GH2M((i1,j1)->F(i1,j1,s1+j),m,n,M,Patience,s1):
    S1:=S1 union {[s1+j,ope1]}:
od:

if nops( {seq(S1[i][2],i=1..nops(S1))})=1 then
   RETURN(S1[1][2]):
fi:

for j from i+7 do

Ope:= GR1PolRat(S1,M,k,m):
 if Ope<>FAIL then
  Ope1:=numer(normal(Ope)):

 if {
seq(
seq(seq(normal(EvalOpeF3(F,Ope1,m,n,k,M,N,K,m0,n0,k0)),m0=20..30),n0=20..30),k0=20..30)
}<>{0} then
 print(Ope, ` failed `):
 RETURN(FAIL):
else

   RETURN(Ope):
 fi:
fi:

ope1:=GH2M((i1,j1)->F(i1,j1,s1+j),m,n,M,Patience,s1):
    S1:=S1 union {[j+s1,ope1]}:

od:
end:




GH3:=proc(F,m,n,k,M,N,K,Patience,s1) local i1,j1,k1: 
GH3M(F,m,n,k,M,Patience,s1),GH3M((i1,j1,k1)->F(j1,i1,k1),n,m,k,N,Patience,s1),
GH3M((i1,j1,k1)->F(k1,j1,i1),k,n,m,K,Patience,s1)
:end:


GH3Pol:=proc(P,x,y,z,m,n,k,M,N,K,Patience,s1) local ope:
if expand(subs({x=y,y=x},P)-P)=0 and
   expand(subs({z=y,y=z},P)-P)=0  then
ope:=GH3M((i,j,k)->CoeffPs(P,[x,y,z],[i,j,k]),m,n,k,M,Patience,s1):
RETURN(ope,subs({m=n,n=m,M=N},ope),subs({m=k,k=m,M=K},ope)):
else
GH3((i,j,k)->CoeffP(P,[x,y,z],[i,j,k]),m,n,k,M,N,K,Patience,s1)
fi:
end:



#OrthantWalk(Pt,k,Steps): the number of walks from the origin to Pt using the steps
#in Steps, staying in the nops(Pt)-dim orthant
#For example, try: OrthantWalk([0,0],{[1,1],[-1,0],[0,-1]},10);
OrthantWalk:=proc(Pt,Steps,k)
local Nei,n,i,Pakhot:
option remember:

Pakhot:=proc(a,b) local i: [seq(a[i]-b[i],i=1..nops(a))]:end:
n:=nops(Pt):

 if k=0 then
   if Pt=[0$n] then
     RETURN(1):
   else
     RETURN(0):
   fi:
fi:

if min(op(Pt))<0 then
   RETURN(0):
fi:

Nei:={seq(Pakhot(Pt,Steps[i]),i=1..nops(Steps))}:
add(OrthantWalk(Nei[i],Steps,k-1),i=1..nops(Nei)):

end:

####Begin Asymptotics
#Aope1(ope,n,N): the asymptotics of the difference operator with poly. coeffs. ope(n,N)
Aope1:=proc(ope,n,N)
local gu:
gu:=expand(numer(normal(ope))):
gu:=coeff(gu,n,degree(gu,n)):
gu:=expand(gu/coeff(gu,N,degree(gu,N))):
factor(gu):
end:

C3:=proc(a,b,c) option remember:

if [a,b,c]=[0,0,0] then
  RETURN(1):
fi:
if a<0 or b<0 or c<0 then
   RETURN(0):
fi:
C3(a-1,b,c)+C3(a,b-1,c)+C3(a,b,c-1)-4*C3(a-1,b-1,c-1):
end:

#Aluf1(pit): is it positive dominant?
Aluf1:=proc(pit)
local aluf,shi,i,makom:


shi:=evalf(abs(pit[1])):
aluf:={evalf(evalc(pit[1]))}:
makom:={1}:

for i from 2 to nops(pit) do

 if evalf(abs(pit[i]))>evalf(shi) then
   shi:=abs(evalf(evalc(pit[i]))):
  aluf:={pit[i]}:
  makom:={i}:
 elif evalf(abs(pit[i]))=shi then
   aluf:=aluf union {evalf(evalc(pit[i]))}:
   makom:=makom union {i}:
 fi:

od:

aluf,shi,makom:

end:



#OneStepAS1(ope1,n,N,alpha,f,S1): Given the partial asymptotic
#expansion of solutions of ope1(n,N) with exponent alpha
#extends it by one term
OneStepAS1:=proc(ope1,n,N,alpha,f,S1)
local x,f1,L,F,A,mu,i,A1:
f1:=subs(n=1/x,f):
L:=degree(f1,x):
F:=f1+A*x^(L+S1):
mu:=add(
subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha*subs(x=x/(1+i*x),F)
,i=ldegree(ope1,N)..degree(ope1,N)):
mu:=normal(mu):
mu:=taylor(mu,x=0,L+11):
mu:=simplify(mu):
for i from L+1 to L+9 while coeff(mu,x,i)=0 do od:
if i=L+10 then
 RETURN(FAIL):
fi:
mu:=coeff(mu,x,i):
A1:=subs(Linear({mu},{A}),A):
if A1=0 then
  RETURN(FAIL):
fi:

subs({x=1/n,A=A1},F):

end:


#OneStepA(ope1,n,N,alpha,f): Given the partial asymptotic
#expansion of solutions of ope1(n,N) with exponent alpha
#extends it by one term
OneStepA:=proc(ope1,n,N,alpha,f)
local S1:
for S1 from 1 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od:
OneStepAS1(ope1,n,N,alpha,f,S1):
end:

#Asy(ope,n,N,K): the asymptotic expansion of solutions 
#to ope(n,N)f(n)=0,where ope(n,N) is  a recurrence operator
#up to the K's term
Asy:=proc(ope,n,N,K)
local gu,pit,lu,makom,alpha,mu,ope1,ku,i,f,x,ka:

gu:=Aope1(ope,n,N):
if degree(gu,N)=0 then
   print(`Not of exponential growth, case not implemented`):
   RETURN(FAIL):
fi:

if not type(expand(normal(ope/gu)),`+`) then
 pit:=[solve(gu,N)]:
 makom:=Aluf1(pit)[3]:

     RETURN(pit[makom[1]]^n*expand((nops(makom)-1)!*binomial(nops(makom)-1+n,nops(makom)-1))):
fi:

pit:=[solve(gu,N)]:
makom:=Aluf1(pit)[3]:

if nops(makom)<>1 then
lu:=evalf(evalc(pit[makom[1]])):
 if coeff(lu,I,1)>(1/10)^(Digits-2) then
 print(`Not only that but Dominant roots are complex`):
 RETURN(FAIL):
fi:

fi:
lu:=evalf(evalc(pit[makom[1]])):

if coeff(lu,I,1)>(1/10)^(Digits-2) then
 print(`Dominant root is complex`):
 RETURN(FAIL):
fi:

lu:=pit[makom[1]]:

ope1:=Yafe(subs(N=lu*N,ope),N)[2]:

mu:=add(
subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)):
mu:=normal(mu):
ku:=factor(coeff(taylor(mu,x=0,nops(makom)+1),x,nops(makom))):

ka:=simplify([solve(ku,alpha)]):

alpha:=max(op(ka)):

if normal(subs(N=1,ope1))=0 then
 RETURN(lu^n*n^alpha):
fi:

f:=1:
for i from 1 to K do
f:=OneStepA(ope1,n,N,alpha,f):
od:

lu^n*n^alpha*f:
end:

#AsyC(ope,n,N,K,Ini,L): the asymptotic expansion 
#complete with the estimated constant, of solutions 
#to ope(n,N)f(n)=0, with initial conditions
#Ini, where ope(n,N) is  a recurrence operator
#up to the K's term. The constant is estimated using up to L terms
#For example, try
#AsyC(N^2-N-1,n,N,1,[1,1],100);
AsyC:=proc(ope,n,N,K,Ini,L)
local lu,gu,C,dig,i:
Digits:=100:

if nops(Ini)<degree(ope,N)-1 then
 print(`Not enough initial conditions`):
 RETURN(FAIL):
fi:

lu:=Asy(ope,n,N,K):
gu:=SeqFromRec(ope,n,N,Ini,L):

gu:=[seq(evalf(gu[i]/subs(n=i,lu)),i=1..nops(gu))]:

dig:=-evalf(ln(abs(gu[nops(gu)]-gu[nops(gu)-1]))/ln(10)):

if dig<3 then
print(`Make the last argument L, bigger`):
 RETURN(FAIL):
fi:

C:=evalf(gu[nops(gu)],trunc(dig)-3):

C*lu:
end:

###EndAsymptotics


#DiagonalStory2(Steps, k,M,MaxC): Given a set of steps, and integers k
#and M
#outputs a story about the number of walks on the diagonal, the asymptotics
#and the ballot version
DiagonalStory2:=proc(Steps,k,M,MaxC) 
local gu,i,ope,n,N,lu1,lu2,t0,gu1a,gu2a,LU,GU,L:
t0:=time():
L:=trunc((MaxC/2+3)^2+10):
gu:=[seq(LaW(Steps,[i,i]),i=1..L)]:
ope:=Findrec(gu,n,N,MaxC):

if ope=FAIL then
 print(`Make the second argument bigger`):
 RETURN(FAIL):
fi:
print(``):

print(`Theorem 1: Let a(n) be the number of lattice walks from`):
print(`(0,0) to (n,n)`):
print(`using the following Fundamental Steps:`):
print(Steps):
print(`By the way, the first few terms (starting at n=1), are`):
print([op(1..10,gu)]):
print(`a(n) satisfies the following linear recurrence: `): 
print(add(coeff(ope,N,i)*a(n+i),i=ldegree(ope,N)..degree(ope,N))=0):
print(`Skecth of proof:`):
print(`Conjecture a linear rerurrence operator with polynomial coeffs.`):
print(` of the form Q(m,n,MN) annihilating F(m,n), the number of`):
print(`walks from the origin to (m,n)`):
print(`Here M and N are the shift operators`):
print(` in the m- and n- variales: Mf(m,n):=f(m+1,n) ; Nf(m,n):=f(m,n+1).`):
print(`Let P(M,N) be the obvious linear partial recurrence operator with`):
print(` constant coeffs. annihilating F(m,n). `):
print(`Note that conjecturing Q(MN,m,n) takes`):
print(`much longer to do, hence we don't want to bother.`):
print(`Next prove the correctness of Q by`):
print(`taking succesive commutators with P`):
print(`Finally set m=n in Q(MN,m,n) and get an ordinary recurrence`):
print(`for the diagonal a(n):=F(n,n). `):
print(``):
print(`The recurrence implies, via the Birkhoff-Trijinski-Proincare method`):
print(`that the asymptotics of a(n) is:`):
lu1:=AsyC(ope,n,N,k,[op(1..degree(ope,N),gu)],M):
gu1a:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],M):
print(lu1):
print(`And in pure floating-point`, evalf(lu1)):


print(`----------------------------`):


gu:=[seq(BaW(Steps,[i,i]),i=1..L)]:
ope:=Findrec(gu,n,N,MaxC):

if ope=FAIL then
 RETURN(FAIL):
fi:

print(`Theorem 2: Let b(n) be the number of lattice walks from`):
print(`(0,0) to (n,n)`):
print(`using the following Fundamental Steps:`):
print(Steps):
print(`Staying in the region m>=n>0 (i.e. ballot walks). `):
print(`By the way, the first few terms (starting at n=1), are`):
print([op(1..10,gu)]):
print(`b(n) satisfies the following linear recurrence `): 
print(add(coeff(ope,N,i)*b(n+i),i=ldegree(ope,N)..degree(ope,N))=0):
print(`Skecth of proof: Analogous to the above.`):
print(``):
print(`This implies, via the Birkhoff-Trijinski-Proincare method, that the`):
print(`the  asymptotics of b(n) is: `):
lu2:=AsyC(ope,n,N,k,[op(1..degree(ope,N),gu)],M):
gu2a:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],M):
print(lu2):
print(`And in floating-point`, evalf(lu2)):
print(`By dividing the asymptotic expression of b(n) by that of a(n)`):
print(` we get the following`):
print(`Corollary: The asymptotics of the probabibility of a walk `):
print(`that winded up at (n,n) staying in m>=n is:`):
LU:=normal(lu2/lu1):
LU:=normal(subs(n=1/n,LU)):
LU:=series(LU,n=0,k+1):
LU:=add(coeff(LU,n,i)*n^i,i=0..k):
LU:=subs(n=1/n,LU):
print(LU+O(1/n^(k+1)) ):
GU:=[seq(evalf(gu2a[i]/gu1a[i]/subs(n=i,LU)),i=nops(gu1a)-9..nops(gu1a))]:
print(`----------------------------`):
print(`Just as a check, the last 10 terms of the sequence of ratios of `):
print(`the probabilities to the above asymptotics is:`):
print(GU):
if abs(1-GU[nops(GU)])<1/100 then
 print(`as you can see, they are pretty close to 1`):
else
 print(`this is rather disappointing, there must be something wrong.`):
fi:
print(`----------------------------`):

print(`The whole thing took `, time()-t0, `seconds of CPU time.`):
end:





#LatticePathStory2(Steps,Patience, s1): the story about Lattice Paths with Fundamental Steps
#Steps: For example, try: LatticePathStory({[0,1],[1,0],[1,1]},10,15);
LatticePathStory2:=proc(Steps, Patience, s1) 
local gu,m,n,M,N,i,j,opeM,opeN,P,mu,ku,t0:
t0:=time():
gu:=GH2((i,j)->LaW(Steps,[i,j]),m,n,M,N,Patience,s1) :

if gu=FAIL or gu[1]=FAIL or gu[2]=FAIL then
print(`Make the second and/or third arguments larger`):
 RETURN(FAIL):
fi:

opeM:=gu[1]: opeN:=gu[2]:

print(`Theorem: Let F(m,n) be the number of lattice paths from (0,0) to (m,n)`):
print(`in the positive quadrant of the two-dimensional square lattice`):
print(`where the Fundamental Steps are`):
print(Steps):
print(`Let M and N be the shift operators in the m and n directions, respectively`):
print(` (M f(m,n):=f(m+1,n), Nf(m,n):=f(m,n+1) ) `):
print(` Then F(m,n) satisfies the following pure recurrences`):
print(add(coeff(opeM,M,i)*F(m+i,n),i=ldegree(opeM,M)..degree(opeM,M))=0):
print(``):
print(add(coeff(opeN,N,i)*F(m,n+i),i=ldegree(opeN,N)..degree(opeN,N))=0):

print(`Proof: We will only do the first recurrence`):
print(`We have to prove that F(m,n) is annihilated by the operator, let's call it Q:=`):
print(opeM):
opeM:=numer(normal(opeM)):
print(`or equivalenty`):
opeM:=numer(normal(opeM)):	
print(opeM):
print(`We know, by the obvious combinatorics, that F(m,n) is annihilated by`):
P:=1-add(1/M^Steps[i][1]/N^Steps[i][2],i=1..nops(Steps)):
P:=numer(normal(P)):
print(P):
mu:=CommSeq(P,opeM,m,n,M,N):

print(`The sequence of successive commutators is`):
print([op(2..nops(mu),mu)]):

ku:=normal(mu[nops(mu)]/P):

if degree(denom(ku),{M,N})<>0 then
  print(`Something is wrong with the proof`):
else
print(`As you can see, the last entry`, mu[nops(mu)], `is a multiple of `, P):
print(`The proof follows by backwards induction, after checking the boundary conditions`):
print(` QED .`):
fi:

print(`------------------------------------------------------`):

print(`The whole thing took `, time()-t0, `seconds of CPU time`):
end:





#ReciPolStory2(POL,x,y,Patience, s1): the story about the Maclaurin coefficients of the rational
#function 1/POL, where P is a polynomial in x,y
#For example, try: ReciPolStory(1-x-y,x,y,10,15);
ReciPolStory2:=proc(POL,x,y, Patience, s1) 
local gu,m,n,M,N,i,j,opeM,opeN,P,mu,ku,t0:
t0:=time():
gu:=GH2Pol(POL,x,y,m,n,M,N,Patience,s1):
if gu=FAIL or gu[1]=FAIL or gu[2]=FAIL then
print(`Make the second and/or third arguments larger`):
 RETURN(FAIL):
fi:

opeM:=gu[1]: opeN:=gu[2]:

print(`Theorem: Let F(m,n) be the coefficient of x^m*y^n in the Maclaurin series of`):
print(1/POL):
print(`F(m,n) satisfies the following pure recurrences`):
print(add(coeff(opeM,M,i)*F(m+i,n),i=ldegree(opeM,M)..degree(opeM,M))=0):
print(``):
print(add(coeff(opeN,N,i)*F(m,n+i),i=ldegree(opeN,N)..degree(opeN,N))=0):

print(`Proof: We will only do the first recurrence`):
print(`Let M and N be the shift operators in the m and n directions, respectively`):
print(` (M f(m,n):=f(m+1,n), Nf(m,n):=f(m,n+1) ) `):

print(`We have to prove that F(m,n) is annihilated by the operator, let's call it Q:=`):
print(opeM):
opeM:=numer(normal(opeM)):
print(`or equivalenty`):
opeM:=numer(normal(opeM)):	
print(opeM):
print(`We know, by the obvious algebra, that F(m,n) is annihilated by`):
P:=subs({x=1/M,y=1/N},POL):
P:=numer(normal(P)):
print(P):
mu:=CommSeq(P,opeM,m,n,M,N):

print(`The sequence of successive commutators is`):
print([op(2..nops(mu),mu)]):

ku:=normal(mu[nops(mu)]/P):

if degree(denom(ku),{M,N})<>0 then
  print(`Something is wrong with the proof`):
else
print(`As you can see, the last entry`, mu[nops(mu)], `is a multiple of `, P):
print(`The proof follows by backwards induction, after checking the boundary conditions`):
print(` QED .`):
fi:

print(`--------------------------------------------------------------------------------`):

print(`The whole thing took `, time()-t0, `seconds of CPU time`):
end:



#DiagonalStory3(Steps, k,M,MaxC): Given a set of steps, and integers k
#and M
#outputs a story about the number of walks to the diagonal, the asymptotics
#and the ballot version
DiagonalStory3:=proc(Steps,k,M,MaxC) 
local gu,i,ope,n,N,lu1,lu2,t0,gu1a,gu2a,LU,GU,L:
t0:=time():
L:=trunc((MaxC/2+3)^2+10):
gu:=[seq(LaW(Steps,[i,i,i]),i=1..L)]:
ope:=Findrec(gu,n,N,MaxC):

if ope=FAIL then
 print(`Make the second argument bigger`):
 RETURN(FAIL):
fi:
print(``):

print(`Theorem 1: Let a(n) be the number of lattice walks from`):
print(`(0,0,0) to (n,n,n)`):
print(`using the following Fundamental Steps:`):
print(Steps):
print(`By the way, the first few terms (starting at n=1), are`):
print([op(1..10,gu)]):
print(`a(n) satisfies the following linear recurrence: `): 
print(add(coeff(ope,N,i)*a(n+i),i=ldegree(ope,N)..degree(ope,N))=0):
print(`Skecth of proof:`):
print(`Conjecture a linear rerurrence operator with polynomial coeffs.`):
print(` of the form Q(m,n,MN) annihilating F(m,n), the number of`):
print(`walks from the origin to (m,n)`):
print(`Here M and N are the shift operators`):
print(` in the m- and n- variales: Mf(m,n):=f(m+1,n) ; Nf(m,n):=f(m,n+1).`):
print(`Let P(M,N) be the obvious linear partial recurrence operator with`):
print(` constant coeffs. annihilating F(m,n). `):
print(`Note that conjecturing Q(MN,m,n) takes`):
print(`much longer to do, hence we don't want to bother.`):
print(`Next prove the correctness of Q by`):
print(`taking succesive commutators with P`):
print(`Finally set m=n in Q(MN,m,n) and get an ordinary recurrence`):
print(`for the diagonal a(n):=F(n,n). `):
print(``):
print(`The recurrence implies, via the Birkhoff-Trijinski-Proincare method`):
print(`that the asymptotics of a(n) is:`):
lu1:=AsyC(ope,n,N,k,[op(1..degree(ope,N),gu)],M):
gu1a:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],M):
print(lu1):
print(`And in pure floating-point`, evalf(lu1)):


print(`----------------------------`):


gu:=[seq(BaW(Steps,[i,i,i]),i=1..L)]:
ope:=Findrec(gu,n,N,MaxC):

if ope=FAIL then
 RETURN(FAIL):
fi:

print(`Theorem 2: Let b(n) be the number of lattice walks from`):
print(`(0,0,0) to (n,n,n)`):
print(`using the following Fundamental Steps:`):
print(Steps):
print(`Staying in the region m>=n>0 (i.e. ballot walks). `):
print(`By the way, the first few terms (starting at n=1), are`):
print([op(1..10,gu)]):
print(`b(n) satisfies the following linear recurrence `): 
print(add(coeff(ope,N,i)*b(n+i),i=ldegree(ope,N)..degree(ope,N))=0):
print(`Skecth of proof: Analogous to the above.`):
print(``):
print(`This implies, via the Birkhoff-Trijinski-Proincare method, that the`):
print(`the  asymptotics of b(n) is: `):
lu2:=AsyC(ope,n,N,k,[op(1..degree(ope,N),gu)],M):
gu2a:=SeqFromRec(ope,n,N,[op(1..degree(ope,N),gu)],M):
print(lu2):
print(`And in floating-point`, evalf(lu2)):
print(`By dividing the asymptotic expression of b(n) by that of a(n)`):
print(` we get the following`):
print(`Corollary: The asymptotics of the probabibility of a walk `):
print(`that winded up at (n,n) staying in m>=n is:`):
LU:=normal(lu2/lu1):
LU:=normal(subs(n=1/n,LU)):
LU:=series(LU,n=0,k+1):
LU:=add(coeff(LU,n,i)*n^i,i=0..k):
LU:=subs(n=1/n,LU):
print(LU+O(1/n^(k+1)) ):
GU:=[seq(evalf(gu2a[i]/gu1a[i]/subs(n=i,LU)),i=nops(gu1a)-9..nops(gu1a))]:
print(`----------------------------`):
print(`Just as a check, the last 10 terms of the sequence of ratios of `):
print(`the probabilities to the above asymptotics is:`):
print(GU):
if abs(1-GU[nops(GU)])<1/100 then
 print(`as you can see, they are pretty close to 1`):
else
 print(`this is rather disappointing, there must be something wrong.`):
fi:
print(`----------------------------`):

print(`The whole thing took `, time()-t0, `seconds of CPU time.`):
end:

#LatticePathStory3(Steps,Patience, s1): the story about Lattice Paths with Fundamental Steps
#Steps: For example, try: LatticePathStory({[0,1],[1,0],[1,1]},10,15);
LatticePathStory3:=proc(Steps, Patience, s1) 
local gu,m,n,M,N,K,i,j,k,opeM,opeN,opeK,P,mu,ku,t0:
t0:=time():
gu:=GH3((i,j,k)->LaW(Steps,[i,j,k]),m,n,k,M,N,K,Patience,s1) :

if gu=FAIL or gu[1]=FAIL or gu[2]=FAIL or gu[3]=FAIL then
print(`Make the second and/or third arguments larger`):
 RETURN(FAIL):
fi:

opeM:=gu[1]: opeN:=gu[2]: opeK:=gu[3]:

print(`Theorem: Let F(m,n,k) be the number of lattice paths from (0,0,0) to (m,n,k)`):
print(`in the positive quadrant of the three-dimensional cubic lattice`):
print(`where the Fundamental Steps are`):
print(Steps):
print(`Let M, N and K be the shift operators in the m, n and k directions, respectively`):
print(` (M f(m,n,k):=f(m+1,n,k), Nf(m,n,k):=f(m,n+1,k) , Kf(m,n,k):=f(m,n,k+1) .) `):
print(` Then F(m,n,k) satisfies the following pure recurrences`):
print(add(coeff(opeM,M,i)*F(m+i,n,k),i=ldegree(opeM,M)..degree(opeM,M))=0):
print(``):
print(add(coeff(opeN,N,i)*F(m,n+i,k),i=ldegree(opeN,N)..degree(opeN,N))=0):
print(``):
print(add(coeff(opeK,K,i)*F(m,n,k+i),i=ldegree(opeK,K)..degree(opeK,K))=0):

print(`Proof: We will only do the first recurrence`):
print(`We have to prove that F(m,n) is annihilated by the operator, let's call it Q:=`):
print(opeM):
opeM:=numer(normal(opeM)):
print(`or equivalenty`):
opeM:=numer(normal(opeM)):	
print(opeM):
print(`We know, by the obvious combinatorics, that F(m,n,k) is annihilated by`):
P:=1-add(1/M^Steps[i][1]/N^Steps[i][2]/K^Steps[i][3],i=1..nops(Steps)):
P:=numer(normal(P)):
print(P):
mu:=CommSeq3(P,opeM,m,n,k,M,N,K):

print(`The sequence of successive commutators is`):
print([op(2..nops(mu),mu)]):

ku:=normal(mu[nops(mu)]/P):

if degree(denom(ku),{M,N,K})<>0 then
  print(`Something is wrong with the proof`):
else
print(`As you can see, the last entry`, mu[nops(mu)], `is a multiple of `, P):
print(`The proof follows by backwards induction, after checking the boundary conditions`):
print(` QED .`):
fi:

print(`------------------------------------------------------`):

print(`The whole thing took `, time()-t0, `seconds of CPU time`):
end:

#ReciPolStory3(POL,x,y,z,Patience, s1): the story about the Maclaurin coefficients of the rational
#function 1/POL, where P is a polynomial in x,y,z
#For example, try: ReciPolStory3(1-x-y-z,x,y,z,10,15);
ReciPolStory3:=proc(POL,x,y,z, Patience, s1) 
local gu,m,n,M,N,K,i,j,k,opeM,opeN,opeK,P,mu,ku,t0:
t0:=time():
gu:=GH3Pol(POL,x,y,z,m,n,k,M,N,K,Patience,s1):
if gu=FAIL or gu[1]=FAIL or gu[2]=FAIL then
print(`Make the second and/or third arguments larger`):
 RETURN(FAIL):
fi:

opeM:=gu[1]: opeN:=gu[2]: opeK:=gu[3]:

print(`Theorem: Let F(m,n,k) be the coefficient of x^m*y^n*z^k in the Maclaurin series of`):
print(1/POL):
print(`F(m,n,k) satisfies the following pure recurrences`):
print(add(coeff(opeM,M,i)*F(m+i,n,k),i=ldegree(opeM,M)..degree(opeM,M))=0):
print(``):
print(add(coeff(opeN,N,i)*F(m,n+i,k),i=ldegree(opeN,N)..degree(opeN,N))=0):
print(``):
print(add(coeff(opeK,K,i)*F(m,n,k+i),i=ldegree(opeK,K)..degree(opeK,K))=0):

print(`Proof: We will only do the first recurrence`):
print(`Let M and N be the shift operators in the m and n directions, respectively`):
print(` (M f(m,n,k):=f(m+1,n,k), Nf(m,n):=f(m,n+1,k), Kf(m,n,k):=f(m,n,k+1) ) `):

print(`We have to prove that F(m,n,k) is annihilated by the operator, let's call it Q:=`):
print(opeM):
opeM:=numer(normal(opeM)):
print(`or equivalenty`):
opeM:=numer(normal(opeM)):	
print(opeM):
print(`We know, by the obvious algebra, that F(m,n,k) is annihilated by`):
P:=subs({x=1/M,y=1/N,z=1/K},POL):
P:=numer(normal(P)):
print(P):
mu:=CommSeq3(P,opeM,m,n,k,M,N,K):

print(`The sequence of successive commutators is`):
print([op(2..nops(mu),mu)]):

ku:=normal(mu[nops(mu)]/P):

if degree(denom(ku),{M,N,K})<>0 then
  print(`Something is wrong with the proof`):
else
print(`As you can see, the last entry`, mu[nops(mu)], `is a multiple of `, P):
print(`The proof follows by backwards induction, after checking the boundary conditions`):
print(` QED .`):
fi:

print(`--------------------------------------------------------------------------------`):

print(`The whole thing took `, time()-t0, `seconds of CPU time`):
end: