HelpOld:=proc(): print(`LD(p,x), OneGame(a,b,p,x), MC(a,b,p,x,K)`): print(`F(a,b,p,x) `): end: Help:=proc(): print(`P(a,b), N(a,b), n(x,y,K), FindGF(x,y), CheckGF(K) `): 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: #N(a,b): the probabilty of the Next Player winning #in a stupid game were you toss a fair coin #if it is Head you move 1 unit towards the end #if it is Tail you move 2 units towards the end #The Next player us a units away and the Previous players b units N:=proc(a,b) option remember: if a<=0 then RETURN(1): fi: if b<=0 then RETURN(0): fi: 1-(P(a-1,b)/2+P(a-2,b)/2): end: #P(a,b): the probabilty of the b- Player winning #if it is the turn of b-player to roll #in a stupid game were you toss a fair coin #if it is Head you move 1 unit towards the end #if it is Tail you move 2 units towards the end P:=proc(a,b) option remember: if a<=0 then RETURN(0): elif b<=0 then RETURN(1): else RETURN(1-(N(a,b-1)+N(a,b-2))/2): fi: end: #n(x,y,K): the truncted generating function of N(a,b) up to degree K n:=proc(x,y,K) local a,b: add(add(N(a,b)*x^a*y^b,a=0..K),b=0..K): end: #p(x,y,K): the truncted generating function of P(a,b) up to degree K p:=proc(x,y,K) local a,b: add(add(P(a,b)*x^a*y^b,a=0..K),b=0..K): end: #Finds the generating function A(x,y)=Sum(Sum(N(a,b)*x^a*y^b,a=0..infinity),b=0..infinity) FindGF:=proc(x,y) local Eq,var,a,b: Eq:={a+(x+x^2)/2*b-1/(1-x)/(1-y)+x+x^2/2, b+(1/2)*(y+y^2)*a-1/(1-x)/(1-y)+1+y/2}; var:=solve(Eq,{a,b}); a:=subs(var,a); end: #CheckGF(K): checks the Generating function given by FindGF up to order K. CheckGF:=proc(K) local A,a,b,A1,x,y: A:=FindGF(x,y): A:=taylor(A,x=0,K+1): for a from 1 to K do A1:=coeff(A,x,a): A1:=taylor(A1,y=0,K+1): print({seq(N(a,b)-coeff(A1,y,b),b=1..K)}): od: end: