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: