#April 18, 2013: general gambling theory
Help:=proc(): print(` GGgame(p,s,N,L)`): 
print(`ManyGGgames(p,s,N,HowMany,L), VGW(p,N,L), VGD(p,N,L) `):
print(`IsAB(L1,L2) , AllVecsLess(L), AllStra(N), BestStrat(p,N) `):
end:

#GGgame(p,s,N,L):inputs prob. p of success (a rational number)
#starting capital, and exit capital N (if you are a winner)
#once you have 0 dollars you must leave the casino
#And a list of integers L, such that 1<=L[i]<=min(i,N-i)
#that tells you that you should gamble on L[i] dollars
#if you have i dollars
#Returns, true (false) if you won (or lost) followed
#by the number of rounds it took
GGgame:=proc(p,s,N,L) local ca,du,toss:

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

ca:=s:
du:=0:

while ca>0 and ca<N do

toss:=LC(p):

du:=du+1:

 if toss then
   ca:=ca+L[ca]:
 else
   ca:=ca-L[ca]:
 fi:
od:

evalb(ca>=N), du:


end:

#ManyGGgames(p,s,N,HowMany,L): repeats GGgame(p,s,N,L) 
#HowMany times, and returns the percentage of times
#the gambler won, and the average duration if it won,
#followed by the average duration if it lost followed
#by the average duration
ManyGGgames:=proc(p,s,N,HowMany,L) local DUw, DUl, DU, W,i,jake:
W:=0:
DUw:=0:
DUl:=0:
DU:=0:

for i from 1 to HowMany do
 jake:=GGgame(p,s,N,L):
  
  if jake[1] then
     W:=W+1:
     DUw:=DUw+jake[2]:
     DU:=DU+jake[2]:
   else
    DUl:=DUl+jake[2]:
    DU:=DU+jake[2]:
   fi:
od:

if W=0 or W=HowMany then
  print(`This is never going to happen`):
   RETURN(FAIL):
fi:

evalf(W/HowMany), evalf(DUw/W), evalf(DUl/(HowMany-W)),
evalf(DU/HowMany):

end:

#VGW(p,N,L): inputs prob. of winning one round, p
#and the exit capital N, outputs the List of
#size N-1 whose i-th entry is the EXACT prob.
#of winning if you currently posses i dollars
#following the strategy L
VGW:=proc(p,N,L) local w,eq,unk,i,sol:
if {seq(evalb(L[i]>=1 and L[i]<=min(i,N-i)),i=1..N-1)}<>{true}
   then
   RETURN(FAIL):
 fi:
#w[i]: prob. of ultimately winning  with current capital i
#in Gambler's ruin name

unk:={ seq(w[i],i=0..N)}:

eq:={w[0]=0, w[N]=1,

seq(w[i]=p*w[i+L[i]]+(1-p)*w[i-L[i]], i=1..N-1)}:

sol:=solve(eq,unk):

subs(sol,  [seq(w[i],i=1..N-1)]):


end:

#VGD(p,N,L): inputs prob. of winning one round, p
#and the exit capital N, outputs the List of
#size N-1 whose i-th entry is the EXACT EXpected Duration
#of winning if you currently posses i dollars
#and follow strategy L
VGD:=proc(p,N,L) local Du,eq,unk,i,sol:
if {seq(evalb(L[i]>=1 and L[i]<=min(i,N-i)),i=1..N-1)}<>{true}
   then
   RETURN(FAIL):
 fi:

#Du[i]: expected duration of
#a gambler's ruin with max. cap N and prob. p of winning 
#a dolllar

unk:={ seq(Du[i],i=0..N)}:

eq:={Du[0]=0, Du[N]=0,

seq(Du[i]=p*Du[i+L[i]]+(1-p)*Du[i-L[i]]+1, i=1..N-1)}:

sol:=solve(eq,unk):

subs(sol,  [seq(Du[i],i=1..N-1)]):


end:



#IsAB(L1,L2) is the list L1 always >= list L2
IsAB:=proc(L1,L2) local i:

evalb({seq(evalb(L1[i]>=L2[i]),i=1..nops(L1))}={true}):


end:


#AllVecsLess(L): The SET of all the vectors of length nops(L)
#of pos. integers such that 1<=M[i]<=L[i]

AllVecsLess:=proc(L) local n, L1, Phil,i, tom:
option remember:
n:=nops(L):

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

L1:=[op(1..n-1,L)] :

Phil:=AllVecsLess(L1):


{seq(seq([op(tom),i], tom in Phil), i=1..L[n])}:

end:

#AllStra(N): The set of all possible strategies with max capital
#N 
AllStra:=proc(N) local i:

AllVecsLess([seq(min(i,N-i),i=1..N-1)]):

end:


#BestStra(p,N)
BestStra:=proc(p,N) local champs,cu,S,cu1,i:
S:=AllStra(N):
champs:={S[1]}:
cu:=VGW(p,N,S[1]):

for i from 2 to nops(S) do

cu1:=VGW(p,N,S[i]):

if cu1=cu then
  champs:=champs union {S[i]}:
elif IsAB(cu1,cu) then
   champs:={S[i]}:
 cu:=cu1:
fi:
od:

champs:

end:















###################old stuff from C22.txt
#April 15, 2013, starting gambling theory

Help22:=proc(): print(` LC(p) , GRgame(p,s,N) `): 
print(`ManyGames(p,s,N,HowMany), VW(p,N) , VD(p,N) `):
end:


#LC(p): simulating a loaded coin whose
#that returns true with prob. p, p rational
LC:=proc(p) local a,b,ra:

a:=numer(p):
b:=denom(p):

ra:=rand(1..b)():

evalb(ra<=a):

end:

#GRgame(p,s,N):inputs prob. p of success (a rational number)
#starting capital, and exit capital N (if you are a winner)
#once you have 0 dollars you must leave the casino
#Returns, true (false) if you won (or lost) followed
#by the number of rounds it took
GRgame:=proc(p,s,N) local ca,du,toss:


ca:=s:
du:=0:

while ca>0 and ca<N do

toss:=LC(p):

du:=du+1:

 if toss then
   ca:=ca+1:
 else
   ca:=ca-1:
 fi:
od:

evalb(ca>=N), du:


end:




#ManyGames(p,s,N,HowMany): repeats GRgame(p,s,N) 
#HowMany times, and returns the percentage of times
#the gambler won, and the average duration if it won,
#followed by the average duration if it lost followed
#by the average duration
ManyGames:=proc(p,s,N,HowMany) local DUw, DUl, DU, W,i,jake:
W:=0:
DUw:=0:
DUl:=0:
DU:=0:

for i from 1 to HowMany do
 jake:=GRgame(p,s,N):
  
  if jake[1] then
     W:=W+1:
     DUw:=DUw+jake[2]:
     DU:=DU+jake[2]:
   else
    DUl:=DUl+jake[2]:
    DU:=DU+jake[2]:
   fi:
od:

if W=0 or W=HowMany then
  print(`This is never going to happen`):
   RETURN(FAIL):
fi:

evalf(W/HowMany), evalf(DUw/W), evalf(DUl/(HowMany-W)),
evalf(DU/HowMany):

end:


#VW(p,N): inputs prob. of winning one round, p
#and the exit capital N, outputs the List of
#size N-1 whose i-th entry is the EXACT prob.
#of winning if you currently posses i dollars
VW:=proc(p,N) local w,eq,unk,i,sol:

#w[i]: prob. of ultimately winning  with current capital i
#in Gambler's ruin name

unk:={ seq(w[i],i=0..N)}:

eq:={w[0]=0, w[N]=1,

seq(w[i]=p*w[i+1]+(1-p)*w[i-1], i=1..N-1)}:

sol:=solve(eq,unk):

subs(sol,  [seq(w[i],i=1..N-1)]):


end:



#VD(p,N): inputs prob. of winning one round, p
#and the exit capital N, outputs the List of
#size N-1 whose i-th entry is the EXACT EXpected Duration
#of winning if you currently posses i dollars
VD:=proc(p,N) local Du,eq,unk,i,sol:

#Du[i]: expected duration of
#a gambler's ruin with max. cap N and porb. p of winning 
#a dolllar

unk:={ seq(Du[i],i=0..N)}:

eq:={Du[0]=0, Du[N]=0,

seq(Du[i]=p*Du[i+1]+(1-p)*Du[i-1]+1, i=1..N-1)}:

sol:=solve(eq,unk):

subs(sol,  [seq(Du[i],i=1..N-1)]):


end: