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: