######################################################################
##WALKSab: Save this file as  WALKSab                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read WALKSab<Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: June 3, 2009
 
print(`Created:  June 3, 2009`):
print(` This is WALKSab `):
print(`That automatically generates, and solves`):
print(`the algebraic equations for the generating functions`):
print(`let's call them F[0,i](x[a],x[-b]), `):
print(`for i=0,1, ..., b-1 `):
print(`for walks in the upper plane starting at the origin`):
print(`with steps (1,a) and (1,-b), such that their graphs`):
print(`are always weakly above the origin and end-up`):
print(`at a point where y=i (and x anything)`):
print(``):
print(`It is one of the packages that accompany the article `):
print(` "An Experimental Mathematics Perspective of the Old, and Still open Question, of when to Stop" `):
print(`by Luis Medina and and Doron Zeilberger,`):
print(`available from both authors' websites.`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):
 print(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(` CheckGFs, Fab, F1ab,OneStep, Tay `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Eqs, GFs, GFsAll, GFspq, Tay1,Tay1d `):



elif nops([args])=1 and op(1,[args])=CheckGFs then
print(`CheckGFs(a,b,N): checks GFs against Tay `):
print(`For example, try:`):
print(`CheckGFs(1,1,20);`):
print(`CheckGFs(3,2,20);`):


elif nops([args])=1 and op(1,[args])=Eqs then
print(`Eqs(x,a,b,F): The set of equations expressed in terms`):
print(`of pairs [LHS, RHS] for walks in the first quadrant`):
print(`with steps [1,a], [1,-b] always staying weakly above the x-axis`):
print(`F[i,j](resp. G[i,j]) denote walks that start at [0,i] `):
print(`and end at the line y=i that are weakly (resp. strictly) `):
print(`above the x-axis.`):
print(`For example, try:`):
print(`Eqs(x,1,1,F);`):
print(`Eqs(x,2,3,F);`):


elif nops([args])=1 and op(1,[args])=Fab then
print(`Fab(a,b,n,i): the number of walks with steps (1,a) (1,-b) from`):
print(`the origin to (n,i) staying weakly above the origin.`):
print(`For example, try:`):
print(`Fab(1,1,10,0);`):

elif nops([args])=1 and op(1,[args])=F1ab then
print(`F1ab(a,b,n): the number of walks with steps (1,a) (1,-b) from`):
print(`the origin with n steps staying weakly above the origin.`):
print(`and at the last step, going down`):
print(`For example, try:`):
print(`F1ab(1,1,10);`):

elif nops([args])=1 and op(1,[args])=GFs then
print(`GFs(x,a,b): The generating functions `):
print(`F[0,j](j=0..b-1) for walks `):
print(`with steps [1,a], [1,-b] always staying weakly above the x-axis`):
print(`and ending at a point where y=j.`):
print(`It also returns, as the second component, the full generating`):
print(`function in x[a], x[-b] for all `):
print(` walks that are always weakly above the x-axis`):
print(` but at the  very last step go below the x-axis and stop.`):
print(`For example, try:`):
print(`GFs(x,1,1);`):
print(`GFs(x,2,3);`):

elif nops([args])=1 and op(1,[args])=GFsAll then
print(`GFsAll(x,a,b): All The generating functions `):
print(`F[i,j] for walks `):
print(`with steps [1,a], [1,-b] always staying weakly above the x-axis`):
print(`and ending at a point where y=j.`):
print(`It also returns, as the second component, the full generating`):
print(`function in x[a], x[-b] for all `):
print(` walks that are always weakly above the x-axis`):
print(` but at the  very last step go below the x-axis and stop.`):
print(`For example, try:`):
print(`GFsAll(x,1,1);`):
print(`GFsAll(x,2,3);`):


elif nops([args])=1 and op(1,[args])=GFspq then
print(`GFspq(a,b,p,q): The probability of eventually exiting the`):
print(`positive real axis if you start at the origin`):
print(`with steps a (prob. p) and -b (prob. q). For example, try: `):
print(`GFspq(1,1,1/2,1/2);`):
print(`GFspq(2,3,1/2,1/2);`):

elif nops([args])=1 and op(1,[args])=OneStep then
print(`OneStep(S,t,n): given a triple S=[Trans,Var,Current]`):
print(`a variable t, and a current subsitution with Taylor`):
print(`series up to t^n, upgrades it to t^(n+1) . `):
print(`For example, try:`):
print(`OneStep(Eqst(x,1,1,F,t),t,1);`):

elif nops([args])=1 and op(1,[args])=Tay then
print(`Tay(x,a,b,t,N): The Taylor expansion up to order N`):
print(`of the GFs(x,a,b)`):
print(`For example, try:`):
print(`Tay(x,1,1,t,10);`):

elif nops([args])=1 and op(1,[args])=Tay1 then
print(`Tay1(a,b,t,N): The Taylor expansion up to order N`):
print(`of the generating function for walks`):
print(`that are above the x-axis until the last step, and then `):
print(`go below and quit, with steps (1,a) and (1,-b). `):
print(`This is done via the algebraic system`):
print(`For example, try:`):
print(`Tay1(1,1,t,10);`):

elif nops([args])=1 and op(1,[args])=Tay1d then
print(`Tay1d(a,b,t,N): The Taylor expansion up to order N`):
print(`of the generating function for walks`):
print(`that are above the x-axis until the last step, and then `):
print(`go below and quit, with steps (1,a) and (1,-b). `):
print(`This is done directly`):
print(`For example, try:`):
print(`Tay1d(1,1,t,10);`):

else
print(`There is no ezra for`,args):
fi:
 
end:




#EqsOld(x,a,b,F,G): The set of equations expressed in terms
#of pairs [LHS, RHS] for walks in the first quadrant
#with steps [1,a], [1,-b] always staying weakly above the x-axis
#F[i,j](resp. G[i,j]) denote walks that start at [0,i] 
#and end at the line y=i that are weakly (resp. strictly) above the x-axis.
#For example, try:
#Eqs(x,1,1,F,G);
#Eqs(x,2,3,F,G);
EqsOld:=proc(x,a,b,F,G) local eq,AlreadyDone,StillToDo,pt,i,j:
eq:={}:
StillToDo:={F[0,0]}:
AlreadyDone:={}:

while StillToDo<>{} do

pt:=StillToDo[1]:
i:=op(1,pt): j:=op(2,pt):

if op(0,pt)=F then
 if i=0 and j=0 then
  AlreadyDone:=  AlreadyDone union {F[0,0]}:
 eq:=eq union {[F[0,0],1+G[0,0]*F[0,0]]}:
 StillToDo:=StillToDo union {G[0,0],F[0,0]} minus AlreadyDone:
 elif i>0 and j=0 then
  AlreadyDone:=  AlreadyDone union {F[i,0]}:
  eq:=eq union {[F[i,0],G[i,0]*F[0,0]]}:
 StillToDo:=StillToDo union {G[i,0],F[0,0]} minus AlreadyDone:
 elif i=0 and j>0 then
  AlreadyDone:=  AlreadyDone union {F[0,j]}:
  eq:=eq union {[F[0,j],F[0,0]*G[0,j]]}:
 StillToDo:=StillToDo union {F[0,0],G[0,j]} minus AlreadyDone:
 elif i>0 and j>0 then
  AlreadyDone:=  AlreadyDone union {F[i,j]}:
  eq:=eq union {[F[i,j],G[i,j]+G[i,0]*F[0,0]*G[0,j]]}:
 StillToDo:=StillToDo union {G[i,j],G[i,0],F[0,0],G[0,j]} minus AlreadyDone:
 else
  ERROR(`Something is wrong`):
 fi:

elif op(0,pt)=G then

 if i=0 and j=0 then
  AlreadyDone:=  AlreadyDone union {G[0,0]}:
 eq:=eq union {[G[0,0],x[a]*G[a,b]*x[-b]]}:
 StillToDo:=StillToDo union {G[a,b]} minus AlreadyDone:
 elif i>0 and j=0 then
  AlreadyDone:=  AlreadyDone union {G[i,0]}:
  eq:=eq union {[G[i,0],G[i,b]*x[-b] ]}:
 StillToDo:=StillToDo union {G[i,b]} minus AlreadyDone:
 elif i=0 and j>0 then
  AlreadyDone:=  AlreadyDone union {G[0,j]}:
  eq:=eq union {[G[0,j],x[a]*G[a,j]]}:
 StillToDo:=StillToDo union {G[a,j]} minus AlreadyDone:
 elif i>0 and j>0 then
  AlreadyDone:=  AlreadyDone union {G[i,j]}:
  eq:=eq union {[G[i,j],F[i-1,j-1]]}:
 StillToDo:=StillToDo union {F[i-1,j-1]} minus AlreadyDone:
 else
  ERROR(`Something is wrong`):
 fi:
else
 ERROR(`Something is wrong`):
fi:

od:

eq,AlreadyDone:

end:





#GFs(x,a,b): The generatinf functions 
#F[0,j](j=0..b-1) for walks 
#with steps [1,a], [1,-b] always staying weakly above the x-axis
#and ending at a point where y=j
#For example, try:
#GFs(x,1,1);
#GFsOld(x,2,3);
GFsOld:=proc(x,a,b) local F,G,gu,eq,var,i,j:
gu:=EqsOld(x,a,b,F,G):
var:=gu[2]:
eq:=gu[1]:
eq:={seq(eq[i][1]-eq[i][2],i=1..nops(eq))}:
var:=solve(eq,var):
[seq(subs(var,F[0,j]),j=0..b-1)]:
end:


#Eqs(x,a,b,F): The set of equations expressed in terms
#of pairs [LHS, RHS] for walks in the first quadrant
#with steps [1,a], [1,-b] always staying weakly above the x-axis
#F[i,j] denote walks that start at [0,i] 
#and end at the line y=i that are weakly (resp. strictly) above the x-axis.
#it also returns the initial conditions
#For example, try:
#Eqs(x,1,1,F);
#Eqs(x,2,3,F);
Eqs:=proc(x,a,b,F) local eq,AlreadyDone,StillToDo,pt,i,j,Ini,a1,i1,j1:
eq:={}:
StillToDo:={F[0,0]}:
AlreadyDone:={}:

while StillToDo<>{} do

pt:=StillToDo[1]:
i:=op(1,pt): j:=op(2,pt):

 if i=0 and j=0 then
  AlreadyDone:=  AlreadyDone union {F[0,0]}:
 eq:=eq union {[F[0,0],1+x[a]*F[a-1,b-1]*x[-b]*F[0,0]]}:
 StillToDo:=StillToDo union {F[a-1,b-1],F[0,0]} minus AlreadyDone:
 elif i>0 and j=0 then
  AlreadyDone:=  AlreadyDone union {F[i,0]}:
  eq:=eq union {[F[i,0],F[i-1,b-1]*x[-b]*F[0,0]]}:
 StillToDo:=StillToDo union {F[i-1,b-1],F[0,0]} minus AlreadyDone:
 elif i=0 and j>0 then
  AlreadyDone:=  AlreadyDone union {F[0,j]}:
  eq:=eq union {[F[0,j],F[0,0]*x[a]*F[a-1,j-1]]}:
 StillToDo:=StillToDo union {F[0,0],F[a-1,j-1]} minus AlreadyDone:
 elif i>0 and j>0 then
  AlreadyDone:=  AlreadyDone union {F[i,j]}:
  eq:=eq union {[F[i,j],F[i-1,j-1]+F[i-1,b-1]*x[-b]*F[0,0]*x[a]*F[a-1,j-1]]}:
 StillToDo:=StillToDo union 
{F[i-1,j-1],F[i-1,b-1],F[0,0],F[a-1,j-1]} minus AlreadyDone:
 else
  ERROR(`Something is wrong`):
 fi:



od:

AlreadyDone:=convert(AlreadyDone,list):

Ini:={}:

for a1 in AlreadyDone do
i1:=op(1,a1): j1:=op(2,a1):
 if i1=j1 then
  Ini:=Ini union {a1=1}:
 else
  Ini:=Ini union {a1=0}:
 fi:
od:


eq,AlreadyDone,Ini:


end:


#GFs(x,a,b): The generating functions 
#F[0,j](j=0..b-1) for walks 
#with steps [1,a], [1,-b] always staying weakly above the x-axis
#and ending at a point where y=j
#For example, try:
#GFs(x,1,1);
#GFs(x,2,3);
GFs:=proc(x,a,b) local F,gu,eq,var,i,j:
gu:=Eqs(x,a,b,F):
var:=convert(gu[2],set):
eq:=gu[1]:
eq:={seq(eq[i][1]-eq[i][2],i=1..nops(eq))}:
var:=solve(eq,var):
[seq(subs(var,F[0,j]),j=0..b-1)],
normal(x[-b]*add(subs(var,F[0,j]),j=0..b-1)):
end:


#Fab(a,b,n,i): the number of walks with steps (1,a) (1,-b) from
#the origin to (n,i) staying weakly above the origin
#For example, try:
#Fab(1,1,10,0);
Fab:=proc(a,b,n,i) 
option remember:

if n<0 or i<0 then
 RETURN(0):
fi:
if n=0 then
  if i=0 then
    RETURN(1):
   else
    RETURN(0):
  fi:
fi:

Fab(a,b,n-1,i-a)+Fab(a,b,n-1,i+b):

end:








#GFspq(a,b,p,q): The probability of eventually exiting the
#positive real axis if you start at the origin
#with steps a (prob. p) and -b (prob. q)
#GFspq(1,1,1/2,1/2);
#GFspq(2,3,1/2,1/2);
GFspq:=proc(a,b,p,q) local F,gu,eq,var,i,j,x,v:
gu:=Eqs(x,a,b,F):
var:=gu[2]:
eq:=gu[1]:
eq:=subs({x[a]=p,x[-b]=q},eq):
eq:={seq(eq[i][1]-eq[i][2],i=1..nops(eq))}:
var:={solve(eq,var)}:

{seq(q*add(subs(v,F[0,j]),j=0..b-1),v in var)}:

end:

#GFsAll(x,a,b): The generating functions 
#F[0,j](j=0..b-1) for walks 
#with steps [1,a], [1,-b] always staying weakly above the x-axis
#and ending at a point where y=j
#For example, try:
#GFsAll(x,1,1);
#GFsAll(x,2,3);
GFsAll:=proc(x,a,b) local F,gu,eq,var,i,v,var1:
gu:=Eqs(x,a,b,F):
var:=gu[2]:
eq:=gu[1]:
eq:={seq(eq[i][1]-eq[i][2],i=1..nops(eq))}:
var1:=solve(eq,var):
{seq([v,subs(var1,v)],v in var)}:
end:




#Eqst(x,a,b,F,t): The set of equations expressed in terms
#of pairs [LHS, RHS] for walks in the first quadrant
#with steps [1,a], [1,-b] always staying weakly above the x-axis
#F[i,j] denote walks that start at [0,i] 
#and end at the line y=i that are weakly (resp. strictly) above the x-axis.
#it also returns the initial conditions
#For example, try:
#Eqst(x,1,1,F,t);
#Eqst(x,2,3,F,t);
Eqst:=proc(x,a,b,F,t) local gu:
gu:=Eqs(x,a,b,F):

[subs({x[a]=t*x[a],x[-b]=t*x[-b]},gu[1]),gu[2],gu[3]]:
end:



#OneStep(S,t,n): given a triple S=[Trans,Var,Current]
#a variable t, and a current subsitution with Taylor
#series up to t^n, upgrades it to t^(n+1)
#For example, try:
#OneStep(Eqst(x,1,1,F,t),t,1);
OneStep:=proc(S,t,n) local gu1,gu2,gu3,g1,kha,i,kha1:
gu1:=S[1]:
gu2:=S[2]:
gu3:=S[3]:

kha:={}:

for g1 in gu1 do
kha1:=expand(subs(gu3,g1[2])):
kha1:=add(expand(coeff(kha1,t,i))*t^i,i=0..n):
kha:=kha union {g1[1]=kha1}:
od:

[gu1,gu2,kha]:
end:




#Tay(x,a,b,t,N): The Taylor expansion up to order N
#of the GFs(x,a,b)
#For example, try:
#Tay(x,1,1,t,10);
Tay:=proc(x,a,b,t,N) local S,F,n,j:
S:=Eqst(x,a,b,F,t):

for n from 1 to N do
 S:=OneStep(S,t,n) :
od:


[[seq(expand(subs(S[3],F[0,j])),j=0..b-1)],
expand(t*x[-b]*add(subs(S[3],F[0,j]),j=0..b-1))]:
end:



#CheckGFs(a,b,N): cheks GFs
CheckGFs:=proc(a,b,N) local gu,i1,j,x,t,mu:
gu:=GFs(x,a,b) :
gu:=subs({x[a]=t*x[a],x[-b]=t*x[-b]},gu[1]),
subs({x[a]=t*x[a],x[-b]=t*x[-b]},gu[2]):

gu:=[[seq(taylor(gu[1][j],t=0,N+2),j=1..nops(gu[1]))],taylor(gu[2],t=0,N+2)]:

gu:=[[seq(add(coeff(gu[1][j],t,i1)*t^i1,i1=0..N),j=1..nops(gu[1]) )],
add(coeff(gu[2],t,i1)*t^i1,i1=0..N)]:

mu:=Tay(x,a,b,t,N+2):

mu:=[[seq(add(coeff(mu[1][j],t,i1)*t^i1,i1=0..N),j=1..nops(mu[1]) )],
add(coeff(mu[2],t,i1)*t^i1,i1=0..N)]:

expand([seq(gu[1][j]-mu[1][j],j=1..nops(mu[1]))]),expand(gu[2]-mu[2]):

end:



#Tay1(a,b,t,N): The Taylor expansion up to order N
#of the GFs(x,a,b)[2]
#For example, try:
#Tay1(1,1,t,10);
Tay1:=proc(a,b,t,N) local S,F,n,j,gu,x,i1:
S:=subs({x[a]=1,x[-b]=1},Eqst(x,a,b,F,t)):

for n from 1 to N do
 S:=OneStep(S,t,n) :
od:

gu:=expand(t*add(subs(S[3],F[0,j]),j=0..b-1)):

add(coeff(gu,t,i1)*t^i1,i1=1..N):
end:


#F1ab(a,b,n): the number of walks with steps (1,a) (1,-b) from
#the origin to (n,i) staying weakly above the origin and then
#going down
#For example, try:
#F1ab(1,1,10);
F1ab:=proc(a,b,n) local i:
add(Fab(a,b,n-1,i),i=0..b-1) :
end:


#Tay1d(a,b,t,N): The Taylor expansion up to order N
#of the GFs(x,a,b)[2], done directly
#For example, try:
#Tay1d(1,1,t,10);
Tay1d:=proc(a,b,t,N) local i:
add(F1ab(a,b,i)*t^i,i=1..N):
end: