HelpM:=proc(): print(`LD(p,x), GGR(p,x,i,N), MG(p,x,i,N,K) `): end:


Help:=proc(): print(`Probs(p,x,N) , ED(p,x,N), SymProbs(p,x,i,N) `): end:

#Given a loaded die whose gen. function is p(x)
#(prob. of having  a step of size i is the coeff. of
#x^i in p(x). For example, for a fair coin it is
#(1/2)*x+(1/2)*x^(-1). For a fair die
#(1/6)*(x+x^2+x^3+x^4+x^5+x^6)
LD:=proc(p,x) local P,N,L,i,L1,j,r,ra:


if {seq(evalb( coeff(p,x,i)>=0)  ,i=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:


if {seq(type(coeff(p,x,i),rational)  ,i=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:

P:=expand(p/subs(x=1,p));

N:=lcm( seq( denom( coeff(P,x,i) ) ,i=ldegree(P,x)..degree(P,x))):


L:=[]:

for i from ldegree(P,x) to degree(P,x) do
 if coeff(P,x,i)<>0 then
   L:=[op(L), [i, N*coeff(P,x,i)]]:
  fi:
od:

ra:=rand(1..N):

L1:=[seq( L[i][2], i=1..nops(L))]:

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

r:=ra():

for j from 1 while L1[j]<r do   od:

L[j][1]:

end:



#GGR(p,x,i,N): Inputs a discrete prob. dist. p(x)
#and a pos. integers i and N, and simulates a game
#where you enter the casino with i dollars
#and play until you have non-positive amount or >=N dollars
#The output is either or N, followed by
#the duration of the game
GGR:=proc(p,x,i,N) local m,c:

c:=0:

m:=i:

while m>0 and m<N do
 m:=m+LD(p,x):
 c:=c+1:
od:

m,c:

end:

#MG(p,x,i,N,K): Playing GGR(p,x,i,N) K times
#how many times we won?/K, followed by the average
#duration of a game
MG:=proc(p,x,i,N,K) local w, c, i1,g:
w:=0: c:=0:

for i1 from 1 to K do
g:=GGR(p,x,i,N):

if g[1]>= N then
 w:=w+1:
fi:

c:=c+g[2]:

od:

evalf(w/K),evalf(c/K):

end:


#inputs a polynomial p of x with pos. coeff. encoding a
#general die in a casino where you have to leave as soon
#as you either have <=0 or >=N dollars. Outputs a list
#of length N-1 whose i^th item is the prob. of winning
#if you enter the casino with i dollars
Probs:=proc(p,x,N) local eq, var, z,i,j,P:


if {seq(evalb( coeff(p,x,i)>=0)  ,i=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:


if {seq(type(coeff(p,x,i),rational)  ,i=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:

P:=expand(p/subs(x=1,p));


var:={ seq(z[i],i=1..N-1)}:


eq:={}:

for i from 1 to N-1 do

eq:=eq union { z[i]- add(z[i+j]*coeff(P,x,j),j=ldegree(P,x)..degree(P,x)) }:

od:

eq:=subs({seq(z[i]=0, i=ldegree(P,x)..0)},eq):

eq:=subs({seq(z[i]=1, i=N..N+degree(P,x))},eq):

var:=solve(eq,var):

[seq(subs(var,z[i]),i=1..N-1)]:

end:




#inputs a polynomial p of x with pos. coeff. encoding a
#general die in a casino where you have to leave as soon
#as you either have <=0 or >=N dollars. Outputs a list
#of length N-1 whose i^th item is the expected duration of the
#game 
#if you enter the casino with i dollars
ED:=proc(p,x,N) local eq, var, z,i,j,P:


if {seq(evalb( coeff(p,x,i)>=0)  ,i=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:


if {seq(type(coeff(p,x,i),rational)  ,i=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:

P:=expand(p/subs(x=1,p));


var:={ seq(z[i],i=1..N-1)}:


eq:={}:

for i from 1 to N-1 do

eq:=eq union { z[i]- add(z[i+j]*coeff(P,x,j),j=ldegree(P,x)..degree(P,x))=1 }:

od:

eq:=subs({seq(z[i]=0, i=ldegree(P,x)..0)},eq):

eq:=subs({seq(z[i]=0, i=N..N+degree(P,x))},eq):

var:=solve(eq,var):

[seq(subs(var,z[i]),i=1..N-1)]:

end:




#inputs a polynomial p of x with pos. coeff. encoding a
#general die in a casino where you have to leave as soon
#as you either have <=0 or >=N dollars. Outputs an explicit
#expression in i and N for the prob. of winning
SymProbs:=proc(p,x,i,N) local P, R,i1,GS,C,Q,var,eq:


if {seq(evalb( coeff(p,x,i1)>=0)  ,i1=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:


if {seq(type(coeff(p,x,i1),rational)  ,i1=ldegree(p,x)..degree(p,x))}<>
{true}  then
 RETURN(FAIL):
fi:

P:=expand(p/subs(x=1,p));

#fair case
if subs(x=1, diff(P,x))=0 then
Q:=normal((P-1)/(x-1)^2):
R:=[solve(Q,x)]:

GS:=add(C[i1]*R[i1]^i,i1=1..nops(R)) + C[nops(R)+1]+C[nops(R)+2]*i :
var:={seq(C[i1],i1=1..nops(R)+2)}: 

else
#unfair case


Q:=normal((P-1)/(x-1)):
R:=[solve(Q,x)]:

GS:=add(C[i1]*R[i1]^i,i1=1..nops(R)) + C[nops(R)+1] :
var:={seq(C[i1],i1=1..nops(R)+1)}:
fi:

eq:={seq(subs(i=j,GS)=0,j=ldegree(P,x)+1..0)}
union
{seq(subs(i=j+N,GS)=1,j=0..degree(P,x)-1)}:

var:=solve(eq,var):

subs(var,GS):

end: