######################################################################
#PURIM:  Save this file as PURIM                                     #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read PURIM<Enter>                                  #
##Then follow the instructions given there.                          #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Nov. 2011
 
print(`Created:  Nov. 2011 .`):
print(` This is  PURIM, `):
print(`a Maple package. It uses an Umbral-Operator approach`):
print(`to study prob. of winning, `):
print(` expected life, and probability distributions  `):
print(`for discrete gambling problems, in the style`):
print(`of Dubins-Savage, Kelly, and Breiman. `):
print(`This Maple package (along with the Maple package HIMURIM)`):
print(` accompanies the article `):
print(` How to Gamble If You're in a Hurry`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
print(`and also 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`):
 print(` is available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/PURIM  .`):


print(`----------------------------------------------------`):
 print(`For a list of the Strategy procedures type ezraStr();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`----------------------------------------------------`):

print(`----------------------------------------------------`):
 print(`For a list of the Checking procedures type ezraCh();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`----------------------------------------------------`):


 print(`For a list of the MAIN procedures type ezra();`):
print(`, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`----------------------------------------------------`):

with(combinat):


ezraCh:=proc()

if args=NULL then
 print(` The  checking procedures are `):
 print(` CheckMsteps, EDdd, PrWdd,`):


else
ezra(args):
fi:

end:

ezraStr:=proc()

if args=NULL then
 print(` The  procedures describing strategies are `):
 print(` AllS, BoS, BrS, IsValidStrat, KS, RS, TS`):

else
ezra(args):
fi:

end:


ezra:=proc()

if args=NULL then
 print(`The MAIN procedures are the following: Msteps, MstepsV, OneSt `):

elif nops([args])=1 and op(1,[args])=AllS then

print(`AllS(N): the set of all strategies with N dollar cap`):
print(`For example, try:`):
print(`AllS(6);`)

elif nops([args])=1 and op(1,[args])=BoS then
print(`BoS(N): The Bold strategy `):
print(`For example, try:`):
print(`BoS(100);`):

elif nops([args])=1 and op(1,[args])=BrS then
print(`BrS(N,f,c): The Breiman strategy with factor f`):
print(`and c For example, try:`):
print(`BrS(100,1/5,1/2);`):


elif nops([args])=1 and op(1,[args])=CheckMsteps then
print(`CheckMsteps(S,p,i,M): checks Msteps, for example, try:`):
print(`CheckMsteps(TS(5),p,3,5);`):

elif nops([args])=1 and op(1,[args])=EDdd then
print(`EDdd(S,p,i,T): Given a strategy S, (with a casino with exit capital`):
print(`nops(S)+1, and prob. of winning a single round p, and an integer i`):
print(`between 1 and N-1`):
print(`computes the expected duration if the game ends in <=M rounds`):
print(`if the current capital is i. For example, try:`):
print(`EDdd([1,1,1],1/2,1,3);`):

elif nops([args])=1 and op(1,[args])=IsValidStrat then
print(`IsValidStrat(L,N): Given a list descriving a gambling strategy,`):
print(`verifies that it is legal. For example, try:`):
print(`IsValidStrat([1,-1,3],4);`):

elif nops([args])=1 and op(1,[args])=KS then
print(`KS(N,f): The Kelly strategy with factor f`):
print(`For example, try:`):
print(`KS(100,1/5);`):


elif nops([args])=1 and op(1,[args])=Msteps then
print(`Msteps(S,p,i,t,s,M): inputs a strategy S ( a list of length N-1, say)`):
print(`a prob. p (numeric, between 0 and 1, or symbolic), an initial capital`):
print(`i (between 1 and N-1), variables t and s, and a positive integer M`):
print(`and you can gamble up to m rounds, and exit when you are either broke`):
print(`or made N dollars, or your time is up.`):
print(`The output is a list of length 6 as follows.`):
print(`(i) a loser (a polynomial in s whose coeff. of s^r, say, `):
print(` tells you the`):
print(` prob. of exiting after r rounds as a loser),`):
print(` (ii) a winner (a polynomial in s whose coeff. of s^r, say, `):
print(`tells you the`):
print(`prob. of exiting after r rounds as a winner .`):
print(`(iii) neither winner or loser, but with a capital between 1 and N-1`):
print(`a polynomial in t whose coeff. of t^k, say, tells you the prob.`):
print(`of still being in the game after m rounds, and having `):
print(`capital exactly k.`):
print(`(iv) the prob. of exiting with the whole capital`):
print(`(v) the expected duration of such a game `):
print(`(vi) the expected gain, assuming that you are allowed to keep`):
print(`the capital that you have when quitting the game after M rounds`):
print(`For example, try:`):
print(`Msteps([1,1,1,1],1/2,2,t,s,10);`):

elif nops([args])=1 and op(1,[args])=MstepsV then
print(`MstepsV(S,p,i,t,s,M): A verbose form of Msteps(S,p,i,t,s,)`):
print(`For example, try:`):
print(`MstepsV([1,1,1,1],1/2,2,t,s,10);`):


elif nops([args])=1 and op(1,[args])=OneSt then
print(`OneSt(S,p,f,t): Given a strategy S, which is `):
print(`a list of positive integers`):
print(`of length N-1 (i.e. N:=nops(S)+1), such that`):
print(`S[i] (that must be smaller than min(i,N-i))`):
print(`tells you how many dollars to stake if you currently have i dollars`):
print(`in a gambling game where the prob. of winning any single round`):
print(`is p (a number between 0 and 1, or left as a symbol),`):
print(`and a polynomial f of the variable t,`):
print(`whose meaning is that the coeff. of t^i (say)(1<i<N) tells you the`):
print(`current prob. of having i dollars, outputs a list of length 3.`):
print(`The third entry is the analogous`):
print(`polynomial after one round of gambling. The first entry is the`):
print(`prob. of exiting a loser (becoming broke). The second entry is the`):
print(`prob. of exiting a winner (getting the target amount N).`):
print(` For example, try:`):
print(`OneSt([1,1,1],1/2,t^2,t);`):

elif nops([args])=1 and op(1,[args])=PrWdd then
print(`PrWdd(S,p,i,M): Given a strategy S, (with a casino with exit capital`):
print(`nops(S)+1, and prob. of winning a single round p, and an integer i`):
print(`between 1 and N-1`):
print(`computes the probability of winning in <=M rounds`):
print(`if the current capital is i. For example, try:`):
print(`PrWdd([1,1,1],3/7,1,3);`):


elif nops([args])=1 and op(1,[args])=RS then
print(`RS(N): A random strategy with N dollars cap`):
print(`For example, try:`):
print(`RS(100);`):

elif nops([args])=1 and op(1,[args])=TS then
print(`TS(N): The timid strategy `):
print(`For example, try:`):
print(`TS(100);`):

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




#Vecs(L): all the vectors v[i] such that 1<=v[i]<=L[i]
#For example, try:
#Vecs([2,3,2]);
Vecs:=proc(L) local n,i,L1,mu,mu1:
option remember:
n:=nops(L):
if not type(L,list) then
 RETURN(FAIL):
fi:
if {seq(type(L[i],posint),i=1..n)}<>{true} then
 RETURN(FAIL):
fi:
if n=1 then
 RETURN({seq([i],i=1..L[1])}):
fi:
L1:=[op(1..n-1,L)]:
mu:=Vecs(L1):
{seq(seq([op(mu1),i],i=1..L[n]),mu1 in mu)}:
end:

#AllS(N): the set of all strategies with N dollar cap
#For example, try:
#AllS(6);
AllS:=proc(N) Vecs(BoS(N)): end:

#BoS(N): The Bold strategy 
#For example, try:
#BoS(100);
BoS:=proc(N) local i: 
[seq(min(i,N-i),i=1..N-1)]: end:

#BrS(N,f,c): The Breiman strategy with factor f
#and c
#For example, try:
#BrS(100,1/5,1/2);
BrS:=proc(N,f,c) local i:
if c>=1 then
[seq(min(N-i,i,trunc(f*i)+1),i=1..N-1)]:
else
[seq(min(trunc(f*i)+1,N-i,i),i=1..trunc(N*c)), seq(min(i,N-i),i=trunc(N*c)+1..N-1)]:
fi:
end:

#KS(N,f): The Kelly strategy with factor f
#For example, try:
#KS(100,1/5);
KS:=proc(N,f) local i:
[seq(min(trunc(f*i)+1,N-i,i),i=1..N-1)]:
end:

#RS(N): A random strategy with N dollars cap
#For example, try:
#RS(100);
RS:=proc(N) local i: 
[seq(rand(1..min(i,N-i))(),i=1..N-1)]: end:

#TS(N): The timid strategy 
#For example, try:
#TS(100);
TS:=proc(N): [1$(N-1)]: end:


#IsValidStrat(L,N): Given a list descriving a gambling strategy,
#verifies that it is legal. For example, try:
#IsValidStrat([1,-1,3],4);
IsValidStrat:=proc(L,N) local i:
if not type(N,integer) or N<1 then
RETURN(false):
fi:

if not type(L,list) then
 RETURN(false):
fi:

if nops(L)<>N-1 then
RETURN(false):
fi:
if {seq(type(L[i],integer),i=1..N-1)}<>{true} then
    RETURN(false):
fi:
if {seq(evalb(L[i]>=0),i=1..N-1)}<>{true} then
    RETURN(false):
fi:

if {seq(evalb(L[i]<=i),i=1..N-1)}<>{true} then
    RETURN(false):
fi:

if {seq(evalb(L[i]<=N-i),i=1..N-1)}<>{true} then
    RETURN(false):
fi:

true:
end:


#OneSt(S,p,f,t): Given a strategy S, which is 
#a list of positive integers
#of length N-1 (i.e. N:=nops(S)+1), such that
#S[i] (that must be smaller than min(i,N-i))
#tells you how many dollars to stake if you currently have i dollars
#in a gambling game where the prob. of winning any single round
#is p (a number between 0 and 1, or left as a symbol),
#and a polynomial f of the variable t,
#whose meaning is that the coeff. of t^i (1<i<N) tells you the
#current prob. of having i dollars, outputs a list of length 3
#whose third entry is the the analogous
#polynomial after one round of gambling, the first entry is the
#prob. of exiting a loser (becoming broke), and the second entry is
#prob. of exiting a winner (getting the target amount N).
#For example, try:
#OneSt([1,1,1],1/2,t^2,t);
OneSt:=proc(S,p,f,t) local gu,i,q,N,lu1,lu2:
N:=nops(S)+1:
if not IsValidStrat(S,N) then
print(`Bad input`, S, `not a valid strategy`):
RETURN(FAIL):
fi:
q:=1-p:
gu:=expand(add(coeff(f,t,i)*t^i*(q*t^(-S[i])+p*t^S[i]),i=1..N-1)):
lu1:=coeff(gu,t,0):
lu2:=coeff(gu,t,N):
gu:=gu-lu1-lu2*t^N:
[lu1,lu2,gu]:
end:



#Msteps(S,p,i,t,s,M): inputs a strategy S ( a list of length N-1, say)
#a prob. p (numeric, between 0 and 1, or symbolic), an initial capital
#i (between 1 and N-1), variables t and s, and a positive integer M
#and you can gamble up to m rounds, and exit when you are either broke
#or made N dollars, or your time is up, having played M games.
#The output is as follows
#(i)a loser (a polynomial in s whose coeff. of s^r, say, tells you the
#prob. of exiting after r rounds as a loser (broke)),
#(ii) a winner (a polynomial in s whose coeff. of s^r, say, tells you the
#prob. of exiting after r rounds as a winner),
#(iii) neither winner or loser, but with a capital between 1 and N-1
#a polynomial in t whose coeff. of t^k, say, tells you the prob.
#of still being in the game after m rounds, and having capital exactly k.
#(iv): the prob. of exiting with the whole capital
#(v) the expected duration of such a game 
#(vi) the expected gain, assuming that you are allowed to keep
#the capital that you have when quitting the game after M rounds
#For example, try:
#Msteps([1,1,1,1],1/2,2,t,s,10);
Msteps:=proc(S,p,i,t,s,M) local gu,lu1,lu2,j,mu1,mu2,N:
N:=nops(S)+1:
gu:=t^i:
lu1:=0:
lu2:=0:
for j from 1 to M do
 gu:=OneSt(S,p,gu,t):
 lu1:=lu1+gu[1]*s^j:
 lu2:=lu2+gu[2]*s^j:
 gu:=gu[3]:
od:

mu1:=subs(s=1,diff(lu1+lu2,s))+M*subs(t=1,gu):
mu2:=subs(s=1,lu2)*N+subs(t=1,diff(gu,t)):
[lu1,lu2,gu,subs(s=1,lu2),mu1,mu2]:
end:

###from HIMURIM (for checking only)
#PrWdd(S,p,i,M): Given a strategy S, (with a casino with exit capital
#nops(S)+1, and prob. of winning a single round p, and an integer i
#between 1 and N-1
#computes the probability of winning in <=M rounds
#if the current capital is i. For example, try:
#PrWdd([1,1,1],1/2,1,3);
PrWdd:=proc(S,p,i,M) local N:
option remember:
N:=nops(S)+1:
if M=0 then
 if i=N then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:
if i=N then
 RETURN(1):
fi:
if i=0 then
 RETURN(0):
fi:
expand((1-p)*PrWdd(S,p,i-S[i],M-1)+p*PrWdd(S,p,i+S[i],M-1)):
end:


#EDdd(S,p,i,M): Given a strategy S, (with a casino with exit capital
#nops(S)+1, and prob. of winning a single round p, and an integer i
#between 1 and N-1
#computes the expected duration if the game ends in <=M rounds
#if the current capital is i. For example, try:
#EDdd(1/2,[1,1,1],1,3);
EDdd:=proc(S,p,i,M) local N:
option remember:
N:=nops(S)+1:
if M=0 then
  RETURN(0):
fi:
if i=N then
 RETURN(0):
fi:
if i=0 then
 RETURN(0):
fi:
expand((1-p)*EDdd(S,p,i-S[i],M-1)+p*EDdd(S,p,i+S[i],M-1)+1):
end:

###end from HIMURIM (for checking only)


#CheckMsteps(S,p,i,M): checks Msteps, for example, try:
#CheckMsteps(TS(5),p,3,5);
CheckMsteps:=proc(S,p,i,M) local gu,t,s:
gu:=Msteps(S,p,i,t,s,M):
evalb(subs(s=1,gu[2])=PrWdd(S,p,i,M)),
evalb(subs(s=1,diff(gu[1]+gu[2],s))+M*subs(t=1,gu[3])=EDdd(S,p,i,M))
end:






#MstepsV(S,p,i,t,s,M): A verbose form of Msteps(S,p,i,t,s,)
MstepsV:=proc(S,p,i,t,s,M) local gu,N:
gu:=Msteps(S,p,i,t,s,M):
N:=nops(S)+1:
print(`Your Winning prospects If You Have to Leave the Casino in`, M, `rounds.`):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`Suppose that you have`,  i, `dollars, and you enter a casino where your probability of winning a round is`, p):
print(`You owe`, N, `dollars to a loan shark, who would kill you if you don't return the debt by`, M, `rounds.`):
print(`any lesser amount would do you no good.`):
print(`Suppose you follow the  strategy below,`):
print(` where the j-th entry tells you what you would stake`):
print(`if you currently have j dollars .`):
print(S):
print(`Your probability of staying alive, i.e. leaving the casino as a winner with`, N, `dollars. is`):
print(evalf(gu[4])):
print(`If you exit a loser, the polynomial whose coeff. of s^r tells you the probability of exiting after exactly r rounds is`):
print(evalf(gu[1])):


print(`If you exit a winner, the polynomial whose coeff. of s^r tells you the probability of exiting after exactly r rounds is`):
print(evalf(gu[2])):

print(`If you quit the casion after`, M, `rounds without being either a winner or loser, the polynomial whose coeff. of t^r`):
print(`tells you the exact probability of having r dollars is`):
lprint(evalf(gu[3])):

print(`The expected duration of such a game is`, evalf(gu[5])):

print(`If the loan shark unexpectedly disappers, and spares your life, then the expected gain, if you are also allowed to keep your`):
print(`current capital is`, evalf(gu[5])):
print(`The whole thing took`, time(), `seconds to compute.`):
end: