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


print(`First Written: April 17, 2007: tested for Maple 10 `):
print(`Version of April 17, 2006:  `):
print():
print(`This is QuarterPlane, A Maple package`):
print(`To study walks on the Quarter Plane`):
print(`accompanying Manuel Kauers and Doron Zeilberger's  article: `):
print(`"The Quasi-Holonomic Ansatz and Restricted Lattice Walks" `):

print():
print(`The most current version is available on WWW at:`):
print(` http://www.math.rutgers.edu/~zeilberg/tokhniot/QuarterPlane .`):
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(`The supporting procedures are: Findrec, GenPureOper,`):
print(`GenOper, GuessOper1, GuessPureOper1 `):

else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` QuarterPlane: A Maple package for studying `):
print(`walks with arbitrary steps, confined to the quarter plane`):

print(`The MAIN procedures are`):
print(` GuessRec0, W  `):




elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): `):
print(`Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. 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]=GenPureOper then
print(`GenPureOper(n,a,b,N,d1): a recurrence operator `):
print(`of the form P(n,a,b,N)`):
print(`of  degree d1`):
print(`For example, try:GenPureOper(n,a,b,N,d1);`):

elif nargs=1 and args[1]=GenOper then
print(`GenOper(n,a,b,N,A,B,d1,d2): a generic recurrence operator `):
print(`of the form P(n,N)+aQ1(n,a,b,N,A,B)+bQ2(n,a,b,N,A,B)`):
print(`with degree of P d1, and degrees of Q1 and Q2 d2`):
print(`followed by the list of coefficients `):
print(`For example, try:GenOper(n,a,b,N,A,B,1,1);`):

elif nargs=1 and args[1]=GuessOper1 then
print(`GuessOper1(n,a,b,N,A,B,d1,d2,Steps): `):
print(`guesses a a recurrence operator `):
print(`of the form P(n,N)+aQ1(n,a,b,N,A,B)+bQ2(n,a,b,N,A,B)`):
print(`with degree of P d1, and degrees of Q1 and Q2 d2`):
print(`annihilating F(n,a,b):= The number of n-step walks `):
print(`from the origin to [a,b] staying in`):
print(`the quarter-plane using the steps in the set of Steps Steps`):
print(`For example, try:`):
print(`GuessOper1(n,a,b,N,A,B,1,1,{[0,1],[1,0],[0,-1],[-1,0]});`):

elif nargs=1 and args[1]=GuessPureOper1 then
print(`GuessPureOper1(n,a,b,N,d1,Steps): `):
print(`guesses a a recurrence operator `):
print(`of the form P(n,a,b,N)`):
print(`of total degree d1 `):
print(`annihilating F(n,a,b):= The number of n-step walks `):
print(`from the origin to [a,b] staying in`):
print(`the quarter-plane using the steps in the set of Steps Steps`):
print(`For example, try:`):
print(`GuessPureOper1(n,a,b,N,1,{[0,1],[1,0],[0,-1],[-1,0]});`):

elif nargs=1 and args[1]=GuessRec0 then
print(`GuessRec0(Steps,a,b,n,N,MaxC): guesses a recurrence operator`):
print(`ope(n,N) of complexity <=MaxC`):
print(`for the sequence: number of ways of going from the origin`):
print(`to [a,b] using data for n<=L, for example, try:`):
print(`GuessRec0({[1,1],[-1,0],[0,-1]},0,0,n,N,4);`):

elif nargs=1 and args[1]=W then
print(`W(n,a,b,Steps): the number of n-step walks in the positive `):
print(`discrete (square-lattice) quarter-plane staring at the`):
print(`origin and ending at [a,b] using the set of steps in Steps`):
print(`For example, try: W(3,2,1,{[1,0],[0,1]});`):
print(`For Kreweras, try: W(9,0,0,{[-1,0],[0,-1],[1,1]});`):
print(`For Gessel, try: W(8,0,0,{[-1,0],[1,0],[-1,-1],[1,1]});`):
else
 print(`There is no such thing as`, args):

fi:

end:


################Start Findrec

with(linalg):
#with(SolveTools);
#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:


#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,i,n1:

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
    if {seq(expand(add(subs(n=n1,coeff(ope,N,i))*f[n1+i],i=0..ORDER)),
           n1=1..nops(f)-ORDER)}={0} then
               RETURN(ope):
  else
 print({seq(expand(add(subs(n=n1,coeff(ope,N,i))*f[n1+i],i=0..ORDER)),
           n1=1..nops(f)-ORDER)}):

      fi:
   fi:
 od:
od:
FAIL:

end:


################End Findrec and its offsprings



#W(n,a,b,Steps): the number of n-step walks in the positive 
#discrete (square-lattice) quarter-plane staring at the
#origin and ending at [a,b] using the set of steps in Steps
#For example, try: W(3,2,1,{[1,0],[0,1]});
W:=proc(n,a,b,Steps) local i:
option remember:
if n<0 then
  RETURN(0):
fi:

if n=0 then
  if a=0 and b=0 then
     RETURN(1):
  else
     RETURN(0):
  fi:
fi:

if a<0 or b<0 then
  RETURN(0):
fi:

add(W(n-1,a-Steps[i][1],b-Steps[i][2],Steps),i=1..nops(Steps)):
end:



#GenOper(n,a,b,N,A,B,d1,d2): a recurrence operator 
#of the form P(n,N)+aQ1(n,a,b,N,A,B)+bQ2(n,a,b,N,A,B)
#with degree of P d1, and degrees of Q1 and Q2 d2
#For example, try:GenOper(n,a,b,N,A,B,1,1);
GenOper:=proc(n,a,b,N,A,B,d1,d2) local ope,var,i1,i2,i3,i4,i5,i6,c1,c2,c3:
ope:=0:
var:={}:
for i1 from 0 to d1 do
for i2 from 0 to d1-i1 do
 ope:=ope+c1[i1,i2]*n^i1*N^i2:
 var:=var union {c1[i1,i2]}:
od:
od:




for i1 from 0 to d2 do
for i2 from 0 to d2-i1 do
for i3 from 0 to d2-i1-i2 do
for i4 from 0 to d2-i1-i2-i3 do
for i5 from 0 to d2-i1-i2-i3-i4 do
for i6 from 0 to d2-i1-i2-i3-i4-i5 do
 ope:=ope+a*c2[i1,i2,i3,i4,i5,i6]*n^i1*N^i2*a^i3*b^i4*A^i5*B^i6:
 var:=var union {c2[i1,i2,i3,i4,i5,i6]}:
od:
od:
od:
od:
od:
od:

for i1 from 0 to d2 do
for i2 from 0 to d2-i1 do
for i3 from 0 to d2-i1-i2 do
for i4 from 0 to d2-i1-i2-i3 do
for i5 from 0 to d2-i1-i2-i3-i4 do
for i6 from 0 to d2-i1-i2-i3-i4-i5 do
 ope:=ope+b*c3[i1,i2,i3,i4,i5,i6]*n^i1*N^i2*a^i3*b^i4*A^i5*B^i6:
 var:=var union {c3[i1,i2,i3,i4,i5,i6]}:
od:
od:
od:
od:
od:
od:
ope,var:
end:



#GuessOper1(n,a,b,N,A,B,d1,d2,Steps): 
#guesses a a recurrence operator 
#of the form P(n,N)+aQ1(n,a,b,N,A,B)+bQ2(n,a,b,N,A,B)
#with degree of P d1, and degrees of Q1 and Q2 d2
#annihilating F(n,a,b):= The number of n-step walks 
#from the origin to [a,b] staying in
#the quarter-plane using the steps in the set of Steps Steps
#For example, try:GuessOper1(n,a,b,N,A,B,1,1,{[0,1],[1,0],[0,-1],[-1,0]});
GuessOper1:=proc(n,a,b,N,A,B,d1,d2,Steps) 
local ope,var,n1,a1,b1,eq,eq1,K,j,mu,dega,degb,degn,coe,var1:

ope:=GenOper(n,a,b,N,A,B,d1,d2):
var:=ope[2]:
ope:=ope[1]:
eq:={}:

for K from 0 while nops(eq)<nops(var)+10 do
 for n1 from 0 to K do
 for a1 from 0 to K-n1 do
 for b1 from 0 to K-n1-a1 do
  eq1:=0:
   for j from 1 to nops(ope) do
     mu:=op(j,ope):
      dega:=degree(mu,A):
      degb:=degree(mu,B):
      degn:=degree(mu,N):
      coe:=coeff(coeff(coeff(mu,A,dega),B,degb),N,degn):
       eq1:=eq1+subs({n=n1,a=a1,b=b1},coe)*W(n1+degn,a1+dega,b1+degb,Steps):
       eq:=eq union {eq1}:
  od:
od:
od:
od:
od:

var1:=solve(eq,var):
subs(var1,ope):

end:



#GenPureOper(n,a,b,N,d1): a recurrence operator 
#of the form P(n,a,b,N)
#of  degree d1
#For example, try:GenPureOper(n,a,b,N,d1);
GenPureOper:=proc(n,a,b,N,d1) local ope,var,i1,i2,i3,i4,c1:
ope:=0:
var:={}:
for i1 from 0 to d1 do
for i2 from 0 to d1-i1 do
for i3 from 0 to d1-i1-i2 do
for i4 from 0 to d1-i1-i2-i3 do
 ope:=ope+c1[i1,i2,i3,i4]*n^i1*N^i2*a^i3*b^i4:
 var:=var union {c1[i1,i2,i3,i4]}:
od:
od:
od:
od:
ope,var:
end:


#GuessPureOper1(n,a,b,N,d1,Steps): 
#guesses a a recurrence operator 
#of the form P(n,a,b,N)
#of total degree d
#annihilating F(n,a,b):= The number of n-step walks 
#from the origin to [a,b] staying in
#the quarter-plane using the steps in the set of Steps Steps
#For example, try:
#GuessPureOper1(n,a,b,N,2,{[0,1],[1,0],[0,-1],[-1,0]});
GuessPureOper1:=proc(n,a,b,N,d1,Steps) 
local ope,var,n1,a1,b1,eq,eq1,K,j,coe,var1,ind,i,gu1,gu:

ope:=GenPureOper(n,a,b,N,d1):
var:=ope[2]:
ope:=ope[1]:
eq:={}:

for K from 0 while nops(eq)<nops(var)+60 do
 for n1 from 0 to 4*K do
 for a1 from 0 to K do
 for b1 from 0 to K do
  eq1:=0:
   for j from 0 to degree(ope,N) do
      coe:=coeff(ope,N,j):
       eq1:=eq1+subs({n=n1,a=a1,b=b1},coe)*W(n1+j,a1,b1,Steps):
  od:
       eq:=eq union {eq1}:
od:
od:
od:
od:
eq:=eq minus {0}:

var1:=solve(eq,var):
ind:={}:

for j from 1 to nops(var1) do
 if op(1,op(j,var1))=op(2,op(j,var1)) then
   ind:=ind union {op(2,op(j,var1))}:
 fi:
od:


ope:=subs(var1,ope):

gu1:={seq(coeff(ope,i,1), i in ind)}:

gu:={}:

for i from 1 to nops(gu1) do
 if diff(gu1[i],N)<>0 then
  gu:=gu union {gu1[i]}:
 fi:
od:
gu:
end:

#GuessRec0(Steps,a,b,n,N,MaxC): guesses a recurrence operator
#ope(n,N), of complexity MaxC
#for the sequence: number of ways of going from the origin
#to [a,b] using data for n<=L, for example, try:
#GuessRec0({[1,1],[-1,0],[0,-1]},0,0,n,N,40);
GuessRec0:=proc(Steps,a,b,n,N,MaxC) local  gu,i,n1,a1,b1,r,j,L,ope:
L:=(2+MaxC)^2+4:

gu:=[seq(W(n1,a,b,Steps),n1=0..L)]:

for i from 1 while gu[i]=0 do od:

a1:=i:
gu:=[op(a1..nops(gu),gu)]:

for i from 2 to nops(gu) while gu[i]=0 do od:
b1:=i-1:


if {seq(seq(gu[j+b1*r],r=0..trunc((nops(gu)-j)/b1)), j=2..b1)}<>{0} then
RETURN(FAIL):
fi:


L:=a1+b1*L:
gu:=[seq(W(n1,a,b,Steps),n1=0..L)]:


gu:=[seq(gu[a1+b1*r],r=0..trunc((nops(gu)-a1-1)/b1))]:



ope:=Findrec(gu,n,N,MaxC):

if ope=FAIL then
 RETURN(FAIL):
else
 RETURN(ope,a1-1,b1):
fi:

end: