######################################################################
##SnakesLadders.txt: Save this file as  SnakesLadders.txt            #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read SnakesLadders.txt<Enter>                      #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: July 2019
 
print(`Created:  July-Aug. , 2019`):
print(` This is SnakesLadders.txt`):
print(`It is one of the Maple packages that accompany the article `):
print(` "A Multi-Computational Exploration of Some Games of Pure Chance" `):
print(`by Thotsaporn "Aek" Thanatipanonda and Doron Zeilberger`):

print(` available from Thanatipanonda's website and Zeilberger's website`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the general directed graph procedures type ezraG();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Simulating procedures type ezraSi();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the STORY procedures type ezraSt();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):


ezraG:=proc()

if args=NULL then
 print(` The General directed graph procedures are: GFDmc, GR1mc  `):
 print(``):

else
ezra(args):
fi:

end:

ezraSi:=proc()

if args=NULL then
 print(` The SIMULATION procedures are: ManyGames, ManyGames2P, OneGame, OneGame2P, OneRoll  `):
 print(``):

else
ezra(args):
fi:

end:

ezraSt:=proc()

if args=NULL then
 print(` The STORY procedures are:  Paper  `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: CCb, CL1, CL2, KefelRg, ProbAheadAppx,`):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Alpha,  GFD,  ProbAhead, PrW, TMdie, TMdieG `):

 print(` `):


elif nops([args])=1 and op(1,[args])=AllStates then
print(`AllStates(CB,st): All the states in the "Count your chicken game with board CB, starting at st. Try:`):
print(`AllStates(CCb1(),[0,0]);`):


elif nops([args])=1 and op(1,[args])=Alpha then
print(`Alpha(f,x,N): Given a probability generating function`):
print(`f(x) (a polynomial, rational function and possibly beyond)`):
print(`returns a list, of length N, whose `):
print(`(i) First entry is the average`):
print(` (ii): Second entry is the variance `):
print(` for i=3...N, the i-th entry is the so-called alpha-coefficients `):
print(` that is the i-th moment about the mean divided by the `):
print(` variance to the power i/2 (For example, i=3 is the skewness `):
print(` and i=4 is the Kurtosis) `):
print(`If f(1) is not 1, than it first normalizes it by dividing`):
print(`by f(1) (if it is not 0) .`):
print(` For example, try: `):
print(` Alpha(((1+x)/2)^100,x,4); `):

elif nops([args])=1 and op(1,[args])=CCb then
print(`CCb(C,D,P,S,T): The "Count Your Chickens board". It is a triple , the labels , followed by the board, followed by the set of blue squares`):
print(`An emptry square is indicated by 0. It is assumed that the last square (nops(gu[2])+1) has all of them. Try:`):
print(`CCb(C,D,P,S,T);`):

elif nops([args])=1 and op(1,[args])=CL1 then
print(`CL1(): The pair of sets [Chutes, Ladders] for the Milton-Bradley game "Chutes and Ladders:". Try:`):
print(`CL1();`):

elif nops([args])=1 and op(1,[args])=CL2 then
print(`CL2(): The pair of sets [Snakes, Ladders] for the Tradition  game "Snakes and Ladders". Try:`):
print(`CL2();`):

elif nops([args])=1 and op(1,[args])=FindNext then
print(`FindNext(CB,i,j): inputs a "Count your Chicken board, and a location i, and an outcome of the spinner  from 1 to 6, outputs`):
print(` the next location of the chicken. Try: `):
print(` FindNext(CCb(C,D,P,S,T),4,4); `):

elif nops([args])=1 and op(1,[args])=GFD then
print(`GFD(P,t): Given a list of length n-1, say, whose i-th entry is a list of the form [destination, probability], and n is the absorbing state.`):
print(`Outputs`):
print(`the list of prob. generating functions, whose i-th entry is the prob. generating function`):
print(` for the duration of the game if starting at the i-th state. Try:`):
print(` GFD(TMdie([[1,1/2],[2,1/2]],30),t); `):
print(` GFD(TMdieG([[1,1/2],[2,1/2]],7,[[[1,3]],[[4,2]]] ),t); `):
print(` GFD(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),t); `):


elif nops([args])=1 and op(1,[args])=GFDmc then
print(`GFDmc(M,t): inputs a Markov chain M with absorbing states with nops(M) vertices, outputs the`):
print(`corresponding prob. generating functions for duration, starting at that vertex. Try:`):
print(`GFDmc(GR1mc(6,1/2),t);`):

elif nops([args])=1 and op(1,[args])=GR1mc then
print(`GR1mc(N,p): the Markov chain with our data structure of the graph with vertices 1, ..., N+1 with absrobing states 1 and N+1`):
print(`with probabilities p of going right. Try:`):
print(`GR1mc(5,1/3);`):

elif nops([args])=1 and op(1,[args])=KefelRg then
print(`KefelRg(R1,R2,t): The Hadamard product of the rational functions R1 and R2. Try:`):
print(`KefelRg(t^3/(1-t-t^2),t^4/(1-t-t^3),t);`):

elif nops([args])=1 and op(1,[args])=ManyGames then
print(`ManyGames(M,i,K,k) `):
print(` Given a Snakes and Ladder game with our format where the starting state is 1 the single absorbing state is nops(M)+1 `):
print(`and M[i] is a list of pairs [desination, prob], simulates K game, returning the  list of length k consisting of`):
print(`the average, variance, and scaled moments up to the k-th. It also returns the theoretical values, for comparison, Try:`):
print(`ManyGames(TMdie([[1,1/2],[2,1/2]],100),1,1000,4); `):
print(` ManyGames(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1,1000,4); `):

elif nops([args])=1 and op(1,[args])=ManyGames2P then
print(`ManyGames2P(M,i,j,K):  simulates K games on a finite "snakes and ladders" game M where the first player is at position i`):
print(`and the second player is at position j, outputs the ratio of times 1 won`):
print(`Try: `):
print(` ManyGames2P(TMdie([[1,1/2],[2,1/2]],30),1,1,1000); `):
print(`ManyGames2P(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1,1,1000);`):

elif nops([args])=1 and op(1,[args])=OneGame then
print(`OneGame(M,i): Given a Snakes and Ladder game with our format where the starting state is 1 the single absorbing state is nops(M)+1`):
print(`and M[i] is a list of pairs [desination, prob], simulates one game, returning the number of moves needed. Try:`):
print(` OneGame(TMdie([[1,1/2],[2,1/2]],30),1); `):
print(`OneGame(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1); `):

elif nops([args])=1 and op(1,[args])=OneGame2P then
print(`OneGame2P(M,i,j):  simulates a game on a finite "snakes and ladders" game M where the first player is at position i`):
print(`and the second player is at position j, outputs the winner (1 for first, 2 for second)`):
print(`Try: `):
print(` OneGame2P(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1,1);`):

elif nops([args])=1 and op(1,[args])=OneRoll then
print(`OneRoll(L): inputs a list of positive rational mumbers that add up to 1 and outputs`):
print(` a random outcome from 1 to nops(L) according to this prob. dist.`):
print(`Try: `):
print(`OneRoll([2/3,1/4,1/12]);`):

elif nops([args])=1 and op(1,[args])=ProbAhead then
print(`ProbAhead(f,t): Given two racers who are racing towards a goal, each with the same prob. gen. function,f, in the variable, t, .`):
print(` (f is a rational function) `):
print(` where the coeff. of t^i in f is the prob. of reaching the destination in i. `):
print(` The prob. of the first racer winning (if he gets to go first), i.e. the prob. of the first getting there>=the second getting there `):
print(`Try: `):
print(`ProbAhead((t+t^2)/2, t);`):

elif nops([args])=1 and op(1,[args])=PrW then
print(`PrW(P,i,j,K): Given a stochastic game (lik Snakes and Ladders) and two locations i and j, where`):
print(`i is the location of the player whose next turn it is, and j is the location of the other player`):
print(`the probaility of the player whose turn it is to win in <=K moves. Try:`):
print(` PrW(TMdie([[1,1/2],[2,1/2]],30),1,1,1000); `):
print(` PrW(TMdieG([[1,1/2],[2,1/2]],20,[[[1,3]],[[4,2]]] ),1,1,1000); `):
print(`PrW(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1,1,1000);`):

elif nops([args])=1 and op(1,[args])=ProbAheadAppx then
print(`ProbAheadAppx(f,g,t,K): Given two racers who are racing towards a goal, with prob. gen. functions, in the variable, t, f and g respecitvely.`):
print(` (f and g are polynomials or rational functions) `):
print(` where the coeff. of t^i in f is the prob. of reaching the destination in i steps, and analogously for g. `):
print(` The prob. of the first racer winning (if he gets to go first), i.e. the prob. of the first getting there>=the second getting there `):
print(` conditioned that the first racer got there is <=K steps `):
print(`Try: `):
print(`ProbAheadAppx((t+t^2)/2, (t+t^2)/2, t,1000);`):

elif nops([args])=1 and op(1,[args])=Paper then
print(`Paper(L,n,t,C,K1,K2,K3): Inputs `):
print(` (i) a (possibly loaded) die L given as a list of pairs [outcome,prob.of coutcome],`):
print(` (ii) a positive integer n,`):
print(` (iii) a pair of Lists C=[Chutes, Ladders] where each of them indicate the list`):
print(`of chutes (snakes) and ladders, and outputs a story about the prob. generating function for the`):
print(`(iv) positive inetegers K1,K2, K3`):
print(`Outputs an article about it giving the average and variance of the duration, followed by the scaled moments up to the K1.`):
print(`It also compares the answer to a simulation with K3 games`):
print(` followed by the approximated probability of the first player winning if there are two players who `):
print(` take turns, conditioned that they last up to K2 moves). Try: `):
print(` Paper([[1,1/2],[2,1/2]],30,t,[[],[]],4,1000,1000); `):
print(` Paper([[1,1/3],[2,2/3]],50,t,[[[6,3]],[[8,15]]],6,1000,3000); `):
print(`Paper([seq([i,1/6],i=1..6)],100,t,CL1(),6,2000,1000);`):

elif nops([args])=1 and op(1,[args])=TMdie then
print(`TMdie(L,n1): Given a loaded die L, [[a1,p1],[a2,p2],..., [ak,pk]], with k outcomes a1,..., ak, with their respective probabilities`):
print(`p1, ... pk (that must add-up to 1), as well as a positive integer n1,`):
print(`forms the transition matrix with states 1,..., n1, where n1 (and beyond is the absorbing state).`):
print(`Try: `):
print(`TMdie([[1,1/2],[2,1/2]],5); `):


elif nops([args])=1 and op(1,[args])=TMdieG then
print(`TMdieG(L,n1,C): Given a loaded die L, [[a1,p1],[a2,p2],..., [ak,pk]], with k outcomes a1,..., ak, with their respective probabilities`):
print(`p1, ... pk (that must add-up to 1), as well as a positive integer n1,`):
print(` forms the transition matrix with states 1,..., n1, where n1 (and beyond is the absorbing state). `):
print(` Try: `):
print(` TMdieG([[1,1/2],[2,1/2]],7,[[[1,3]],[[4,2]]] ); `):
print(` TMdieG([seq([i,1/6],i=1..6)],100,CL1()); `):

else
print(`There is no ezra for`,args):
fi:
 
end:

##########From Cfinite

#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:

eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:




#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

for d from 1 to trunc(nops(L)/2)-2 do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:


#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):
L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(f):
fi:
end:

#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:

#Kefel(S1,S2): inputs two C-finite sequences S1,S2, and  outputs
#their product, For example, try:
#Kefel([[1,2],[3,4]],[[2,5],[-1,6]]);
Kefel:=proc(S1,S2) local d1,d2,L1,L2,L,i,K:

if nops(S1[1])<>nops(S1[2]) or nops(S2[1])<>nops(S2[2]) then
print(`Bad input`):
RETURN(FAIL):
fi:

d1:=nops(S1[1]): d2:=nops(S2[1]):
K:=2*(d1*d2)+4:

L1:=SeqFromRec(S1,K):
L2:=SeqFromRec(S2,K):

L:=[seq(L1[i]*L2[i],i=1..K)]:

GuessRec(L):

end:

#RtoC(R,t): Given a rational function R(t) finds the
#C-finite description of its sequence of coefficients
#For example, try:
#RtoC(1/(1-t-t^2),t);
RtoC:=proc(R,t) local S1,S2,D1,R1,d,f1,i,a:
R1:=normal(R):
D1:=denom(R1):
d:=degree(D1,t):
if degree(numer(R1),t)>=d then
print(`The degree of the top must be less than the degree of the bottom`):
 RETURN(FAIL):
fi:

a:=coeff(D1,t,0):
S2:=[seq(-coeff(D1,t,i)/a,i=1..degree(D1,t))]:
f1:=taylor(R,t=0,d+1):
S1:=[seq(coeff(f1,t,i),i=0..d-1)]:
[S1,S2]:
end:

#KefelR(R1,R2,t): inputs rational functions of t,
#R1 and R2,
#and outputs the rational function 
#whose coeffs. is their Hadamard product, for example try:
#KefelR(1/(1-t-t^2),1/(1-t-t^3),t);
KefelR:=proc(R1,R2,t) local S1,S2,S:

if R1=0 or R2=0 then
  RETURN(0):
fi:

S1:=RtoC(R1,t):
S2:=RtoC(R2,t):
S:=Kefel(S1,S2):
normal(CtoR(S,t)):
end:

#BUr(R,t): the break-up of R into polynomial part and rational function
BUr:=proc(R,t) local R1,T,B,P1:
R1:=normal(R):
T:=numer(R1):
B:=denom(R1):
P1:=quo(T,B,t):
R1:=normal(R1-P1):
[P1,R1]:
end:


#KefelRg(R1,R2,t): The Hadamard product of the rational functions R1 and R2
KefelRg:=proc(R1,R2,t) local gu1,gu2, P1, R1r,P2,R2r,mu,i,muT,R1T,R2T:

gu1:=BUr(R1,t):
gu2:=BUr(R2,t):

P1:=gu1[1]:
R1r:=gu1[2]:
P2:=gu2[1]:
R2r:=gu2[2]:

mu:=add(coeff(P1,t,i)*coeff(P2,t,i)*t^i,i=0..min(degree(P1,t),degree(P2,t))):

mu:=mu+add(coeff(taylor(R1r,t=0,degree(P2,t)+1),t,i)*coeff(P2,t,i)*t^i,i=0..degree(P2,t)):
mu:=mu+add(coeff(taylor(R2r,t=0,degree(P1,t)+1),t,i)*coeff(P1,t,i)*t^i,i=0..degree(P1,t)):
mu:=mu+KefelR(R1r,R2r,t):

muT:=taylor(mu,t=0,51):

R1T:=taylor(R1,t=0,51):

R2T:=taylor(R2,t=0,51):

if {seq(coeff(muT,t,i)-coeff(R1T,t,i)*coeff(R2T,t,i),i=1..50)}<>{0} then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:

mu:


end:

##########End From Cfinite


##########From HISTABRUT

#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:

#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:

##########End From HISTABRUT



#CL1(): The pair of sets [Chutes, Ladders] for the Milton-Bradley game "Chutes and Ladders. Try:
#CL1();
CL1:=proc():
[
[[16,6],[47,26],[49,11],[56,53],[62,19],[64,60],[87,24],[93,73],[95,75],[98,78]],
[[1,38], [4,14],[9,31],[21,42],[28,84],[36,44],[51,67],[71,91],[80,100]]
];
end:

#CL2(): The pair of sets [Snakes, Ladders] for the Traditions game "Snakes and Ladders. Try:
#CL2();
CL2:=proc():
[
[[16,5],[50,8],[63,20], [57,25],[98,46],[98,76],[95,90]],
[[2,44], [6,13], [9,31],[28,84],[59,61],[67,93],[70,73],[79,100]]
];
end:







#Tsamek1(L): Given a  of list of destinations with their probability [j,Pr[i->j]]
#ony mentions those destination with non-zero prob. Try:
#Tsamek1([[2,1/2],[3,1/2],[4,0],[5,0]]));
Tsamek1:=proc(L) local gu, i:
gu:=[]:

for i from 1 to nops(L) do
 if L[i][2]<>0 then
  gu:=[op(gu), L[i]]:
 fi:
od:
gu:
end:



#Tsamek(L): Given a list of lists L where L[i] is the list of destinations with their probability [j,Pr[i->j]]
#shrinks it as to only mention  destination with non-zero probabilities. Try:
#Tsamek([[[2,1/3],[3,2/3],[4,0]]]);
Tsamek:=proc(L) local i:
[seq(Tsamek1(L[i]),i=1..nops(L))]:
end:

#TMdie(L,n1): Given a loaded die L, [[a1,p1],[a2,p2],..., [ak,pk]], with k outcomes a1,..., ak, with their respective probabilities
#p1, ... pk (that must add-up to 1), as well as a positive integer n1,
#forms the transition matrix with states 1,..., n1, where n1 (and beyond is the absorbing state).
#Try:
#TMdie([[1,1/2],[2,1/2]],5);
TMdie:=proc(L,n1) local i,j,T:



for i from 1 to n1-1 do
 for j from 1 to n1 do
  T[i,j]:=0:
 od:
od:

for i from 1 to n1-1 do
 for j from 1 to nops(L) do
  if i+L[j][1]<n1 then
     T[i,i+L[j][1]]:=T[i,i+L[j][1]]+L[j][2]:
  else
  T[i,n1]:=T[i,n1]+L[j][2]:
  fi:

 od:
od:



Tsamek([seq([seq([j,T[i,j]],j=1..n1)],i=1..n1-1)]):

end:


#TMdieG(L,n1,C): Given a loaded die L, [[a1,p1],[a2,p2],..., [ak,pk]], with k outcomes a1,..., ak, with their respective probabilities
#p1, ... pk (that must add-up to 1), as well as a positive integer n1,
#forms the transition matrix with states 1,..., n1, where n1 (and beyond is the absorbing state).
#Try:
#TMdieG([[1,1/2],[2,1/2]],7,[[[1,3]],[[4,2]]] );
TMdieG:=proc(L,n1,C) local gu,K,C1,C2,j,i,mu,gu1,j1,mu1:
gu:=TMdie(L,n1):

C1:=C[1]:
C2:=C[2]:

for i from 1 to n1 do
 K[i]:=i:
od:

for i from 1 to nops(C1) do
 K[C1[i][1]]:=C1[i][2]:
od:


for i from 1 to nops(C2) do
 K[C2[i][1]]:=C2[i][2]:
od:

mu:=[]:

for i from 1 to nops(gu) do
 gu1:=gu[i]:
 mu1:=[]:
 for j from  1 to nops(gu1) do
  mu1:=[op(mu1),[subs({seq(j1=K[j1],j1=1..n1)},gu1[j][1]),gu1[j][2]]]:
 od:
mu:=[op(mu),mu1]:
od:

Tsamek(mu):

end:




#GFD(P,t): Given a list of length n-1, say, whose i-th entry is a list of the form [destination, probability], and n is the absorbing state.
#Outputs
#the list of prob. generating functions, whose i-th entry is the prob. generating function
# for the duration of the game if starting at the i-th state. Try:
#GFD(TMdie([[1,1/2],[2,1/2]],30),t); 

GFD:=proc(P,t) local i,n,f,eq,var,j,i1:

option remember:

if not type(P,list) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

n:=nops(P)+1:

var:={seq(f[i],i=1..n)}:
eq:={seq(f[i]=t*add(P[i][j][2]*f[P[i][j][1]],j=1..nops(P[i])),i=1..n-1),f[n]=1}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

[seq(factor(normal(subs(var,f[i1]))),i1=1..n-1)] :

end:




#ProbAheadAppx(f,g,t,K): Given two racers who are racing towards a goal, with prob. gen. functions, in the variable, t, f and g respecitvely.
#(f and g are polynomials or rational functions)
#where the coeff. of t^i in f is the prob. of reaching the destination in i steps, and analogously for g.
#The prob. of the first racer winning (if he gets to go first), i.e. the prob. of the first getting there>=the second getting there
#conditioned that the first racer got there is <=K steps
#Try:
#ProbAheadAppx((t+t^2)/2, (t+t^2)/2, t,1000);
ProbAheadAppx:=proc(f,g,t,K) local i,j,f1,g1:
f1:=taylor(f,t=0,K+1):
g1:=taylor(g,t=0,K+1):
add(coeff(f1,t,i)*add(coeff(g1,t,j),j=i..K),i=0..K):
end:





##start simulation procedures

#OneRoll(L): inputs a list of positive rational mumbers that add up to 1 and outputs
# a random outcome from 1 to nops(L) according to this prob. dist.
#Try
#OneRoll([2/3,1/4,1/12]);
OneRoll:=proc(L) local i,M,L1,j,ra:

if convert(L,`+`)<>1 then
 RETURN(FAIL):
fi:

if min(op(L))<0 then
 RETURN(FAIL):
fi:

M:=lcm(seq(denom(L[i]),i=1..nops(L))):

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

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

M:=L1[nops(L1)]:
ra:=rand(1..M)():

for j from 1 to nops(L1) while ra>L1[j] do od:
j:

end:

#OneGame(M,i): 
#Given a Snakes and Ladder game with our format where the starting state is 1 the single absorbing state is nops(M)+1
#and M[i] is a list of pairs [desination, prob], simulates one game, returning the number of moves needed. Try:
#OneGame(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1);

OneGame:=proc(M,i) local n,j,co,lu,L,i1:
n:=nops(M)+1:
j:=i:
co:=0:

while j<n do

 co:=co+1:

 lu:=M[j]:

 L:=[seq(lu[i1][2],i1=1..nops(lu))]:
 j:=lu[OneRoll(L)][1]:
od:


co:

end:

#ManyGames(M,i,K,k)
#Given a Snakes and Ladder game with our format where the starting state is 1 the single absorbing state is nops(M)+1
#and M[i] is a list of pairs [desination, prob], simulates K game, returning the  list of length k consisting of
#the average, variance, and scaled moments up to the k-th. Try:
#ManyGames(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1,1000,4); 

ManyGames:=proc(M,i,K,k) local t,gu,j:

gu:=add(t^OneGame(M,i),j=1..K):

evalf(Alpha(gu,t,k)), evalf(Alpha(GFD(M,t)[1],t,k)):

end:


#ProbAhead(f,t): Given two racers who are racing towards a goal, each with the same prob. gen. function,f, in the variable, t, 
#(f is a rational function) 
# where the coeff. of t^i in f is the prob. of reaching the destination in i. 
# The prob. of the first racer winning (if he gets to go first), i.e. the prob. of the first getting there>=the second getting there 
#Try: 
#ProbAhead((t+t^2)/2, t,1000);
ProbAhead:=proc(f,t) local gu:

gu:=KefelRg(f,f,t):
(subs(t=1,gu)+1)/2:
end:


#GR1mc(N,p): the Markov chain with our data structure of the graph with vertices 1, ..., N+1 with absrobing states 1 and N+1
#with probabilities p of going right. Try:
#GR1mc(5,1/3);
GR1mc:=proc(N,p) local i,T:
T[1]:={}:
T[N+1]:={}:

for i from 2 to N do
 T[i]:={[i+1,p],[i-1,1-p]}:
od:

[seq(T[i],i=1..N+1)]:

end:


#GFDmc(M,t): inputs a Markov chain M with absorbing states with nops(M) vertices, outputs the
#corresponding prob. generating functions for duration, starting at that vertex. Try:
#GFDmc(GR1mc(6,1/2),t);
GFDmc:=proc(M,t) local eq,var,X,i,j:
var:={seq(X[i],i=1..nops(M))}:

eq:={}:

for i from 1 to nops(M) do
 if M[i]={} then
  eq:=eq union {X[i]=1}:
 else
  eq:=eq union {X[i]=t*add( X[M[i][j][1]]*M[i][j][2],j=1..nops(M[i]))}:
 fi:
od:

var:=solve(eq,var):

[seq(normal(subs(var,X[i])),i=1..nops(M))]:

end:




#PrW(P,i,j,K): Given a stochastic game (lik Snakes and Ladders) and two locations i and j, where
#i is the location of the player whose next turn it is, and j is the location of the other player
#the probaility of the player whose turn it is to win in <=K moves. Try:
#PrW(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1,1,300);
PrW:=proc(P,i,j,K) local gu,t:
gu:=GFD(P,t):
if not (1<=i and i<=nops(P) and 1<=j and j<=nops(P)) then
 print(`bad input`):
 RETURN(FAIL):
fi:

evalf(ProbAheadAppx(gu[i],gu[j],t,K)):
end:


#OneGame2P(M,i,j):  simulates a game on a finite "snakes and ladders" game M where the first player is at position i
#and the second player is at position j, outputs the winner (1 for first, 2 for second)
#Try
#OneGame2P(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1,1);
OneGame2P:=proc(M,i,j) local n,i1,j1,lu,L,r:
n:=nops(M)+1:
i1:=i:
j1:=j:

while i1<n and j1<n do
 lu:=M[i1]:

 L:=[seq(lu[r][2],r=1..nops(lu))]:
 i1:=lu[OneRoll(L)][1]:

if i1>=n then
  RETURN(1):
fi:

 lu:=M[j1]:

 L:=[seq(lu[r][2],r=1..nops(lu))]:


 j1:=lu[OneRoll(L)][1]:

if j1>=n then
  RETURN(2):
fi:

od:

FAIL:

end:


#ManyGames2P(M,i,j,K):  simulates K games on a finite "snakes and ladders" game M where the first player is at position i
#and the second player is at position j, outputs the ratio of times 1 won
#Try
#ManyGames2P(TMdieG([seq([i,1/6],i=1..6)],100,CL1()),1,1,1000);
ManyGames2P:=proc(M,i,j,K) local co,r:

co:=0:

for r from 1 to K do
 if OneGame2P(M,i,j)=1 then
   co:=co+1:
 fi:
od:
co/K:
end:

####end simulation procedures

#Paper(L,n,t,C,K1,K2,K3): print inputs a (possibly loaded) die L given as a list of pairs [outcome,prob.of coutcome],
#a positive integer n, and a pair of Lists C=[Chutes, Ladders] where each of them indicate the list
#of chutes (snakes) and ladders, and outputs a story about the prob. generating function for the
#duration of a solitaire game, followed by the expectation duration, up to the K1-th moment about the mean
#followed by the approximated probability of the first player winning if there are two players who
#take turns, conditioned that they last up to K2 moves). 
#It also compares the answer to a simulation with K3 games
#Try:
#Paper(TMdie([[1,1/2],[2,1/2]],10),10,t,[[],[]],4,1000);
#Paper(TMdieG([seq([i,1/6],i=1..6)],100,t,CL1(),6,2000);
Paper:=proc(L,n,t,C,K1,K2,K3) local gu,i,ka,ka1,lu,t0,lu1:

t0:=time():

if {seq(type(L[i][2],numeric),i=1..nops(L))}={true} then

if min(seq(L[i][2],i=1..nops(L)))<0 then
 print(L, ` is not a prob. distribution`):
 RETURN(FAIL):
fi:

fi:

if add(L[i][2],i=1..nops(L))<>1 then
  print(L, ` is not a prob. distribution`):
 RETURN(FAIL):
fi:

if not (type(C,list) and nops(C)=2 and type(C[1],list) and type(C[2],list)) then
 print(C, `is illegal `):
 RETURN(FAIL):
fi:


if C<>[[],[]] then
ka:={seq(op(C[1][i]),i=1..nops(C[1])), seq(op(C[2][i]),i=1..nops(C[2]))}:

if {seq(type(ka1,integer),ka1 in ka)}<>{true} then
 print(C, `is illegal `):
 RETURN(FAIL):
fi:

if min(op(ka))<1 or max(op(ka))>n then
 print(C, `is illegal `):
 RETURN(FAIL):
fi:
fi:




print(`Statistical Analysis of a certain Chutes-and-Ladders Game `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Consider the following game with a 1-dimensional board with`, n, `squares numbered from  1 to`, n):
print(`You start at square 1 and proceed by rolling a die (or spinner) with the following possible outcomes. `):
print(``):

for i from 1 to nops(L) do
 print(`It lands on`, L[i][1], `with probability `, L[i][2]):
od:

print(`At each turn, you roll the die (or spin the spinner), and advance forward the number of squares indicated.`):
print(`Your goal is to reach the last square, square`, n, ` it does not matter if in the last round you went beyond it, either way the game ended.`):
print(``):

if C<>[[],[]] then

print(`However there are a certain number of Chutes (snakes) and Ladders`):
print(``):
print(`If your new position is at the bottom of a ladder, you immediately move to the top, and if it landed`):
print(`at the top of the chute (snake), your slide down to its bottom`):

print(``):
print(`The  Chutes (Snakes) are as follows`):
print(``):

for i from 1 to nops(C[1]) do
 print(`From `, C[1][i][1], ` down to `,  C[1][i][2]):
od:

print(``):
print(`The  Ladders are as follows`):
print(``):

for i from 1 to nops(C[2]) do
 print(`From `, C[2][i][1], ` up to `,  C[2][i][2]):
od:

fi:

gu:=GFD(TMdieG(L,n,C),t)[1]:

print(`The probabilty generating function for the random variable, "number of turns to get to the end"  is`):
print(``):
print(gu):
print(``):

print(`and in Maple notation it is`):
print(``):
lprint(gu):
print(``):

lu:=evalf(Alpha(gu,t,K1)):

print(`The expected duration of the (solitaire) game is`):
print(``):
print(evalf(lu[1])):
print(``):

print(`The variance of the duration of the (solitaire) game is`):
print(``):
print(evalf(lu[2])):
print(``):

for i from 3 to K1 do
print(`The scaled `, i, `-th moment about the mean is`):
print(``):
print(evalf(lu[i])):
print(``):
od:

print(`Summarizing the, expectation, variance and the scaled moments up to the`, K1, `-th , are `):
print(``):
print(lu):
print(``):
print(`Let's compare it to simulating `, K3, `random games , the corresponding list is`):
print(``):
lu1:=ManyGames(TMdieG(L,n,C),1,K3,K1)[1]:
print(lu1):
print(``):
print(`[Of course this changes (hopefully only slightly, by the law of large numbers)  each time] `):
print(``):

print(``):
print(`So far the game was Solitaire. Suppose that there are two players, and they keep taken turns.`):
print(`Assuming that the game did not last beyond`, K2, ` moves ` ):
print(``):
print(`The probability, that the first player is going to win is`):

print( evalf(ProbAheadAppx(gu,gu,t,K2))):
print(``):
print(`Just for fun, let's compare it to the result of simulating`, K3, `such games, and see the fraction of times the first player won`):
print(``):
print(evalf(ManyGames2P(TMdieG(L,n,C),1,1,K3)));

print(`-------------------------`):
print(``):
print(`This ends this article, that took `, time()-t0, ` seconds to generate.`):

end: