Help:=proc(): print(`LD(p,x), OneGame(a,b,p,x), MC(a,b,p,x,K)`): 
print(`F(a,b,p,x) `):
end:


#LD(p,x): Simulates one throw of a die whose prob. gen. function
#is p(x). For example LD((x+x^2+x^3+x^4+x^5+x^6)/6,x);
LD:=proc(p,x) local M,P,i,L,L2,ra,r:

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

M:=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,M*coeff(P,x,i)]]:
 fi:
od:

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

ra:=rand(1..M):
r:=ra():

for i from 1 to nops(L2) while L2[i]< r do od:

L[i][1]:

end:


#OneGame(a,b,p,x): playing SCM backgammon with a general die
#given by p(x) where the guy who is about to roll is a units
#away from 0 and the other guy is b units 
OneGame:=proc(a,b,p,x) local a1,b1,G,cu:

G:=[[a,b]]:
cu:=[a,b]:
while cu[1]>0 and cu[2]>0 do

  cu:=[cu[1]-LD(p,x),cu[2]]:

   G:=[op(G),cu]:

    if cu[1]<=0 then
       RETURN(G,1):
     fi:

 cu:=[cu[1],cu[2]-LD(p,x)]:

 if cu[2]<=0 then
    RETURN(G,2):
 fi:
G:=[op(G),cu]:
od:





end:


#MC(a,b,p,x,K): estimating the prob. of the first player to win
#in a SMC backgammon with prob. dist. p(x) doing it K times
MC:=proc(a,b,p,x,K) local c,i,jeff:
c:=0:

for i from 1 to K do
jeff:=OneGame(a,b,p,x)[2]:
 
  if jeff=1 then
     c:=c+1:
  fi:
od:

c/K:

end:


#F(a,b,p,x): the EXACT probability of winning a SMC game
#when you are a units away, and the opponet is b units
#and it is your turn from the end, and prob. dist. p(x)
F:=proc(a,b,p,x) local dennis,P:
option remember:

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

if b<=0 then
 RETURN(0):
fi:

if a<=0 then
 RETURN(1):
fi:


1- add(coeff(P,x,i)*F(b,a-i,p,x),i=ldegree(P,x)..degree(P,x)):

end: