######################################################################
##STADJE: Save this file as  STADJE                                  #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read STADJE<Enter>                                 #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#First posted: June 23, 2009
 
print(`First posted:  June 23, 2009`):
print(` This is STADJE `):
print(`It is one of the packages that accompany the article `):
print(`"An Experimental Mathematics Perspective on the Old and still`):
print(`open, question of when to stop?"`):
print(`by Luis A. Medina and and Doron Zeilberger`):
print(`available from Zeilberger's website`):
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: BR,FindR, GettingOut `):
 print(`gExGain, gPrEsc, gW,`):
 print(` NetGain, PrSeqs, Strip, Zinn `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  `):
 print(` ExGain,MiddleAB, LineAB, `):
 print(` PrEsc, PrEscF, PrTouch, PrTouchF, `):
 print(`  SipurLineAB, SipurMiddleAB, SipurPr ,  `):
 print(` W, Ws  `):
 print(` `):



elif nops([args])=1 and op(1,[args])=BR then
print(`BR(R): the boundary of the region R expressed`):
print(`as a set of pairs [ptInRegion,ptNotInRegion]`):
print(`For example, try:`):
print(`BR([[1,3]$3]);`):


elif nops([args])=1 and op(1,[args])=ExGain then
print(`ExGain(a,b,K): the expected gain of keep going`):
print(`if you currently have`):
print(`a heads and b tails right now`):
print(`in less than (a+b)*K steps`):
print(`and you stop as soon as you get a better score`):
print(`or die at age (a+b)*K, where you are allowed to collect`):
print(`1/2 in case your score is less than that.`):
print(`For example, try:`):
print(`ExGain(5,3,100);`):


elif nops([args])=1 and op(1,[args])=FindR then
print(`FindR(a,b,K): the region y>=b/a*x`):
print(`For example, try:`):
print(`FindR(5,3,10);`):



elif nops([args])=1 and op(1,[args])=GettingOut then
print(`GettingOut(a,b): all the pairs [i1,j1] such that [a*i+i1,b*i+j1]->[a*i+i1+1,b*i+j1] is the exiting edge`):
print(`in the gambler's escape from y>=b/a*x `):
print(`For example, try:`):
print(`GettingOut(5,3);`):



elif nops([args])=1 and op(1,[args])=gExGain then
print(`gExGain(a0,b0,R): the expected gain of travelling in region R`):
print(`starting at a0,b0 where the payoff is`):
print(`max(1/2,x/(x+y)) as soon as you get out.`):
print(`For example, try:`):
print(`gExGain(1,1,[[1,5]$5]);`):


elif nops([args])=1 and op(1,[args])=gPrEsc then
print(`gPrEsc(a0,b0,R): the probability of escaping the region R`):
print(`starting at a0,b0`):
print(`For example, try:`):
print(`gPrEsc(1,1,[[1,5]$5]);`):


elif nops([args])=1 and op(1,[args])=gW then
print(`gW(a0,b0,x,y,R): the number of walks with unit steps [1,0],[0,1]`):
print(`starting at a0,b0 and`):
print(`staying in the region R. For example, try:`):
print(`gW(1,1,5,5,[[1,5]$5]);`):


elif nops([args])=1 and op(1,[args])=LineAB then
print(`LineAB(a,b,K,n): the first K terms of the sequence`):
print(`W(a*i,b*i,a,b) (walks from the origin to (a*i,b*i)`):
print(` staying in the region`):
print(`y>=b/a*x), followed by the conjectured asympotics in n.`):
print(`For example, try:`):
print(`LineAB(3,2,100,n);`):

elif nops([args])=1 and op(1,[args])=MiddleAB then
print(`MiddleAB(a,b,K,n): the first K terms of the sequence`):
print(`W(i,i,a,b) (walks from the origin to (i,i) staying in the region`):
print(`y>=b/a*x), followed by the conjectured asympotics in n.`):
print(`For example, try:`):
print(`MiddleAB(3,2,100,n);`):

elif nops([args])=1 and op(1,[args])=NetGain then
print(`NetGain(a,b,K): the relative net gain of keep going`):
print(`if you`):
print(`a heads and b tails right now`):
print(`in less than (a+b)*K steps`):
print(`For example, try:`):
print(`NetGain(5,3,100);`):

elif nops([args])=1 and op(1,[args])=PrEsc then
print(`PrEsc(a,b,K): the probability of escaping the region y>=b/a*x `):
print(`in less than (a+b)*K steps.`):
print(`For example, try:`):
print(`PrEsc(5,3,100);`):

elif nops([args])=1 and op(1,[args])=PrEscF then
print(`PrEscF(a,b,eps): the probability of escaping the region y>=b/a*x `):
print(`in floating point with error eps`):
print(`For example, try:`):
print(`PrEscF(5,3,10^(-11));`):


elif nops([args])=1 and op(1,[args])=PrTouch then
print(`PrTouch(a,b,K): the probability of escaping the region y>=b/a*x `):
print(`or at least touching it, besides the origin`):
print(`in less than (a+b)*K steps.`):
print(`For example, try:`):
print(`PrTouch(5,3,100);`):

elif nops([args])=1 and op(1,[args])=PrTouchF then
print(`PrTouchF(a,b,eps): the probability of escaping the region y>=b/a*x `):
print(`or at least touching it, besides the origin`):
print(`in floating point, with error <eps`):
print(`For example, try:`):
print(`PrTouchF(5,3,10^(-11));`):

elif nops([args])=1 and op(1,[args])=PrSeqs then
print(`PrSeqs(a,b,K): the summand-sequences in the`):
print(`probability of escaping the region y>=b/a*x `):
print(`in less than (a+b)*K steps.`):
print(`For example, try:`):
print(`PrSeqs(5,3,100);`):



elif nops([args])=1 and op(1,[args])=SipurMiddleAB then
print(`SipurMiddleAB(M,K,L,n): the sequences for`):
print(`[W(i,i,a,b),i=1..L] `):
print(`(walks from the origin to (i,i) staying in the region`):
print(`y>=b/a*x), followed by the conjectured asympotics in n`):
print(`for starting at the origin for M>=a>=b>=1 and gcd(a,b)=1`):
print(`it displays the first L terms. `):
print(`SipurMiddleAB(3,100,20,n);`):

elif nops([args])=1 and op(1,[args])=SipurPr then
print(`SipurPr(K,eps): `):
print(`all the probabilities of escaping the region y>=b/a*x `):
print(`and or y>b/a*x  (except at the beginning, of course)`):
print(`with unit steps to the right and up in the discrete lattice`):
print(`starting at the origin for K>=a>=b>=1 and gcd(a,b)=1`):
print(`with error eps`):
print(`For example, try:`):
print(`SipurPr(4,10^(-11));`):

elif nops([args])=1 and op(1,[args])=Strip then
print(`Strip(a,b,alpha,L,K): the strip with parameters alpha and L`):
print(`of length K`):
print(`emanating from [a,b] (a>=b)`):
print(`For example, try:`):
print(`Strip(5,3,1,10,200);`):

elif nops([args])=1 and op(1,[args])=W then
print(`W(x,y,a,b): the number of walks with unit steps [1,0],[0,1]`):
print(`staying in the region y>=b/a*x. For example, try:`):
print(`W(5,10,1,1);`):


elif nops([args])=1 and op(1,[args])=Ws then
print(`Ws(x,y,a,b): the number of walks with unit steps [1,0],[0,1]`):
print(`staying in the region y>b/a*x. For example, try:`):
print(`Ws(5,10,1,1);`):

elif nops([args])=1 and op(1,[args])=Zinn then
 print(` Zinn's method `):
 print(`Zinn(IntegerSequence): Given an increasing sequence a(n) of positive integers , expressed`):
 print(`in terms of a list, estimates the theta and mu such that `):
 print(` a(n) is asympt. to n^(theta)*mu^n. The output is `):
 print(` theta, mu `):

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



#W(x,y,a,b): the number of walks with unit steps [1,0],[0,1]
#staying in the region y>=b/a*x. For example, try:
#W(5,10,1,1);
W:=proc(x,y,a,b) option remember:

if x<0 then
 0:
elif x=0 then
 1:
elif y<b/a*x then
 RETURN(0):
else
RETURN(W(x-1,y,a,b)+W(x,y-1,a,b)):
fi:
end:


#GettingOut(a,b): all the pairs [i1,j1] such that [a*i+i1,b*i+j1]->[a*i+i1+1,b*i+j1] is the exiting edge
#in the gambler's escape from y>=b/a*x 
#For example, try:
#GettingOut(5,3);
GettingOut:=proc(a,b) local j1:
[seq([trunc(a*j1/b),j1],j1=0..b-1)]:
end:



#PrEsc(a,b,K): the probability of escaping the region y>=b/a*x 
#in less than (a+b)*K steps
#For example, try:
#PrEsc(5,3,100);
PrEsc:=proc(a,b,K) local i,gu,g:
gu:=GettingOut(a,b):
add(add(W(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1), g in gu),i=0..K):
end:


#PrEscF1(a,b,K): the probability of escaping the region y>=b/a*x 
#in less than (a+b)*K steps in floating point
#For example, try:
#PrEscF1(5,3,100);
PrEscF1:=proc(a,b,K) local i,gu,g:
gu:=GettingOut(a,b):
add(add(W(a*i+g[1],b*i+g[2],a,b)/2.^((a+b)*i+g[1]+g[2]+1), g in gu),i=0..K):
end:


#ExGain(a,b,K): the expected gain of keep going
#if you
#a heads and b tails right now
#in less than (a+b)*K steps
#For example, try:
#ExGain(5,3,100);
ExGain:=proc(a1,b1,K) local i,gu,g,shever,a,b:
shever:=a1/b1:
a:=numer(shever):
b:=denom(shever):
gu:=GettingOut(a,b):
add(add(
W(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1)*
(a1+a*i+g[1]+1)/(a1+b1+(a+b)*i+g[1]+g[2]+1), g in gu),i=0..K)
+(1-PrEsc(a,b,K))/2:
end:


#ExGainOld(a,b,K): the expected gain of keep going
#if you
#a heads and b tails right now
#in less than (a+b)*K steps
#For example, try:
#ExGainOld(5,3,100);
ExGainOld:=proc(a1,b1,K) local i,gu,g,shever,a,b:
shever:=a1/b1:
a:=numer(shever):
b:=denom(shever):
gu:=GettingOut(a,b):
add(add(
W(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1)*
(a1+a*i+g[1]+1)/(a1+b1+(a+b)*i+g[1]+g[2]+1), g in gu),i=0..K):
end:





#NetGain(a,b,K): the relative net gain of keep going
#if you a heads and b tails right now
#in less than (a+b)*K steps
#For example, try: (if it is neg. stop, if it is positive, stop)
#NetGain(5,3,100);
NetGain:=proc(a,b,K) local gu:
gu:=ExGain(a,b,K): (gu-a/(a+b))/(a/(a+b)):
end:




#PrSeqs(a,b,K): the probability sequences (summands)
#of escaping the region y>=b/a*x 
#in less than (a+b)*K steps
#For example, try:
#PrSeqs(5,3,100);
PrSeqs:=proc(a,b,K) local i,gu,g:

gu:=GettingOut(a,b):

[seq(
[seq(W(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1),i=0..K)],
g in gu)]:
end:





#Strip(a,b,alpha,L,K): the strip with parameters alpha and L
#of length K
#emanating from [a,b] (a>=b)
#For example, try:
#Strip(5,3,1,10,200);
Strip:=proc(a,b,alpha,L,K) local gu,R,R1,i,j:
gu:=[seq(ceil(b/a*i),i=0..a-1)]:

R:=[b$(a-1),seq(seq(b*j+gu[i],i=1..nops(gu)),j=1..trunc(K/nops(gu)))]:
R1:=[seq(trunc(alpha*R[i])+L,i=1..nops(R))]:
[seq([R[i],R1[i]],i=1..nops(R))]:


end:




#gW(a0,b0,x,y,R): the number of walks with unit steps [1,0],[0,1]
#starting at a0,b0 and
#staying in the region R. For example, try:
#gW(1,1,5,5,[[1,5]$5]);
gW:=proc(a0,b0,x,y,R) option remember:

if x<a0 or y<b0 then
 0:
elif x=a0 and y=b0 then
 1:
elif x>nops(R) or y<R[x][1] or y>R[x][2] then
 RETURN(0):
else
RETURN(gW(a0,b0,x-1,y,R)+gW(a0,b0,x,y-1,R)):
fi:
end:




#BR(R): the boundary of the region R expressed
#as a set of pairs [ptInRegion,ptNotInRegion]
#For example, try:
#BR([[1,3]$3]);


BR:=proc(R) local gu,k,j,i:

gu:={}:

for i from 1 to nops(R)-1 do
 if R[i][1]<R[i+1][1] then
  gu:=gu union {[[i,R[i][1]],[i+1,R[i][1]]]}:
 fi:
od:

k:=nops(R):

gu:=gu union {seq([[k,j],[k+1,j]],j=R[k][1]..R[k][2])}:

for i from 1 to nops(R) do
 gu:=gu union {[[i,R[i][2]],[i,R[i][2]+1]]}:
od:

gu:

end:



#gPrEsc(a0,b0,R): the probability of escaping the region R
#starting at a0,b0
#For example, try:
#gPrEsc(1,1,[[1,5]$5]);
gPrEsc:=proc(a0,b0,R) local B,b,gu:
if b0<R[a0][1] or b0>R[a0][2] then
 RETURN(FAIL):
fi:
B:=BR(R):
gu:=0:
for b in B do
 gu:=gu+gW(a0,b0,op(b[1]),R)/2^(b[2][1]+b[2][2]-a0-b0):
od:

gu:

end:


#gExGain(a0,b0,R): the expected gain of travelling in region R
#starting at a0,b0 where the payoff is
#max(1/2,x/(x+y)) as soon as you get out
#For example, try:
#gExGain(1,1,[[1,5]$5]);
gExGain:=proc(a0,b0,R) local B,b,gu:
if b0<R[a0][1] or b0>R[a0][2] then
 RETURN(FAIL):
fi:
B:=BR(R):
gu:=0:
for b in B do
 gu:=gu+gW(a0,b0,op(b[1]),R)/2^(b[2][1]+b[2][2]-a0-b0)*
max(1/2,b[2][1]/(b[2][1]+b[2][2])):
od:

gu:

end:


#FindR(a,b,K): the region y>=b/a*x
#For example, try:
#FindR(5,3,10);
FindR:=proc(a,b,K) local i:

[seq([ceil(b*i/a),K],i=1..K)]:

end:









#PrTouch(a,b,K): the probability of achieveing region y>=b/a*x 
#in less than (a+b)*K steps
#For example, try:
#PrTouch(5,3,100);
PrTouch:=proc(a,b,K) local i,gu,g:
gu:=GettingOut(a,b):
gu:=[op(2..nops(gu),gu)]:
add(add(Ws(a*i+g[1],b*i+g[2],a,b)/2^((a+b)*i+g[1]+g[2]+1), g in gu),i=0..K)
+
add(Ws(a*i-1,b*i,a,b)/2^((a+b)*i),i=1..K)+1/2:
:
end:


#PrTouchF1(a,b,K): the probability of achieveing region y>=b/a*x 
#in less than (a+b)*K steps, in Floating point
#For example, try:
#PrTouchF1(5,3,100);
PrTouchF1:=proc(a,b,K) local i,gu,g:
gu:=GettingOut(a,b):
gu:=[op(2..nops(gu),gu)]:
add(add(Ws(a*i+g[1],b*i+g[2],a,b)/2.^((a+b)*i+g[1]+g[2]+1), g in gu),i=0..K)
+
add(Ws(a*i-1,b*i,a,b)/2.^((a+b)*i),i=1..K)+1/2:
:
end:



#Ws(x,y,a,b): the number of walks with unit steps [1,0],[0,1]
#staying in the region y>b/a*x. For example, try:
#Ws(5,10,1,1);
Ws:=proc(x,y,a,b) option remember:

if x<0 then
 0:
elif x=0 then
 1:
elif y<=b/a*x then
 RETURN(0):
else
RETURN(Ws(x-1,y,a,b)+Ws(x,y-1,a,b)):
fi:
end:


#PrEscF(a,b): the probability of escaping the region y>=b/a*x 
#in floating point with the current accuracy
#For example, try:
#PrEscF(5,3);
PrEscF:=proc(a,b,eps) local K,gu1,gu2:

gu1:=PrEsc(a,b,100):
gu2:=PrEsc(a,b,110):

for K from 2 to 20 while abs(gu1-gu2)>eps do
gu1:=gu2:
gu2:=PrEsc(a,b,100+20*K):
od:
if K=21 then
FAIL:
else
evalf(gu2):
fi:


end:


#PrTouchF(a,b): the probability of passing or touching the line y=b/a*x 
#in floating point with the current accuracy
#For example, try:
#PrTouchF(5,3);
PrTouchF:=proc(a,b,eps) local K,gu1,gu2:

gu1:=PrTouch(a,b,100):
gu2:=PrTouch(a,b,110):

for K from 2 to 20 while abs(gu1-gu2)>eps do
gu1:=gu2:
gu2:=PrTouch(a,b,100+20*K):
od:
if K=21 then
FAIL:
else
evalf(gu2):
fi:


end:


#SipurPr(K,eps): 
#all the probabilities of escaping the region y>=b/a*x 
#and or y>b/a*x  (except at the beginning, of course)
#with unit steps to the right and up in the discrete lattice
#starting at the origin for K>=a>=b>=1 and gcd(a,b)=1
#with error eps
#For example, try:
#SipurPr(4,10^(-11));
SipurPr:=proc(K,eps) local a,b,t0:
t0:=time():
print(`This is the story of all the probabilities of`):
print(`(i) escaping the region y>=b/a*x  `):
print(`(ii) escaping  y>b/a*x  (except at the beginning, of course)`):
print(`with unit steps to the right and up in the discrete lattice`):
print(`each with prob. 1/2 `):
print(`starting at the origin`):
print(`for 2<=b<=a<=`, K):

for a from 2 to K do
for b from 1 to a do
 if gcd(a,b)=1 then
  print(`a=`,a,`b=`,b, `(i)=`, PrEscF(a,b,eps),`(ii)=`, PrTouchF(a,b,eps)):
 fi:
od:
od:
print(`------------------------`):
print(`This took`, time()-t0, `seconds. `):
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:
 



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



#MiddleAB(a,b,K,n): the first K terms of the sequence
#W(i,i,a,b) (walks from the origin to (i,i) staying in the region
#y>=b/a*x), followed by the conjectured asympotics in n
#For example, try:
#MiddleAB(3,2,100,n);
MiddleAB:=proc(a,b,K,n) local gu,i,gu1,err,D1,C:
gu:=[seq(W(i,i,a,b),i=1..K)]:
gu1:=[seq(evalf(gu[i]/4^i*sqrt(i)),i=1..nops(gu))]:
err:=abs(gu1[trunc(K/2)]-gu1[K]):
D1:=trunc(log(1/err)/log(10.))-1:

if D1>2 then
C:=evalf(gu1[nops(gu1)],D1):
RETURN(gu,C*4^n/sqrt(n)):
else
RETURN(gu,FAIL):
fi:


end:


#LineAB(a,b,K,n): the first K terms of the sequence
#W(a*i,b*i,a,b) (walks from the origin to (i,i) staying in the region
#y>=b/a*x), followed by the conjectured asympotics in n
#For example, try:
#LineAB(3,2,100,n);
LineAB:=proc(a,b,K,n) local gu,i,gu1,err,D1,C,c:
gu:=[seq(W(a*i,b*i,a,b),i=1..K)]:
c:=(a+b)^(a+b)/a^a/b^b:
gu1:=[seq(evalf(gu[i]/c^i*i^(3/2)),i=1..nops(gu))]:
err:=abs(gu1[trunc(K/2)]-gu1[K]):
D1:=trunc(log(1/err)/log(10.))-1:



if D1>1 then
C:=evalf(gu1[nops(gu1)],D1):
RETURN(gu,C*c^n/n^(3/2)):
else
RETURN(gu,FAIL):
fi:
end:

#SipurMiddleAB(M,K,L,n): the sequences for
#[W(i,i,a,b),i=1..L] (walks from the origin to (i,i) staying in the region
#y>=b/a*x), followed by the conjectured asympotics in n
#for starting at the origin for M>=a>=b>=1 and gcd(a,b)=1
#it displays the first L terms 
#SipurMiddleAB(3,100,20,n);
SipurMiddleAB:=proc(M,K,L,n) local a,b,t0,gu:
t0:=time():
print(`Below are the sequences`):
print(`enumerating the number of ways of walking from the origin`):
print(`with positive unit steps, in the square lattice`):
print(`to the point (i,i), all the while staying in the`):
print(`the region y>=b/a*x  `):
print(`followed by the conjectured asymptotics (if found) `):
print(`for 2<=b<=a<=`, M):

for a from 2 to M do
for b from 1 to a do
 if gcd(a,b)=1 then
 gu:=MiddleAB(a,b,K,n):
 if gu[2]<>FAIL then
  print(`a=`,a,`b=`,b, `seq.:` [op(1..L,gu[1])]):
  print(`asymp.:`, gu[2]):
else
  print(`a=`,a,`b=`,b, `seq.:` [op(1..L,gu[1])]):
 fi:
fi:
od:
od:
print(`------------------------`):
print(`This took`, time()-t0, `seconds. `):
end:


#SipurLineAB(M,K,L,n): the sequences for
#[W(i,i,a,b),i=1..L] (walks from the origin to (i,i) staying in the region
#y>=b/a*x), followed by the conjectured asympotics in n
#for starting at the origin for M>=a>=b>=1 and gcd(a,b)=1
#it displays the first L terms 
#SipurLineAB(3,100,20,n);
SipurLineAB:=proc(M,K,L,n) local a,b,t0,gu:
t0:=time():
print(`Below are the sequences`):
print(`enumerating the number of ways of walking from the origin`):
print(`with positive unit steps, in the square lattice`):
print(`to the point (a*i,b*i), all the while staying in the`):
print(`the region y>=b/a*x  `):
print(`followed by the conjectured asymptotics (if found) `):
print(`for 2<=b<=a<=`, M):

for a from 2 to M do
for b from 1 to a do
 if gcd(a,b)=1 then
 gu:=LineAB(a,b,K,n):
 if gu[2]<>FAIL then
  print(`a=`,a,`b=`,b, `seq.:` [op(1..L,gu[1])]):
  print(`asymp.:`, gu[2]):
else
  print(`a=`,a,`b=`,b, `seq.:` [op(1..L,gu[1])]):
 fi:
fi:
od:
od:
print(`------------------------`):
print(`This took`, time()-t0, `seconds. `):
end: