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


print(`First Written: Nov., 2007: tested for Maple 10 `):
print(`Version of Nov. 1, 2007:  `):
print():
print(`This is WalkCarefully, A Maple package`):
print(`To study walks on the positive orthant`):
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/WalkCarefully .`):
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`):
print(`  Comp, Findrec, GoodWalks`):
else
 ezra(args):
fi:
end:

ezra:=proc()
if args=NULL then

print(` WalkCarefully: A Maple package for  `):

print(`The MAIN procedures are`):
print(`NuGoodWalks , SeqDiag, SeqGenDiag, SeqTotGoodWalks, TotGoodWalks `):
print(`Zinn `):



elif nargs=1 and args[1]=Comp then
print(`Comp(n,k): all the vectors of length k, with entries`):
print(`that non-negative integers, that add up to n`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to `):
print(`find a linear recurrence equation with`):
print(`poly coffs.`):
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]=GoodWalks then
print(`GoodWalks(a,Restrictions): The set of walks from`):
print(`the origin to a (where a is a list of non-negative`):
print(`integers, in the k(:=nops(a))-dimensional positive  `):
print(`(cubic) lattice with unit positive steps and`):
print(`staying in the region add(b[i]*x[i],i=1..k)+b[k+1]>=0`):
print(`for each b in Restrictions. For example for`):
print(`the (Bousquet-Melou)-Petkovsek walk, try`):
print(`GoodWalks([3,4],{[-1,2,1],[2,-1,1]});`):
print(`For the Kreweras walks try:`):
print(`GoodWalks([3,4,5],{[-1,0,1,0],[0,-1,1,0]});`):
print(`For the Gessel walks try:`):
print(`GoodWalks([3,4,5,2],{[1,-1,1,-1,0],[0,0,1,-1,0]});`):


elif nargs=1 and args[1]=NuGoodWalks then
print(`NuGoodWalks(a,Restrictions): The number of walks from`):
print(`the origin to a (where a is a list of non-negative`):
print(`integers, in the k(:=nops(a))-dimensional positive  `):
print(`(cubic) lattice with unit positive steps and`):
print(`staying in the region add(b[i]*x[i],i=1..k)+b[k+1]>=0`):
print(`for each b in Restrictions. For example for`):
print(`the (Bousquet-Melou)-Petkovsek walk, try`):
print(`NuGoodWalks([3,4],{[-1,2,1],[2,-1,1]});`):
print(`For the Kreweras walks try:`):
print(`NuGoodWalks([3,4,5],{[-1,0,1,0],[0,-1,1,0]});`):
print(`For the Gessel walks try:`):
print(`NuGoodWalks([3,4,5,2],{[1,-1,1,-1,0],[0,0,1,-1,0]});`):

elif nargs=1 and args[1]=SeqDiag then
print(`SeqDiag(N,k,Restrictions): The first N terms of the sequence`):
print(`number of good lattice-walks in the k-dim positive orthant`):
print(`with unit positive steps staying within the restrictions Restrictions`):
print(`to the point [i,..,i], for i=1..N`):
print(`For example, for 2D ballot walks to the diagonal, try:`):
print(`SeqDiag(10,2,{[1,-1,0]});`):
print(``):
print(`For the (Bousquet-Melou)-Petkovsek walk, try`):
print(`SeqDiag(30,2,{[-1,2,1],[2,-1,1]});`):
print(`For the Kreweras walks try:`):
print(`SeqDiag(30,3,{[-1,0,1,0],[0,-1,1,0]});`):
print(`For the Gessel walks try:`):
print(`SeqDiag(30,4,{[1,-1,1,-1,0],[0,0,1,-1,0]});`):

elif nargs=1 and args[1]=SeqGenDiag then
print(`SeqGenDiag(sta,N,Restrictions): The first N terms of the sequence`):
print(`number of good lattice-walks in the k(=nops(sta))-`):
print(`dim positive orthant`):
print(`with unit positive steps staying within the restrictions Restrictions`):
print(`to the point sta+[i$k], for i=1..N`):
print(`For example, for 2D ballot walks to the diagonal, try:`):
print(`SeqGenDiag([0,0],10,{[1,-1,0]});`):
print(``):
print(`For the (Bousquet-Melou)-Petkovsek walk, try`):
print(`SeqGenDiag([0,0],30,{[-1,2,1],[2,-1,1]});`):
print(`For the Kreweras walks try:`):
print(`SeqGenDiag([0,0,0],30,{[-1,0,1,0],[0,-1,1,0]});`):
print(`For the Gessel walks try:`):
print(`SeqGenDiag([0,0,0,0],30,{[1,-1,1,-1,0],[0,0,1,-1,0]});`):


elif nargs=1 and args[1]=SeqTotGoodWalks then
print(`SeqTotGoodWalks(N,k,Restrictions): The first N steps`):
print(`in the sequence enumerating the number of n-step walks`):
print(`in the k-dim positive orthant with unit steps under`):
print(`the retrictions of Restrictions. For example do:`):
print(`SeqTotGoodWalks(10,2,{[1,-1,0]});`):


elif nargs=1 and args[1]=TotGoodWalks then
print(`TotGoodWalks(n,k,Restrictions): The number of walks from`):
print(`the origin with n steps`):
print(` in the k-dimensional positive  `):
print(`(cubic) lattice with unit positive steps and`):
print(`staying in the region add(b[i]*x[i],i=1..k)+b[k+1]>=0`):
print(`for each b in Restrictions. For example for`):
print(`the (Bousquet-Melou)-Petkovsek walk, try`):
print(`TotGoodWalks(7,2,{[-1,2,1],[2,-1,1]});`):
print(`For the Kreweras walks try:`):
print(`TotGoodWalks(12,3,{[-1,0,1,0],[0,-1,1,0]});`):
print(`For the Gessel walks try:`):
print(`TotGoodWalks(14,4,{[1,-1,1,-1,0],[0,0,1,-1,0]});`):

elif nargs=1 and args[1]=Zinn then
print(`Zinn(L): given a list L, estimates mu and theta`):
print(`such that L[n] is asympt. to mu^n*n^theta `):




else
 print(`There is no such thing as`, args):

fi:

end:


#GoodWalks(a,Restrictions): The set of walks from
#the origin to a (where a is a list of non-negative
#integers, in the k(:=nops(a))-dimensional positive  
#(cubic) lattice with unit positive steps and
#staying in the region add(b[i]*x[i],i=1..k)+b[k+1]>=0
#for each b in Restrictions. For example for
#the (Bousquet-Melou)-Petkovsek walk, try
#GoodWalks([3,4],{[-1,2,1],[2,-1,1]});
#For the Kreweras walks try:
#GoodWalks([3,4,5],{[-1,0,1,0],[0,-1,1,0]});
#For the Gessel walks try:
#GoodWalks([3,4,5,2],{[1,-1,1,-1,0],[0,0,1,-1,0]});

GoodWalks:=proc(a,Restrictions) local k,i,r,gu,gu1,b,g:
option remember:

if  Restrictions<>{} and {seq(nops(r), r in Restrictions)}<>{nops(a)+1} then
 ERROR(`Bad input`):
fi:

k:=nops(a):

if a=[0$k] then
 RETURN({[]}):
fi:

if min(op(a))<0 then
  RETURN({}):
fi:

if Restrictions<>{} and not 
{seq(evalb(add(r[i]*a[i],i=1..k)+r[k+1]>=0), r in Restrictions)}={true}
                                                            then 
 RETURN({}):
fi:



gu:={}:

for i from 1 to k do
b:=[op(1..i-1,a),a[i]-1,op(i+1..k,a)]:

gu1:=GoodWalks(b,Restrictions):

gu:=gu union {seq([op(g),i],g in gu1)}:
od:
gu:

end:


#NuGoodWalks(a,Restrictions): The number of walks from
#the origin to a (where a is a list of non-negative
#integers, in the k(:=nops(a))-dimensional positive  
#(cubic) lattice with unit positive steps and
#staying in the region add(b[i]*x[i],i=1..k)+b[k+1]>=0
#for each b in Restrictions. For example for
#the (Bousquet-Melou)-Petkovsek walk, try
#NuGoodWalks([3,4],{[-1,2,1],[2,-1,1]});
#For the Kreweras walks try:
#NuGoodWalks([3,4,5],{[-1,0,1,0],[0,-1,1,0]});
#For the Gessel walks try:
#NuGoodWalks([3,4,5,2],{[1,-1,1,-1,0],[0,0,1,-1,0]});

NuGoodWalks:=proc(a,Restrictions) local k,i,r:
option remember:

if  Restrictions<>{} and {seq(nops(r), r in Restrictions)}<>{nops(a)+1} then
 ERROR(`Bad input`):
fi:

k:=nops(a):

if a=[0$k] then
 RETURN(1):
fi:

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

if Restrictions<>{} and not 
{seq(evalb(add(r[i]*a[i],i=1..k)+r[k+1]>=0), r in Restrictions)}={true}
                                                            then 
 RETURN(0):
fi:



add(NuGoodWalks([op(1..i-1,a),a[i]-1,op(i+1..k,a)],Restrictions),i=1..k):

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



#Comp(n,k): all the vectors of length k, with entries
#that non-negative integers, that add up to n
Comp:=proc(n,k) local gu,gu1,i,g:
option remember:

if k=0 then
 if n=0 then
   RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:

gu:={}:
for i from 0 to n do
 gu1:=Comp(n-i,k-1):
  gu:=gu union {seq([op(g),i],g in gu1)}:
od:
gu:
end:

#TotGoodWalks(n,k,Restrictions): total number of n-step walks
#in the k-dim positive orthant of the cubic lattice
#with unit steps obeying the conditions specificed by Restictions
##for each b in Restrictions. For example for
#the (Bousquet-Melou)-Petkovsek walk, try
#TotGoodWalks(7,2,{[-1,2,1],[2,-1,1]});
#For the Kreweras walks try:
#TotGoodWalks(7,3,{[-1,0,1,0],[0,-1,1,0]});
#For the Gessel walks try:
#TotGoodWalks(7,4,{[1,-1,1,-1,0],[0,0,1,-1,0]});

TotGoodWalks:=proc(n,k,Restrictions) local gu,g,r:
option remember:

if  Restrictions<>{} and {seq(nops(r), r in Restrictions)}<>{k+1} then
 ERROR(`Bad input`):
fi:

gu:=Comp(n,k):

add(NuGoodWalks(g,Restrictions), g in gu):
end:


#SeqTotGoodWalks(N,k,Restrictions): The first N steps
#in the sequence enumerating the number of n-step walks
#in the k-dim positive orthant with unit steps under
#the retrictions of Restrictions. For example do:
#SeqTotGoodWalks(10,2,{[1,-1,0]});

SeqTotGoodWalks:=proc(N,k,Restrictions) local n:
[seq(TotGoodWalks(n,k,Restrictions),n=1..N)]:
end:


sn:=proc(resh,n):
-1/log(op(n+1,resh)*op(n-1,resh)/op(n,resh)^2):
end:

Zinn:=proc(resh)
local s1,s2,n:
n:=nops(resh)-2:
s1:=sn(resh,n):
s2:=sn(resh,n-1):
evalf(2*(s1+s2)/(s1-s2)^2),
evalf(sqrt(op(n+1,resh)/op(n-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
end:
 
 
 

#SeqDiag(N,k,Restrictions): The first N terms of the sequence
#number of good lattice-walks in the k-dim positive orthant
#with unit positive steps staying within the restrictions Restrictions
#to the point [i,..,i] for i=1..N
#For example, try:
#SeqDiag(10,2,{[1,-1,0]});
SeqDiag:=proc(N,k,Restrictions) local i:
if Restrictions<>{} and 
 {seq(nops(Restrictions[i]),i=1..nops(Restrictions))}<>{k+1} then
   ERROR(`bad input`):
fi:

[seq(NuGoodWalks([i$k],Restrictions),i=1..N)]:
end:


#SeqGenDiag(sta,N,k,Restrictions): The first N terms of the sequence
#number of good lattice-walks in the k-dim positive orthant
#with unit positive steps staying within the restrictions Restrictions
#to the point sta+[i,..,i] for i=1..N
#For example, try:
#SeqGenDiag([0,0],10,2,{[1,-1,0]});
SeqGenDiag:=proc(sta,N,Restrictions) local i,j,k:
k:=nops(sta):
if Restrictions<>{} and 
 {seq(nops(Restrictions[i]),i=1..nops(Restrictions))}<>{k+1} then
   ERROR(`bad input`):
fi:

[seq(NuGoodWalks([seq(sta[j]+i,j=1..k)],Restrictions),i=1..N)]:
end: