######################################################################
##ChowRobbins: Save this file as  ChowRobbins                        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read ChowRobbins<Enter>                            #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Luis A. Medina and Doron Zeilberger, Rutgers University,#
#[medina, zeilberg] at math dot rutgers dot edu                      #
######################################################################
 
#Created: June 3, 2009
 
print(`Created:  June 3, 2009`):
print(` This is the Maple package ChowRobbins. `):
print(`It is one of the packages that accompany the article `):
print(`"An Experimental Mathematics Perspective on the Old, and`):
print(`Still Open, Question of When to Stop" `):
print(`by Luis Medina and Doron Zeilberger`):
print(`available from Medina's and Zeilberger's 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 dealing with the criteria`):
print(`of maximixing the prob. of supassing a given goal`):
print(`type ezraGoal();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

 print(`For a list of the procedures dealing with general prob. dist.`):
print(`type ezraPD();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

 print(`For a list of the MAIN 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(` Bd, CRpay, IsBoundaryPt, IsInteriorPt, PayOff, W `):

else
ezra(args):
fi:

end:

ezraPD:=proc()

if args=NULL then
 print(` The procedures with general prob. distributions are: `):
 print(`Av, CRt, CRtm, CRtFast, Mom, PrExc, Sd `):

else
ezra(args):
fi:

end:

ezraGoal:=proc()

if args=NULL then
 print(` The procedures with exceeding a goal are: `):
 print(` ComparativeGambling, CRg, CRgt `):

else
ezra(args):
fi:

end:

ezraPreComputed:=proc()

if args=NULL then
 print(`The pre-comuted data procedures are:`):
 print(` Beta10000, Beta50000, CutOffs2000 `):
else
 ezra(args):
fi:
end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: BetaSeq,BetaSeqFloat,`):
 print(` CR, CRFast, CRm, CRmFloat, CRFastFloat, CutOffs, Sipur`):
 print(` SipurCG,TurningPoint `):
 print(` `):



elif nops([args])=1 and op(1,[args])=Av then
print(`Av(F,t): The average of the prob. dist. given by the p.g.f. F(t).`):
print(`For example try Av((1+t)/2,t);`):


elif nops([args])=1 and op(1,[args])=Bd then
print(`Bd(N,beta): the set of boundary points of the discrete region `):
print(`{[i,j]; 0<=i,0<=j, i+j<=N, i-j<=beta[i+j]}`):
print(`For example, try:`):
print(`Bd(4,[1,2,3]);`):

elif nops([args])=1 and op(1,[args])=Beta10000 then
print(`Beta10000():`):
print(`The (sure!, using exact arithmetic, with rational number`):
print(`sequence of pre-computed sequene of stopping times for n=1..50`):
print(`with N=10000 `):
print(`Warning: it is not rigorously proved that this is the same as`):
print(`the ultimate sequence, N=inifinty, but it is very likely so`):
print(`Do: Beta10000(); `):

elif nops([args])=1 and op(1,[args])=Beta50000 then
print(`Beta50000()`):
print(`The (probably, using floating-point`):
print(`sequence of pre-computed sequene of stopping times for n=1..186`):
print(`with N=50000 `):
print(`Warning: it is not rigorously proved that this is the same as`):
print(`the ultimate sequence, N=inifinty, but it is rather likely so`):
print(`Do: Beta50000(); `):

elif nops([args])=1 and op(1,[args])=BetaSeq then
print(`BetaSeq(N,m) the sequence of beta(N,n) the`):
print(`for n=1...m of the`):
print(`sequence such that you should stop as soon`):
print(`as you reach an n such that`):
print(`the number of #Heads- #Tails is>=beta(N,n)`):
print(`if you are allowed N coin-tosses.`):
print(`For example, try: BetaSeq(1000,10);`):

elif nops([args])=1 and op(1,[args])=BetaSeqFloat then
print(`BetaSeqFloat(N,m):like BetaSeq(N,m)`):
print(`but using Floating point.`):
print(`For example, try: BetaSeqFloat(1000,10);`):


elif nops([args])=1 and op(1,[args])=ComparativeGambling then
print(`ComparativeGambling(h,t,N,g): compares`):
print(`the the max. expectation and the Goal-oriented`):
print(`approaches to gambling if you have h heads and`):
print(`t tails, with longevity N and goal g.`):
print(`For example, try:`):
print(`ComparativeGambling(2,1,200,3/4);`):

elif nops([args])=1 and op(1,[args])=CR then
print(`CR(i,j,N): In the coin-tossing game`):
print(`if right now you have i Heads and j Tails`):
print(`and you are allowed to toss a coin up to N times`):
print(`and at any time you have the option to collect i/(i+j) or`):
print(`to keep playing, hoping to do better, but taking a chance`):
print(`of doing worse.`):
print(`What is your expected gain?`):
print(`It uses backwards induction. `):
print(`The output is the expected gain`):
print(`Warning: a much faster version is CRFast .`):
print(`For example, try:`):
print(`CR(2,1,54); `):


elif nops([args])=1 and op(1,[args])=CRg then
print(`CRg(i,j,N,g): In the coin-tossing game`):
print(`if right now you have i Heads and j Tails`):
print(`and you are allowed to toss a coin up to N times`):
print(`and at any time you have the option to collect i/(i+j) or`):
print(`to keep playing, hoping to do better, but taking a chance`):
print(`of doing worse.`):
print(`Your goal is to maximize your chance of getting>=g `):
print(`What is your probability for that?`):
print(`Of course if i/(i+j)>=g then it is 1, and you quit.`):
print(`It uses backwards induction. `):
print(`The output is the prob. of exiting the game with>=g`):
print(`Warning: a much faster version is CRFast .`):
print(`For example, try:`):
print(`CRg(2,1,54,3/4); `):

elif nops([args])=1 and op(1,[args])=CRgt then
print(`CRgt(i,j,N,g,t): Like CRg(i,j,N,g) `):
print(`but second component gives the prob. gen. function of gain`):
print(`For example, try:`):
print(`CRgt(2,1,54,3/4,t); `):


elif nops([args])=1 and op(1,[args])=CRpay then
print(`CRpay(i,j,N): same as CR(i,j,N)`):
print(`using PayOff `):
print(`For example, try:`):
print(`CRpay(2,1,54); `):

elif nops([args])=1 and op(1,[args])=CRt then
print(`CRt(i,j,N,t): like CR(i,j,N)`):
print(`but returns in addition to the expectation`):
print(`also the prob. gen. function`):
print(`Warning: a much faster version is CRFast .`):
print(`For example, try:`):
print(`CRt(2,1,54,t); `):

elif nops([args])=1 and op(1,[args])=CRtFast then
print(`CRtFast(N,K1,K2,t) the table of CRt(i,j,N,t)`):
print(`from i+j=K1 to i+j=K2. In particular, to get the whole`):
print(`table., type CRFast(N,N,0); `):
print(`For example, try: `):
print(`CRtFast(10,10,0,t);`):
print(`CRtFast(10,1,0,t);`):

elif nops([args])=1 and op(1,[args])=CRm then
print(`CRm(i,j,N): faster version of CR(i,j,N).`):
print(`For example, try:`):
print(`CRm(2,1,54); `):

elif nops([args])=1 and op(1,[args])=CRtm then
print(`CRtm(i,j,N): faster version of CRt(i,j,N,t).`):
print(`For example, try:`):
print(`CRtm(2,1,54,t); `):

elif nops([args])=1 and op(1,[args])=CRmFloat then
print(`CRmFloat(i,j,N): faster version of CRf(i,j,N).`):
print(`For example, try:`):
print(`CRmFloat(2,1,54); `):

elif nops([args])=1 and op(1,[args])=CRFast then
print(`CRFast(N,K1,K2) the table of CR(i,j,N)`):
print(`from i+j=K1 to i+j=K2. In particular, to get the whole`):
print(`table., type CRFast(N,N,0); `):
print(`For example, try: `):
print(`CRFast(10,10,0);`):
print(`CRFast(10,1,0);`):

elif nops([args])=1 and op(1,[args])=CRFastFloat then
print(`CRFastFloat(N,K1,K2): a (faster) floating-point`):
print(`version of CRFast(N,K1,K2)`);
print(`For example, try: `):
print(`CRFastFloat(10,10,0);`):
print(`CRFastFloat(10,1,0);`):


elif nops([args])=1 and op(1,[args])=CutOffs then
print(`CutOffs(N,m) the sequence of cutsoffs (h,t) `):
print(`for each h+t=k, k<=m, (the smallest number of heads where `):
print(`it is a GO with longevity N, and the first N1 where it is a GO.`):
print(`For example, try: CutOffs(1000,10);`):

elif nops([args])=1 and op(1,[args])=CutOffs2000 then
print(`CutOffs2000(): The sequence of size 100 whose i-th entry`):
print(`is the pair [[h,t],HowOld]  where [h,t] is the position with`):
print(`smallest h such that h+t=i and for N=20000 it is a GO position`):
print(`(so it is definitely a GO position for N=infinity, but there`):
print(`is a (very slight) chance that [h+1,i-h-1] is also a GO position`):
print(`for N=infinity (but for N=2000, it is definitly a STOP) `):
print(`HowOld is the first N such that for N=HowOld it is GO`):
print(`but for N<HowOld it is a STOP. `):
print(`It took  55552.319 seconds of CPU time to generate it.`):
print(`Do: CutOffs2000(); `):


elif nops([args])=1 and op(1,[args])=IsBoundaryPt then
print(`IsBoundaryPt(pt,N,beta): Is the point pt on`):
print(`the boundary of the discrete region`):
print(`(0<=i+j<=N, i-j<=beta[i+j])`):
print(`For example try: IsBoundaryPt([2,0],3,[1,2]);`):

elif nops([args])=1 and op(1,[args])=IsInteriorPt then
print(`IsInteriorPt(pt,N,beta): Is the point pt in`):
print(`the interior of the discrete region`):
print(`(0<=i+j<=N, i-j<=beta[i+j])`):
print(`For example try: IsInteriorPt([2,0],3,[1,2]);`):

elif nops([args])=1 and op(1,[args])=Mom then
print(`Mom(F,t,i): The ith-moment about the mean `):
print(` of the prob. dist. given by the p.g.f. F(t).`):
print(`For example try Mom((1+t)/2,t,2);`):


elif nops([args])=1 and op(1,[args])=PayOff then
print(`PayOff(pt,N,beta): the pay-off, at point pt, with the boundary`):
print(`given by N,beta. For example, try:`):
print(`PayOff([0,0],5,[1,2,3,4]);`):

elif nops([args])=1 and op(1,[args])=PrExc then
print(`PrExc(F,t,g): Given a prob. gen. function F of t`):
print(`finds the probability of getting >=g. For example, try`):
print(`PrExc((1+t+t^2)/3,t,1);`):

elif nops([args])=1 and op(1,[args])=Sd then
print(`Sd(F,t): The standard deviation`):
print(` of the prob. dist. given by the p.g.f. F(t).`):
print(`For example try Sd((1+t)/2,t);`):

elif nops([args])=1 and op(1,[args])=Sipur then
print(`Sipur(N,m,r): All the info about expected gain`):
print(`for positions with i+j<=m living up to N tosses`):
print(`with resolution r for prob.`):
print(`For example, try:`):
print(`Sipur(200,3,.1);`):


elif nops([args])=1 and op(1,[args])=SipurCG then
print(`SipurCG(N,m,r): All the info about comparative gambling with`):
print(`various goals`):
print(`for positions with i+j<=m living up to N tosses`):
print(`with resolution r for prob.`):
print(`For example, try:`):
print(`SipurCG(200,3,.1);`):

elif nops([args])=1 and op(1,[args])=TurningPoint then
print(`TurningPoint(i,j,N0): The first N<=N0 (if it exists)`):
print(`such that i Heads and j Tails starts being GO rather`):
print(`than Stop. For example, try:`):
print(`TurningPoint(2,1,100);`):

elif nops([args])=1 and op(1,[args])=W then
print(`W(pt1,pt2,N,beta): The number of walks from pt1 to pt2`):
print(`staying in the region i+j<N i-j<beta[i+j] and quitting`):
print(`as soon as it hits the boundary i+j=N or i-j=beta[i+j].`):
print(`For example try: W([0,0],[1,1],3,[1,2]);`):

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

ez:=proc(): print(`CR(i,j,N), CRfa(N,m), CRfaf(N,m) `): 
print(`Tria(N), CutOff(N), CutOffU(M,L1,L2), CRh(i,j,N)  `):
print(` CRhFast(N,K1,K2), CRFast(N,K1,K2) `):
print(`CutOffFast(N), Pd(i,j,N,t), Moments(i,j,N,M) `):
print(`Tran(i,j,N), TranH(i,j,N), PolToPd(P,t), PolToCd(P,t) `):
print(`CRp(i,j,N,g), CRf(i,j,N), CRm(i,j,N), CRmF(i,j,N) `):
print(` CRmH(i,j,N) , TriaH(N) , IsStop(i,j,N), TriaS(N,S)  `):
print(`GuessRat1(gu,L,n,d), GuessRat(gu,L,n) `):
print(`GuessAB(N,A,B,K)`):
end:


#Aitken(L): The conjectured limit of the sequence L
#from its three last elements, using Aitken's fromula
Aitken:=proc(L)
local n,a,b,c:
n:=nops(L):

if n<3 then
 RETURN(FAIL):
fi:

a:=L[n-2]: b:=L[n-1]: c:=L[n]:

c-(c-b)^2/(c-2*b+a):
end:

#CR(i,j,N): backwards induction for the coin-tossing game
#when you have to quit after N steps
CR:=proc(i,j,N)  option remember:

if i+j>N then
 RETURN(FAIL):
fi:

if i+j=N then
 RETURN(max(1/2,i/N)):

elif i+j>0 then
 RETURN(max( (CR(i+1,j,N)+CR(i,j+1,N))/2, i/(i+j) )):
else
 RETURN(max( (CR(i+1,j,N)+CR(i,j+1,N))/2)):
fi:

end:



#CRf(i,j,N): backwards induction for the coin-tossing game
#when you have to quit after N steps
CRf:=proc(i,j,N)  option remember:

if i+j>N then
 RETURN(FAIL):
fi:

if i+j=N then
 RETURN(max(.5,evalf(i/N))):

elif i+j>0 then
 RETURN(evalf(max( (CR(i+1,j,N)+CR(i,j+1,N))/2, i/(i+j) ))):
else
 RETURN(evalf(max( (CR(i+1,j,N)+CR(i,j+1,N))/2))):
fi:

end:







#Tria(N): the triangle with N
Tria:=proc(N) local j,i1:
[seq([seq(CR(i1,j-i1,N),i1=0..j)],j=0..N)]:
end:


#TriaH(N): the triangle with N
TriaH:=proc(N) local j,i1:
[seq([seq(CRh(i1,j-i1,N),i1=0..j)],j=0..N)]:
end:

#CutOff1(L): the largest j such that CR(i,j,nops(L)-1) d
CutOff1:=proc(L) local m,i:
m:=nops(L)-1:

 if m=0 then
  RETURN(0):
 fi:

for i from 1 while L[i]<>(i-1)/m do od:

i-1:

end:



CutOff:=proc(N) local gu,i:
gu:=Tria(N):
[seq(CutOff1(gu[i]),i=1..nops(gu))]:
end:


#CutOffU(M,L1,L2): the first M terms of the ultimate cutoffs
#Cutoff(infinity)
CutOffU:=proc(M,L1,L2) local gu,i,j:


for i from 2*M while nops({seq([op(1..M+1,CutOff(i+j))],j=0..L1)})<>1 do od:

gu:=[op(1..M+1,CutOff(i+L1))]:


if nops({seq([op(1..M+1,CutOff(i+j))],j=0..L2)})<>1 then
 RETURN(FAIL):
else
 RETURN(gu):
fi:

end:







#CRh(i,j,N): backwards induction for the honest
#version of the coin-tossing game
#when you have to quit after N steps with #H/N
CRh:=proc(i,j,N)  option remember:

if i+j>N then
 RETURN(FAIL):
fi:

if i+j=N then
 RETURN(i/N):

elif i+j>0 then
 RETURN(max( (CRh(i+1,j,N)+CRh(i,j+1,N))/2, i/(i+j) )):
else
 RETURN(max( (CRh(i+1,j,N)+CRh(i,j+1,N))/2)):
fi:

end:

#CRhFast(N,K1,K2) the table of CRh(i,j,N) only for i+j<=K
CRhFast:=proc(N,K1,K2) local gu,k,i,mu:
gu:=[seq(i/N,i=0..N)]:

for k from N-1 by -1 to K1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:
od:

mu:=[gu]:

for k from K1-1 by -1 to max(1,K2) do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:
mu:=[op(mu),gu]:
od:

if K2=0 then
gu:=[seq(max( (gu[i]+gu[i+1])/2),i=1..nops(gu)-1)]:
mu:=[op(mu),gu]:
fi:
mu:
end:



#CRFast(N,K1,K2) the table of CR(i,j,N) only for i+j<=K
CRFast:=proc(N,K1,K2) local gu,k,i,mu:
option remember:
gu:=[seq(max(1/2,i/N),i=0..N)]:

for k from N-1 by -1 to K1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:
od:

mu:=[gu]:

for k from K1-1 by -1 to max(1,K2) do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:
mu:=[op(mu),gu]:
od:

if K2=0 then
gu:=[seq(max( (gu[i]+gu[i+1])/2),i=1..nops(gu)-1)]:
mu:=[op(mu),gu]:
fi:
mu:
end:


#CRFastFloat(N,K1,K2) the table of CR(i,j,N) only for i+j<=K
#in Floating-point
CRFastFloat:=proc(N,K1,K2) local gu,k,i,mu:
option remember:
gu:=evalf([seq(max(.5,i/N),i=0..N)]):

for k from N-1 by -1 to K1 do
gu:=evalf([seq(max( (gu[i]+gu[i+1])/2. ,(i-1)/k),i=1..nops(gu)-1)]):
od:

mu:=[gu]:

for k from K1-1 by -1 to max(1,K2) do
gu:=evalf([seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]):
mu:=[op(mu),gu]:
od:

if K2=0 then
gu:=evalf([seq(max( (gu[i]+gu[i+1])/2),i=1..nops(gu)-1)]):
mu:=[op(mu),gu]:
fi:
mu:
end:






#CutOffFast(N): Like CutOff(N):
CutOffFast:=proc(N) local gu,k,i,mu:
gu:=[seq(max(1/2,i/N),i=0..N)]:
mu:=[CutOff1(gu)]:
for k from N-1 by -1 to 1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:
mu:=[op(mu),CutOff1(gu)]:
od:

gu:=[seq(max( (gu[i]+gu[i+1])/2),i=1..nops(gu)-1)]:
mu:=[op(mu),CutOff1(gu)]:

[seq(mu[nops(mu)-i+1],i=1..nops(mu))]:
end:





#Pd(i,j,N,t): the probability distribution (in terms of a g.f. in t)
#if you quit after N throws, and right now you have i Heads and j Tails
Pd:=proc(i,j,N,t)  option remember:

if i+j=N then
 RETURN(t^CR(i,j,N)):
fi:

if i+j=0 then
 RETURN( (Pd(i+1,j,N,t)+Pd(i,j+1,N,t))/2):
fi:
if CR(i,j,N)=i/(i+j) then
  RETURN(t^(i/(i+j))):
else
 RETURN( (Pd(i+1,j,N,t)+Pd(i,j+1,N,t))/2):
fi:

end:


#Moments(i,j,N,M): the expectation, variance, and
#list of normalized moments,up to the Mth.
#for coin-tossing if you have i Heads and j Tails
Moments:=proc(i,j,N,M) local gu,f,t,av,va,i1:
f:=Pd(i,j,N,t):
av:=CR(i,j,N):
gu:=[av]:
f:=f/t^av:

f:=t*diff(f,t):
f:=t*diff(f,t):
va:=subs(t=1,f):

if va=0 then
 RETURN([av,0$(M-1)]):
fi:

gu:=[op(gu),va]:

for i1 from 1 to M-2 do
 f:=t*diff(f,t):
 gu:=[op(gu),evalf(subs(t=1,f)/va^((i1+2)/2))]:
od:

gu:

end:






#Tran(i,j,N): the smallest number of coin-tosses such that
#i Heads and j Tails starts to be a GO <=N. If with N it is a stop
#it returns FAIL
Tran:=proc(i,j,N) local M:

for M from i+j to N while CRm(i,j,M)=i/(i+j) do od:

if M=N+1 then
  RETURN(FAIL):
else
 RETURN(M):
fi:



end:


#PolToPd(P,t): Given a polynomial P(t) translates it to a discrete
#prob. distibution. For example, try: PolToPd((t^(1/2)+t^(2/3))/2,t);
PolToPd:=proc(P,t) local P1,i,mu,ev,po,c,T:
P1:=expand(P):

mu:=[]:
for i from 1 to nops(P1) do
ev:=op(i,P1):
c:=subs(t=1,ev):
po:=subs(t=1,normal(diff(ev,t)/ev)):
mu:=[op(mu),po]:
T[po]:=c:
od:

mu:=sort(mu):


[seq([mu[i],T[mu[i]]],i=1..nops(mu))]:

end:


#PolToCd(P,t): Given a polynomial P(t) translates it to a 
#com. discrete
#prob. distibution. For example, try: PolToPd((t^(1/2)+t^(2/3))/2,t);
PolToCd:=proc(P,t) local gu,i,j:
gu:=PolToPd(P,t):

[seq([gu[i][1], add(gu[j][2],j=1..i)],i=1..nops(gu))]:

end:




#CRp(i,j,N,g): backwards induction for the coin-tossing game
#when you have to quit after N steps and your criterion is
#to maximize your chances of getting at least g
CRp:=proc(i,j,N,g)  local gu1,gu2:
option remember:

if g<=1/2 then
print(`g must be bigger than 1/2 `):
fi:

if i+j>N then
 RETURN(FAIL):
fi:

if i+j=N then
if i/N>=g then
  RETURN([1,i/N]):
else
  RETURN([0,max(1/2,i/N)]):
fi:

elif i+j>0 then

 if i/(i+j)>=g then
   RETURN([1,i/(i+j)]):
else
  gu1:=CRp(i+1,j,N,g):   gu2:=CRp(i,j+1,N,g):
  RETURN([(gu1[1]+gu2[1])/2,(gu1[2]+gu2[2])/2] ):
 fi:

elif i+j=0 then

  gu1:=CRp(i+1,j,N,g):   gu2:=CRp(i,j+1,N,g):
  RETURN([(gu1[1]+gu2[1])/2,(gu1[2]+gu2[2])/2] ):

else
 ERROR(`Something is wrong`):
fi:

end:


#TranH(i,j,N): the smallest number of coin-tosses such that
#i Heads and j Tails starts to be a GO <=N. If with N it is a stop
#it returns FAIL
TranH:=proc(i,j,N) local M:

for M from i+j to N while CRh(i,j,M)=i/(i+j) do od:

if M=N+1 then
  RETURN(FAIL):
else
 RETURN(M):
fi:



end:



CRm:=proc(i,j,N)
option remember:
if i+j>0 then
CRFast(N,i+j,i+j)[1][i+1]
else
(CRm(0,1,N)+CRm(1,0,N))/2:
fi:

end:

CRmFloat:=proc(i,j,N)
option remember:
CRFastFloat(N,i+j,i+j)[1][i+1]
end:






CRmH:=proc(i,j,N):
CRhFast(N,i+j,i+j)[1][i+1]
end:







#IsStop(i,j,N): is (i,j) stop?
IsStop:=proc(i,j,N) 
evalb(CRm(i,j,N)=i/(i+j)):
end:




#TriaS(N,S): the triangle with N
TriaS:=proc(N,S) local j,i1,gu,gu1:
gu:=[]:

for j from 1 to N do
gu1:=[]:
 for i1 from 0 to j do
  if IsStop(i1,j-i1,N) then
    gu1:=[op(gu1),S]:
  else
   gu1:=[op(gu1), CR(i1,j-i1,N)]:
 fi:
 od:
gu:=[op(gu),gu1]:
od:
gu:
end:



#GuessRat1(gu,L,n,d): given a list gu 
#gueses a rational function R(n) of degree <=d , such that gu[i]=R(i)
#for i>=L
#For example, try: GuessRat1([seq(i/(i+1),i=1..20)],3,n,1);
GuessRat1:=proc(gu,L,n,d) local eq,var,a,i,b,R:

if nops(gu)-L<=2*d+2+5 then
 print(`make list bigger`):
 FAIL:
fi:

var:={seq(a[i],i=0..d-1), seq(b[i],i=0..d)}:
R:=add(a[i]*n^i,i=0..d)/add(b[i]*n^i,i=0..d):
R:=subs(a[d]=1,R):
eq:={ seq(subs(n=i,R)-gu[i],i=L..nops(gu))}:

var:=solve(eq,var):

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

normal(subs(var,R)):

end:


#GuessRat(gu,L,n): given a list gu 
#gueses a rational function R(n) of degree <=d , such that gu[i]=R(i)
#for i>=L
#For example, try: GuessRat([seq(i/(i+1),i=1..20)],3,n);
GuessRat:=proc(gu,L,n) local R,d:

for d from 1 to (nops(gu)-L-7)/2 do
 R:=GuessRat1(gu,L,n,d):
 if R<>FAIL then
  RETURN(R):
 fi:
od:
FAIL:
end:


 
#GuessAB(N,A,B,K): guesses a rational function for CR(N-A,N-B,2*N)
#for numeric A and B and symbolic N, (K is the size of thedata set) 
#nd the place where
#it starts to apply
#For example, try:
#GuessAB(N,1,1,20);
GuessAB:=proc(N,A,B,K) local gu,C,R,i,N1:
C:=max(A,B):
gu:=[seq(CRm(N1-A,N1-B,2*N1),N1=C+3..C+K)]:


R:=GuessRat(gu,K/2,N):



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

R:=normal(subs(N=N-C-2,R)):

for i from C+K/2 by -1 to C while subs(N=i,R)=CRm(i-A,i-B,2*i) do od:

RETURN(R,i):


end:





#BetaSeq(N,m) the sequence of beta(N,n) the
#for n=1...m of the
#sequence such that you should stop as soon
#as you reach an n such that
# #Heads-#Tails is>=beta(N,n)
#if you are allowed N coin-tosses.
#For example, try: BetaSeq(1000,10);
BetaSeq:=proc(N,m) local gu,k,i,j,T:
if m>=N then
 ERROR(`Second argument must be strictly smaller than first argument`):
fi:

gu:=[seq(max(1/2,i/N),i=0..N)]:

for k from N-1 by -1 to m+1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:
od:

for k from m by -1 to  1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:

for j from 0 to nops(gu)-1 while gu[j+1]<>j/k do od:
 T[k]:=j:
od:


[seq(2*T[j]-j,j=1..m)]:


end:


#BetaSeqFloat(N,m) the sequence of beta(N,n) the
#like BetaSeq(N,m) but using floating point
#For example, try: BetaSeqFloat(1000,10);
BetaSeqFloat:=proc(N,m) local gu,k,i,j,T:
gu:=[seq(max(.5,evalf(i/N)),i=0..N)]:

for k from N-1 by -1 to m+1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2. ,evalf((i-1)/k)),i=1..nops(gu)-1)]:
od:

for k from m by -1 to  1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2. ,evalf((i-1)/k)),i=1..nops(gu)-1)]:

for j from 0 to nops(gu)-1 while abs(gu[j+1]-evalf(j/k))>10^(-Digits+4) do od:
 T[k]:=j:
od:


[seq(2*T[j]-j,j=1..m)]:

end:


#TurningPoint(i,j,N0): The first N<=N0 (if it exists)
#such that i Heads and j Tails starts being GO rather
#than Stop. For example, try:
#TurningPoint(2,1,100);
TurningPoint:=proc(i,j,N0) local N:

for N from i+j+1 to N0 while CR(i,j,N)=i/(i+j) do od:
N:

end:




#CutOffs(N,m) the sequence of cutsoffs (h,t) with h+t<=m
#for each h+t=k (k<=m), and the first N1 where it is a GO
#with life-span N
#For example, try: CutOffs(1000,10);
CutOffs:=proc(N,m) local gu,k,i,j,aluf,mu:
gu:=[seq(max(1/2,i/N),i=0..N)]:
mu:=[]:

for k from N-1 by -1 to m+1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:
od:

for k from m by -1 to 1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:

for j from 0 to nops(gu)-1 while gu[j+1]<>j/k do od:
aluf:=[j-1,k-j+1]:
mu:=[op(mu),[aluf,TurningPoint(j-1,k-j+1,N)]]:
od:
[seq(mu[nops(mu)-i+1],i=1..nops(mu))]:
end:



#TurningPointM(i,j,N0): The first N<=N0 (if it exists)
#such that i Heads and j Tails starts being GO rather
#than Stop. For example, try:
#TurningPointM(2,1,100);
TurningPointM:=proc(i,j,N0) local N:

for N from i+j+1 to N0 while CRm(i,j,N)=i/(i+j) do od:
N:

end:




#CutOffsM(N,m) the sequence of cutsoffs (h,t) with h+t<=m
#for each h+t=k (k<=m), and the first N1 where it is a GO
#with life-span N
#For example, try: CutOffs(1000,10);
CutOffsM:=proc(N,m) local gu,k,i,j,aluf,mu:
gu:=[seq(max(1/2,i/N),i=0..N)]:
mu:=[]:

for k from N-1 by -1 to m+1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:
od:

for k from m by -1 to 1 do
gu:=[seq(max( (gu[i]+gu[i+1])/2 ,(i-1)/k),i=1..nops(gu)-1)]:

for j from 0 to nops(gu)-1 while gu[j+1]<>j/k do od:
aluf:=[j-1,k-j+1]:
mu:=[op(mu),[aluf,TurningPointM(j-1,k-j+1,N)]]:
od:
[seq(mu[nops(mu)-i+1],i=1..nops(mu))]:
end:

Beta10000:=proc():
[1, 2, 3, 2, 3, 2, 3, 2, 3, 4, 3, 4, 3, 4, 3, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4,

    5, 4, 5, 6, 5, 6, 5, 6, 5, 6, 5, 6, 5, 6, 5, 6, 7, 6, 7, 6, 7, 6, 7, 6]:
end:


Beta50000:=proc():
[1, 2, 3, 2, 3, 2, 3, 2, 3, 4, 3, 4, 3, 4, 3, 4, 5, 4, 5, 4, 5, 4, 5, 4, 5, 4,

    5, 4, 5, 6, 5, 6, 5, 6, 5, 6, 5, 6, 5, 6, 5, 6, 7, 6, 7, 6, 7, 6, 7, 6, 7,

    6, 7, 6, 7, 6, 7, 6, 7, 8, 7, 8, 7, 8, 7, 8, 7, 8, 7, 8, 7, 8, 7, 8, 7, 8,

    7, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9, 8, 9,

    10, 9, 10, 9, 10, 9, 10, 9, 10, 9, 10, 9, 10, 9, 10, 9, 10, 9, 10, 9, 10,

    9, 10, 9, 10, 9, 10, 11, 10, 11, 10, 11, 10, 11, 10, 11, 10, 11, 10, 11,

    10, 11, 10, 11, 10, 11, 10, 11, 10, 11, 10, 11, 10, 11, 12, 11, 12, 11, 12,

    11, 12, 11, 12, 11, 12, 11, 12, 11, 12, 11, 12, 11, 12, 11, 12, 11, 12, 11,

    12, 11, 12, 11, 12, 11, 12]:
end:




#CutOffs2000(): The sequence of size 100 whose i-th entry
#is the pair [[h,t],HowOld]  where [h,t] is the position with
#smallest h such that h+t=i and for N=20000 it is a GO position
#(so it is definitely a GO position for N=infinity, but there
#is a (very slight) chance that [h+1,i-h-1] is also a GO position
#for N=infinity (but for N=2000, it is definitly a STOP) 
#HowOld is the first N such that for N=HowOld it is GO
#but for N<HowOld it is a STOP. 
#It took  55552.319 seconds of CPU time to generate it.
#Do: CutOffs2000(); 
CutOffs2000:=proc():
[[[0, 1], 2], [[1, 1], 3], [[2, 1], 51], [[2, 2], 5], [[3, 2], 7], [[3, 3], 7],

    [[4, 3], 9], [[4, 4], 9], [[5, 4], 11], [[6, 4], 35], [[6, 5], 13],

    [[7, 5], 23], [[7, 6], 15], [[8, 6], 21], [[8, 7], 17], [[9, 7], 21],

    [[10, 7], 1421], [[10, 8], 23], [[11, 8], 91], [[11, 9], 25], [[12, 9], 57],

    [[12, 10], 25], [[13, 10], 47], [[13, 11], 27], [[14, 11], 43],

    [[14, 12], 29], [[15, 12], 43], [[15, 13], 31], [[16, 13], 43],

    [[17, 13], 277], [[17, 14], 43], [[18, 14], 139], [[18, 15], 43],

    [[19, 15], 103], [[19, 16], 45], [[20, 16], 87], [[20, 17], 45],

    [[21, 17], 79], [[21, 18], 47], [[22, 18], 75], [[22, 19], 49],

    [[23, 19], 73], [[24, 19], 1679], [[24, 20], 71], [[25, 20], 423],

    [[25, 21], 71], [[26, 21], 249], [[26, 22], 69], [[27, 22], 185],

    [[27, 23], 69], [[28, 23], 155], [[28, 24], 71], [[29, 24], 137],

    [[29, 25], 71], [[30, 25], 125], [[30, 26], 73], [[31, 26], 119],

    [[31, 27], 73], [[32, 27], 113], [[32, 28], 75], [[33, 28], 109],

    [[34, 28], 833], [[34, 29], 107], [[35, 29], 477], [[35, 30], 107],

    [[36, 30], 343], [[36, 31], 105], [[37, 31], 275], [[37, 32], 105],

    [[38, 32], 235], [[38, 33], 105], [[39, 33], 211], [[39, 34], 105],

    [[40, 34], 193], [[40, 35], 105], [[41, 35], 181], [[41, 36], 105],

    [[42, 36], 171], [[42, 37], 105], [[43, 37], 165], [[43, 38], 107],

    [[44, 38], 159], [[45, 38], 1039], [[45, 39], 155], [[46, 39], 679],

    [[46, 40], 153], [[47, 40], 513], [[47, 41], 151], [[48, 41], 419],

    [[48, 42], 149], [[49, 42], 361], [[49, 43], 147], [[50, 43], 321],

    [[50, 44], 147], [[51, 44], 293], [[51, 45], 147], [[52, 45], 271],

    [[52, 46], 145], [[53, 46], 255], [[53, 47], 145]]:
end:



####Begin general prob. dist.

#Av(F,t): The average of the prob. dist. given by the p.g.f. F(t)
#For example try Av((1+t)/2,t);
Av:=proc(F,t):
if subs(t=1,F)<>1 then
ERROR(`Not prob. dist.`):
fi:
subs(t=1,diff(F,t)):
end:


#Sd(F,t): The standard deviation
# of the prob. dist. given by the p.g.f. F(t)
#For example try Sd((1+t)/2,t);
Sd:=proc(F,t):
if subs(t=1,F)<>1 then
ERROR(`Not prob. dist.`):
fi:
sqrt(subs(t=1,diff(F,t$2))+Av(F,t)-Av(F,t)^2):
end:


#Mom(F,t,i): the i-th moment of the prob. dist. given by 
#the prob. gen. fun. F(t). 
#For example try Mom((1+t)/2,t,5);
Mom:=proc(F,t,i) local gu,i1:
if subs(t=1,F)<>1 then
ERROR(`Not prob. dist.`):
fi:
gu:=F/t^Av(F,t):
for i1 from 1 to i do
gu:=expand(t*diff(gu,t)):
od:
subs(t=1,gu):
end:



#CRt(i,j,N,t): backwards induction for the coin-tossing game
#when you have to quit after N steps, also gives the
#prob. dist.
CRt:=proc(i,j,N,t)  option remember:

if i+j>N then
 RETURN(FAIL):
fi:

if i+j=N then
 RETURN([max(1/2,i/N), t^max(1/2,i/N)]):

elif i+j>0 then
   if (CRt(i+1,j,N,t)[1]+CRt(i,j+1,N,t)[1])/2 >i/(i+j) then
     RETURN([(CRt(i+1,j,N,t)[1]+CRt(i,j+1,N,t)[1])/2,
    expand((CRt(i+1,j,N,t)[2]+CRt(i,j+1,N,t)[2])/2)]):
   else
   RETURN([i/(i+j),t^(i/(i+j))]):
  fi:

else
 RETURN(
[(CRt(1,0,N,t)[1]+CRt(0,1,N,t)[1])/2,
expand((CRt(1,0,N,t)[2]+CRt(0,1,N,t)[2])/2)] ):
fi:

end:



#CRtFast(N,K1,K2,t) the table of CRt(i,j,N,t) only for i+j<=K
CRtFast:=proc(N,K1,K2,t) local gu,k,i,mu,lu:
option remember:
gu:=[seq([max(1/2,i/N),t^max(1/2,i/N)],i=0..N)]:

for k from N-1 by -1 to K1 do
lu:=[]:
for i from 1 to nops(gu)-1 do
 if (gu[i][1]+gu[i+1][1])/2>(i-1)/k then
   lu:=[op(lu),[(gu[i][1]+gu[i+1][1])/2,(gu[i][2]+gu[i+1][2])/2]]:
 else
  lu:=[op(lu),[(i-1)/k,t^((i-1)/k)] ]:
 fi
od:
gu:=lu:
od:

mu:=[gu]:

for k from K1-1 by -1 to max(1,K2) do
lu:=[]:
for i from 1 to nops(gu)-1 do
 if (gu[i][1]+gu[i+1][1])/2>(i-1)/k then
   lu:=[op(lu),[(gu[i][1]+gu[i+1][1])/2,(gu[i][2]+gu[i+1][2])/2]]:
 else
   lu:=[op(lu),[(i-1)/k,t^((i-1)/k)] ]:
  fi:
od:
gu:=lu:
mu:=[op(mu),gu]:
od:

if K2=0 then
gu:=[seq(
[(gu[i][1]+gu[i+1][1])/2,(gu[i][2]+gu[i+1][2])/2],
i=1..nops(gu)-1 )]:
mu:=[op(mu),gu]:
fi:
mu:
end:

CRtm:=proc(i,j,N,t)
option remember:
if i+j>0 then
CRtFast(N,i+j,i+j,t)[1][i+1]
else
[(CRtm(0,1,N,t)[1]+CRtm(1,0,N,t)[1])/2,(CRtm(0,1,N,t)[2]+CRtm(1,0,N,t)[2])/2]:
fi:


end:


#PrExc(F,t,g): Given a prob. gen. function F of t
#finds the probability of getting >=g. For example, try
#PrExc((1+t+t^2)/3,t,1);
PrExc:=proc(F,t,g) local F1,gu,i,F11,c1,c2:
F1:=expand(F):

if type(F1,`^`) then
  F11:=op(2,F1):
  if F11>=g then
    RETURN(1):
  else
   RETURN(0):
  fi:
fi:

  

if not type(F1,`+`) then
RETURN(FAIL):
fi:

gu:=0:

for i from 1 to nops(F1) do
 F11:=op(i,F1):
 if not (type(F11,`*`) or type(F11,numeric)) then
    RETURN(FAIL):
 fi:

if not type(F11,numeric) then

 c1:=op(1,F11):
 c2:=op(2,F11):



 if not type(c1,numeric)  then
   RETURN(FAIL):
 fi:

if type(c2,`^`)  then
 if op(2,c2)>=g then
   gu:=gu+c1:
 fi:
elif c2=t then
   gu:=gu+c1:
else
 RETURN(FAIL):
fi:


fi:

od:

gu:
end:




####End general prob. dist.



####Begin goal-oriented
#CRg(i,j,N,g): backwards induction for the coin-tossing game
#with maximizing your chances to exceed g
#when you have to quit after N steps
CRg:=proc(i,j,N,g)  option remember:

if g<=1/2 then
print(`The last argument must be> 1/2 `):
 RETURN(FAIL):
fi:

if i=0 and j=0 then
 RETURN((CRg(i+1,j,N,g)+ CRg(i,j+1,N,g))/2):
fi:

if i+j>N then
 RETURN(FAIL):
fi:

if i+j=N then
  if i/N>=g then
     RETURN(1):
   else
     RETURN(0):
   fi:

else
 if i/(i+j)>=g then
     RETURN(1):
 else
 RETURN((CRg(i+1,j,N,g)+ CRg(i,j+1,N,g))/2):
 fi:

fi:

end:


####Begin goal-oriented
#CRgt(i,j,N,g,t): backwards induction for the coin-tossing game
#with maximizing your chances to exceed g
#also returns the prob. gen. fun.
#when you have to quit after N steps
CRgt:=proc(i,j,N,g,t)  local gu1,gu2:
option remember:

if g<=1/2 then
print(`The last argument must be> 1/2 `):
 RETURN(FAIL):
fi:

if i+j>N then
 RETURN(FAIL):
fi:


if i+j=0 then
gu1:=CRgt(1,0,N,g,t):
gu2:=CRgt(0,1,N,g,t):
RETURN([(gu1[1]+gu2[1])/2,(gu1[2]+gu2[2])/2]):
fi:

if i+j=N then
  if i/N>=g then
     RETURN([1,t^(i/N)]):
   else
     RETURN([0,1]):
   fi:

else
 if i/(i+j)>=g then
     RETURN([1,t^(i/(i+j))]):
 else

gu1:=CRgt(i+1,j,N,g,t):
gu2:=CRgt(i,j+1,N,g,t):

 RETURN([(gu1[1]+gu2[1])/2, (gu1[2]+gu2[2])/2]):

 fi:

fi:

end:


#ComparativeGambling(h,t,N,g): compares
#the the max. expectation and the Goal-oriented
#approaches to gambling if you have h heads and
#t tails, with longevity N and goal g
#For example, try:
#ComparativeGambling(2,1,200,3/4);
ComparativeGambling:=proc(i,j,N,g) local t,gu1,gu2:
if i>0 and j=0 then
 print(`Of course you should quit!`):
print(`ExpectedGain (FirstWay):1`,
`ExpectedGain (SecondWay):1`):
print(`Prob.Exceeding Goal(MaxExp): 1`,`Prob.Exc.Goal(ExcGoal):1`):
RETURN():
fi:

if g<=1/2 then
print(`Last argument must be>1/2 `):
RETURN(FAIL):
fi:
gu1:=CRt(i,j,N,t):
gu2:=CRgt(i,j,N,g,t):
print
(`ExpectedGain (FirstWay):`, evalf(gu1[1]),
`ExpectedGain (SecondWay):`,evalf(Av(gu2[2],t))):
print(`Prob.Exceeding Goal(FirstWay):`, evalf(PrExc(gu1[2],t,g)), `Prob.Exc.Goal(SecondWay):`, evalf(gu2[1])    ):


end:


####End goal-oriented

### strart stories



#Sipur(N,m,r): All the info about expected gain
#for positions with i+j<=m living up to N tosses
#with resolution r for prob.
#For example, try:
#Sipur(200,3,.1);

Sipur:=proc(N,m,r) local t,i,k,S,gu,h,t0:
t0:=time():
S:=BetaSeq(N,m):

print(`This is the story of expected gain, and stat. information`):
print(`for the Chow-Robbins coin-tossing game with life duration`, N):
print(`for the first`, m, `tosses `):
print(`with prob. of getting at least a certain amount with resolution`, r):

gu:=CRtm(0,0,N,t):

print(`The expected gain  at the very beginning is`, evalf(gu[1])):
print(`The standard deviation is`, evalf(Sd(gu[2],t))):

print(`The prob. of getting at least`, .5+r*i, ` for i from 1 to`, trunc(.5/r)):
print(`are: `):
print( evalf([seq([.5+r*i, PrExc(gu[2],t,.5+r*i)],i=1..trunc(.5/r))])):


for k from 1 to m do
print(`If the total number of tosses is`, k, `then whenever the`):
print(`number of heads minus the number of tails is >=`, S[k] ):
print(`in other words, the number of heads is >=`, (S[k]+k)/2 ):
print(`it is a stop, and you collect for sure h/(h+t) `):
print(`for the remaining case we have`):

  for h from 0 to (S[k]+k)/2-1 do
   
gu:=CRtm(h,k-h,N,t):

print(`The expected gain with`,h, `heads and `, k-h, `tails is`, evalf(gu[1])):
print(`The standard deviation is`, evalf(Sd(gu[2],t))):

print(`The prob. of getting at least .5+r*i, for i from 1 to`, trunc(.5/r)):
print(`are: `):
print( evalf([seq([.5+r*i, PrExc(gu[2],t,.5+r*i)],i=1..trunc(.5/r))])):

od:

od:


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

end:



#SipurCG(N,m,r): All the info about comparative gambling with
#various goals
#for positions with i+j<=m living up to N tosses
#with resolution r for prob.
#For example, try:
#SipurCG(200,3,.1);
SipurCG:=proc(N,m,r) local k,h,t0,i1,g:
t0:=time():

print(`This is a comparative study of two gambling methodologies`):
print(`for playing the Chow-Robbins Coin-Tossing game,`):
print(`where you toss a fair coin, and you are allowed to stop any time`):
print(`collecting, as payoff, the average number of heads.`):
print(``):
print(`The first way is the usual one, maximizing  the expected gain.`):
print(`The second way is that you have to match or exceed a certain`):
print(` goal, and wishes to maximize the prob. of making it.`):
print(`(e.g. if you need a certain score to get admitted to your college`):
print(`of choice, and you need your SAT score to be at least 1500,`):
print(`or you want to buy something`):
print(`that costs a certain amount, and if you get one penny less`):
print(`, you can't.)`):
print(``):
print(`Below is the comparative study for the first`, m, `tosses `):
print(`where we can only toss the coin up to`,  N, `times `):
print(`(if we make it so far, we collect 1/2 if the number of tails is `):
print(`larger than the number of heads). `):
print(`------------------------------------------`):

print(`At the very beginning `):

for i1 from 1 to trunc(.5/r) do
g:=.5+r*i1:
 print(`if you want to exceed the goal of `, g):
ComparativeGambling(0,0,N,g):
od:


print(`---------------`):
for k from 1 to m do
print(`If the total number of tosses so far is`, k, `then things are as follows.`):

  for h from 0 to k-1 do
print(`----------------`):
print(`If the number of heads is`, h, `and the number of tails is`, k-h):

for i1 from 1 to trunc(.5/r) do
g:=.5+r*i1:
 print(`if you want to exceed the goal of `, g):
ComparativeGambling(h,k-h,N,g):
od:

od:
od:

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

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

end:






#IsBoundaryPt(pt,N,beta): Is the point pt on
#the boundary of the discrete region
#(0<=i+j<=N, i-j<=beta[i+j])
#For example try: IsBoundaryPt([2,0],3,[1,2]);
IsBoundaryPt:=proc(pt,N,beta) local i,j:

i:=pt[1]: j:=pt[2]:

if i+j=0 then
 RETURN(false):
fi:


if i+j=N then
 RETURN(true):
fi:

if i-j=beta[i+j] then
RETURN(true):
fi:

false:

end:



#IsInteriorPt(pt,N,beta): Is the point pt in
#the interior of the discrete region
#(0<=i+j<=N, i-j<=beta[i+j])
#For example try: IsInteriorPt([2,0],3,[1,2]);
IsInteriorPt:=proc(pt,N,beta) local i,j:

i:=pt[1]: j:=pt[2]:

if i+j=0 then
 RETURN(true):
fi:

if i<0 then
 RETURN(false):
fi:

if i+j=N then
 RETURN(false):
fi:

if i-j>=beta[i+j] then
RETURN(false):
fi:

if i-j<beta[i+j] then
RETURN(true):
fi:

false:

end:



#W(pt1,pt2,N,beta): The number of walks from pt1 to pt2
#staying in the region i+j<N i-j<beta[i+j] and quitting
#as soon as it hits the boundary i+j=N or i-j=beta[i+j]
#For example try: W([0,0],[1,1],3,[1,2]);
W:=proc(pt1,pt2,N,beta) local c,d,pt2a,pt2b,gu:
option remember:
if nops(beta)<>N-1 then
 ERROR(`Bad input`):
fi:

if not IsInteriorPt(pt1,N,beta) then
 ERROR(`Bad input`):
fi:


c:=pt2[1]:d:=pt2[2]:


if pt2=pt1 then
  RETURN(1):
fi:


pt2a:=[c-1,d]:pt2b:=[c,d-1]:

gu:=0:

if  IsInteriorPt(pt2a,N,beta) then
 gu:=gu+ W(pt1,pt2a,N,beta):
fi:

if  IsInteriorPt(pt2b,N,beta) then
 gu:=gu+ W(pt1,pt2b,N,beta):
fi:

gu:

end:


#Bd(N,beta): the set of boundary points of the discrete region 
#{[i,j]; 0<=i,0<=j, i+j<=N, i-j<=beta[i+j]} .
#For example, try:
#Bd(4,[1,2,3]);
Bd:=proc(N,beta) local gu,i,j,k,i1,j1:

gu:=[]:

for k from 1 to N-1 do
i:=(k+beta[k])/2:
j:=(k-beta[k])/2:

if not ( type(i, integer)  and type(j, integer) ) then
ERROR(`Bad input `):
fi:

gu:=[op(gu),[i,j]]:
od:

i1:=gu[k-1][1]:

[op(gu), seq([i1-j1,N-i1+j1],j1=0..i1)]:

end:




#PayOff(pt,N,beta): the pay-off, at point pt, with the boundary
#given by N,beta. For example, try:
#PayOff([0,0],5,[1,2,3,4]);
PayOff:=proc(pt,N,beta) local gu,g:
if nops(beta)<>N-1 then
 ERROR(`Bad input`):
fi:
gu:=Bd(N,beta):
add(W(pt,g,N,beta)*max(1/2,g[1]/(g[1]+g[2]))/2^(g[1]+g[2]-pt[1]-pt[2]),
 g in gu):
end:




#CRpay(i,j,N): Same as CR(i,j,N)
#using PayOff.
#For example, try:
#CRpay(0,0,10);
CRpay:=proc(i,j,N)  local b:
b:=BetaSeq(N,N-1):

if IsInteriorPt([i,j],N,b) then
 RETURN(PayOff([i,j],N,b)):
else
 RETURN(max(1/2,i/(i+j))):
fi:

end: