######################################################################
##GenPileGames.txt: Save this file as  GenPileGames.txt              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read GenPileGames.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 , 2019`):
print(` This is GenPileGames.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(`and also available from 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 Checking procedures type ezraCheck();, 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 Simulation procedures type ezraSi();, 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):

_MAXORDER:=20;

ezraCheck:=proc()

if args=NULL then
 print(` The checking procedures are:  Checkf1mk, CheckPGF1k, CheckPGFk1, SOS1kRopCheck  `):
 print(``):

else
ezra(args):
fi:

end:

ezraCata:=proc()

if args=NULL then
 print(` The Catalan procedures are: AlphaCata, Cata, fCata,  MSnCata `):
 print(``):

else
ezra(args):
fi:

end:

ezraSi:=proc()

if args=NULL then
 print(` The Simulations procedures are:  OneGame, OneRoll, ManyGames, ManyGamesCata, ManyGames1k `):
 print(``):

else
ezra(args):
fi:

end:

ezraSt:=proc()

if args=NULL then
 print(` The STORY procedures are:  Paper1k, PaperCata, Paperk1, PaperN, PaperSOS11Rfair, PaperSOS1kR, PaperSOS1kRS   `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  Alphak1N, Avk1N , Av21S, Av21SS, Av31S, Av31SS, EQk1, ExF2m1`):
 print(` fkm1,  fkm1Pre, GF1k, tDkm1, tDmfn , MSn1k, MyPashet, NSgf, PGF1k, PGFk1, PW , SeqFromRec, SOS11fairF, SOS1k, SOS1kRop `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:  Alpha1k, Alphak1,   NS ,  PrAheadN, SOS11Rfair, SOS1kR, SOS1kRS `):

 print(` `):



elif nops([args])=1 and op(1,[args])=Alpha1k then
print(`Alpha1k(p,k,n,N):  The average, variance ,.., followed the scaled moemnts `):
print(`to the N-th moment for the r.v. time until you have n dollars for the first time if prob. of winning a dollar is p >k/(k+1) `):
print(`and losing k dollars is 1-p,  for SYMBOLIC p and n . Note that k can be numeric (a positive integer) or symbolic. Try:`):
print(`Alpha1k(p,k,n,4);`):


elif nops([args])=1 and op(1,[args])=AlphaCata then
print(`AlphaCata(p,n,N):  The average, variance ,.., followed the scaled moemnts `):
print(`to the N-th moment for the r.v. time until you have n dollars for the first time if prob. of winning a dollar is p >1/2 for symbolic p and n. Try:`):
print(`AlphaCata(p,n,4); `):


elif nops([args])=1 and op(1,[args])=Alphak1 then
print(`Alphak1(k,p,R): inputs a pos. integer k and a symbol  p assumed to be less than k/(k+1) and ouptuts the `):
print(`list of length R, whose first entry is the expectation, followed by the variance, followed by the scaled`):
print(`moments up to R of the r.v. duration until getting a positive capital for the first time. Note that k MUST be numeric (a positive integer) `):
print(` where Pr(-1)=p Pr(k)=1-p. Try:`):
print(`Alphak1(2,p,4);`):

elif nops([args])=1 and op(1,[args])=Alphak1N then
print(`Alphak1N(k,p,R): Like Alphak1(k,p,R) but p must be symbolic.  Try:`):
print(`Alphak1N(2,1/2,4);`):

elif nops([args])=1 and op(1,[args])=Av21S then
print(`Av21S(p): inputs a symbol p assumed less than 2/3 and ouptuts the expected duration in terms of p`):
print(`of a gambling game with Prob. of losing a dollar p, and winning 2 dollars 1-p. It ends the first time the`):
print(`player has positive capital. `):
print(`It is done ab initio. Because Maple is not very good at simplifying radicals, it gives something more`):
print(`complicated than necessary. Try:`):
print(`Av21S(p);`):

elif nops([args])=1 and op(1,[args])=Av21SS then
print(`Av21SS(p): A post-processing (using our own simplification procedure) of Av21S(p) (q.v. ). Try: `):
print(`Av21SS(p);`):

elif nops([args])=1 and op(1,[args])=Av31S then
print(`Av31S(p): inputs a symbol p assumed less than 3/4 and ouptuts the expected duration in terms of p`):
print(`of a gambling game with Prob. of losing a dollar p, and winning 3 dollars 1-p. It ends the first time the`):
print(`player has positive capital. `):
print(`It is done ab initio. Try:`):
print(`Av31S(p);`):

elif nops([args])=1 and op(1,[args])=Av31SS then
print(`Av31SS(p): the pre-computed expression given by A31S(p) (q.v.). JUST for the record. Try: `):
print(`Av31SS(p);`):

elif nops([args])=1 and op(1,[args])=Avk1N then
print(`Avk1N(k,p): inputs a pos. integer k and a number p less than k/(k+1) and ouptuts the expected duration of getting above `):
print(`the x-axis Pr(-1)=p Pr(k)=1-p. `):
print(`It is a special case of Momsl1N(k,p,R) that also gives the variance and the scaled moments up to the R-th`):
print(`Try:`):
print(`Avk1N(2,1/2);`):

elif nops([args])=1 and op(1,[args])=Cata then
print(`Cata(x): the generating function for the Catalan sequence. Try:`):
print(` Cata(x); `):



elif nops([args])=1 and op(1,[args])=Checkf1mk then
print(`Checkf1mk(p): Checks f1mk(f,k,p,t): for k=1 against the Catalan case.`):
print(` Try: `):
print(`Checkg1mk(p); `):


elif nops([args])=1 and op(1,[args])=CheckPGF1k then
print(`CheckPGF1k(k,N): Checks PGF1k(k,p,t,N). Try:`):
print(`CheckGF1k(3,10);`):

elif nops([args])=1 and op(1,[args])=CheckPGFk1 then
print(`CheckPGFk1(k,N): Checks PGFk1(k,p,t,N). Try:`):
print(`CheckGFk1(3,10);`):

elif nops([args])=1 and op(1,[args])=EQk1 then
print(`EQk1(k,p,t,g): all the option for [value of root of g(1),value of g'(1)] for (-1,k) walks. Try:`):
print(`EQk1(2,p,t,g);`):

elif nops([args])=1 and op(1,[args])=ExF2m1 then
print(`ExF2m1(n): The exact formula (obtained by semi-rigorous experimental mathematics) for the expected number of`):
print(`moves in a "one step forward two steps backwards random walk" until reaching a location >=n for  the first time.`):
print(`Try:`):
print(`ExF2m1(10);`):


elif nops([args])=1 and op(1,[args])=f1mk then
print(`f1mk(f,k,p,t): the derivative of the p.g.f, in terms of f=f(t), for the duration until you have 1 Dollar with Prob(Winning a dollar)=p`):
print(`Pr(losing -k dollars)=1-p, in terms of f=f(t), k, t, and p. It is assumed that p>k/(k+1). DONE DIRECTLY. Try:`):
print(`f1mk(f,k,p,t);`):

elif nops([args])=1 and op(1,[args])=f1mkPre then
print(`f1mkPre(f,k,p,t): the derivative of the p.g.f, in terms of f=f(t), for the duration until you have 1 Dollar with Prob(Winning a dollar)=p`):
print(`Pr(losing -k dollars)=1-p, in terms of f=f(t), k, t, and p. It is assumed that p>k/(k+1). Try:`):
print(`f1mk(f,k,p,t);`):

elif nops([args])=1 and op(1,[args])=fkm1 then
print(`fkm1(f,k,t): the derivative of the p.g.f, in terms of f=f(t), for the duration until you have k Dollar with Prob(Winning k dollars)=1-p`):
print(`Pr(losing 1 dollar)=p, in terms of f=f(t), k, t, and p. It is assumed that p<1/(k+1). PRE-COMPUTED. Try:`):
print(`fkm1(f,k,t);`):

elif nops([args])=1 and op(1,[args])=fkm1Pre then
print(`fkm1Pre(f,k,t): the derivative of the p.g.f, in terms of f=f(t), for the duration until you have k Dollar with Prob(Winning k dollars)=1-p`):
print(`Pr(losing 1 dollar)=p, in terms of f=f(t), k, t, and p. It is assumed that p<1/(k+1).FROM FIRST PRINCIPLE. Try:`):
print(`fkm1Pre(f,k,t);`):

elif nops([args])=1 and op(1,[args])=fCata then
print(`fCata(p,t): the prob. generating function for the duration of the game.Try:`):
print(` fCata(p,t); `):


elif nops([args])=1 and op(1,[args])=GF1k then
print(`GF1k(k0,p,t,K,m0): the first K terms of the prob. generating function for waks with Pr(1)=p Pr(-k)=1-p until reaching m0 dollars`):
print(`Using the formula obtianed from Lagrange inversion`):
print(`Try: `):
print(`GF1k(1,p,t,10,1);`):

elif nops([args])=1 and op(1,[args])=MSn1k then
print(`MSn1k(p,k,n,N): `):
print(`The average, variance ,.., up to the N-th moment (about the mean)`):
print(`for the r.v. time until you have n dollars for the first time if prob. of winning a dollar is p >k/(k+1)`):
print(`and losing k dollars is 1-p,`):
print(`for symbolic p and n and k. Try`):
print(`MSn1k(p,k,n,4);`):

elif nops([args])=1 and op(1,[args])=ManyGames then
print(`ManyGames(G,n1,MaxRounds,N,k): `):
print(`Given a prob. distrubution of steps and a positive integer n1, simulates N games with max. number of rounds MaxRounds`):
print(` with the goal of getting to >=n1 . It returns the average, variance, up to the k-th scaled moment, Also returns the ratio of successes.Try:  `):
print(`ManyGames([[1,2/3],[-1,1/3] ],1,100,1000,4);`):

elif nops([args])=1 and op(1,[args])=ManyGames1k then
print(`ManyGames1k(p,k,n1,K1,K2,m), When p>1/2,  ManyGames([[1,p],[-1,1-p]],n1,K1,K2,m) followed by the theoretical value.`):
print(`Try: `):
print(` ManyGames1k(3/4,1,1,100,1000,4); `):

elif nops([args])=1 and op(1,[args])=ManyGamesCata then
print(`ManyGamesCata(p,n1,K1,K2,m), When p>1/2,  ManyGames([[1,p],[-1,1-p]],n1,K1,K2,m) followed by the theoretical value.`):
print(`When p<1/2, then it also returns the porb. winning a dollar in the simulation vs. the theoretical value of p/(1-p) `):
print(`Try: `):
print(` ManyGamesCata(3/4,1,100,1000,4); `):
print(` ManyGamesCata(1/4,1,1000,1000,4); `):


elif nops([args])=1 and op(1,[args])=MSnCata then
print(`MSnCata(p,N,n): The average, variance ,.., up to the N-th moment about the mean, for the r.v. time until you have n dollars`):
print(` for the first time if prob. of winning a dollar is p >1/2, and losing a dollar is 1-p. `):
print(`for symbolic p and n. Try:`):
print(`MSnCata(p,4,n); `):


elif nops([args])=1 and op(1,[args])=MyPashet then
print(`MyPashet(L): Our simplification of radicals. Try: `):
print(` MyPashet((sqrt(3)+1)/(sqrt(3)-1),3);`):


elif nops([args])=1 and op(1,[args])=NS then
print(`NS(P,K,R,n0): inputs  a probability distributions P of the form [i,p[i]] where i `):
print(` are i is an integer `):
print(` and p[i] is the prob. of getting i dollars, at any given round (all independent)`):
print(` and the game ends as soon as you have a positive amount`):
print(`computes the probability of getting to n0 dollars in <=K moves, followed by the list of length R whose first`):
print(`entry is the expected duration, the second the variance, followed by the scaled moments. Try:`):
print(`NS([[1,2/3],[-1,1/3]],50,4,1);`):

elif nops([args])=1 and op(1,[args])=NSgf then
print(`NSgf(P,K,t,n0): inputs a probability distributions P of the form {[i,p[i]]} where i is an integer `):
print(` and p[i] is the prob. of getting i dollars, and the game ends as soon as you have n0 dollars.`):
print(`Computes the first K terms of the probability generating function, in t, for the first time it gets above ground`):
print(`NSgf([[1,2/3],[-1,1/3]],50,t,1);`):

elif nops([args])=1 and op(1,[args])=OneGame then
print(`OneGame(G,n1,MaxRounds): `):
print(`Given a prob. distrubution of steps (all forward) and a positive integer n1, simulates one game with the goal of getting to >=n1 . `):
print(`the maximum number of rounds is MaxRounds. Try: `):
print(`OneGame([[1,2/3],[-1,1/3] ],1,100);`):


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])=Paper1k then
print(`Paper1k(k0,p0,n0,CutOff,R,N): A verbose version of Alpha1k(p,k,n,R) (q.v.). `):
print(`with, in the numerical checking, <=CutOff rounds,`):
print(`and in the simulation running the game  N times. Try:`):
print(`Paper1k(2,3/4,1,400,6,1000);`):

elif nops([args])=1 and op(1,[args])=PaperCata then
print(`PaperCata(p0,n0,CutOff,R,N): A paper about simple random walks stating up to the R-th moment with numerical and simulation checking. Try:`):
print(`with, in the numerical checking, <=CutOff rounds,`):
print(`and in the simulation running the game  N times. Try:`):
print(`PaperCata(2/3,3,500,4,2000);`):

elif nops([args])=1 and op(1,[args])=Paperk1 then
print(`Paperk1(k0,p0,CutOff,R,N): A verbose version of Alphak1(p,k,R) (q.v.) `):
print(`with, in the numerical checking, <=CutOff rounds,`):
print(`and in the simulation running the game  N times`):
print(`Paperk1(2,1/2,400,6,1000);`):


elif nops([args])=1 and op(1,[args])=PaperN then
print(`PaperN(P,m0,K,R): inputs P, a probability distribution with positive expecation, a positive integer m0, a large integer K`):
print(`and a positive (not too big) integer R, outputs an article about following the walk in P and`):
print(`terminating as soon as the capital is >=0 or it exceeded K rounds. It tells you`):
print(` (i) the probability of finishing in <=K rounds  `):
print(` (ii) the expectation, variance, up to the scaled R-th moment  about the mean of the random variable duration `):
print(` (iii) If there are two players who take turns walking according to this probability distribution `):
print(` the probability of the first player winning conditioned on them both terminating in <=K moves. Try: `):
print(`PaperN([[2,1/2],[-1,1/2]],1,300,4);`):

elif nops([args])=1 and op(1,[args])=PaperSOS11Rfair then
print(`PaperSOS11Rfair(m,K1): inputs a discrete variable name m and a positive integer K1, outputs an article`):
print(`about the recurrence satisfied by the sum-of-squares of the coefficients of the probability generating function`):
print(`of the random-variable, duration until reaching m dollars for the first time, if at each round you win or lose a dollar`):
print(`each with probability 1/2 `):
print(` try:`):
print(`PaperSOS11Rfair(m,1000);`):

elif nops([args])=1 and op(1,[args])=PaperSOS1kR then
print(`PaperSOS1kR(k0,p,m,f): verbose form of SOS1kR(k0,p,m,f). Try:`):
print(`PaperSOS1kR(1,1/2,m,f);`):
print(`PaperSOS1kR(1,2/3,m,f);`):
print(`PaperSOS1kR(1,3/4,m,f);`):
print(`PaperSOS1kR(2,2/3,m,f);`):
print(`PaperSOS1kR(2,3/4,m,f);`):
print(`PaperSOS1kR(3,3/4,m,f);`):
print(`PaperSOS1kR(3,4/5,m,f);`):

elif nops([args])=1 and op(1,[args])=PaperSOS1kRS then
print(`PaperSOS1kRS(k0,p,m,f): verbose form of SOS1kRS(k0,p,m,f). Try:`):
print(`PaperSOS1kRS(1,p,m,f);`):
print(`PaperSOS1kRS(2,p,m,f);`):
print(`PaperSOS1kRS(3,4/5,m,f);`):
print(`PaperSOS1kRS(3,p,m,f);`):



elif nops([args])=1 and op(1,[args])=PGF1k then
print(`PGF1k(k,p,t,N): The begining of the formal power series, in t, that is the probability generating function of the duration of `):
print(`a "gambler's ruin game with unlimited credit"`):
print(`where P(1)=p, P(-k)=1-p where you win a dollar with prob. p and lose k dollars with probability 1-p and quit as soon as you have one dollar`):
print(`in your pocket, through term  term (k+1)*N+1. Try:`):
print(`PGF1k(2,p,t,10);`):

elif nops([args])=1 and op(1,[args])=PGFk1 then
print(`PGFk1(k,p,t,N): The begining of the formal power series, in t, that is the probability generating function of the duration of `):
print(`a "gambler's ruin game with unlimited credit"`):
print(`where P(-1)=p, P(k)=1-p where you win a dollar with prob. p and lose k dollars with probability 1-p and quit as soon as you have positive capital`):
print(`in your pocket, through term  term (k+1)*N+k. Try:`):
print(`PGFk1(3,p,t,10);`):

elif nops([args])=1 and op(1,[args])=PrAheadN then
print(`PrAheadN(P,m0,K): inputs a probability distrubtion P with positive expectation and a positive integer m0`):
print(`computes the probability of the first player winning if two players follow P and the first to come to m0`):
print(`or above is declared the winner. Try:`):
print(`PrAheadN([[1,1/3],[4,1/3],[-2,1/3]],2,100);`):

elif nops([args])=1 and op(1,[args])=PW then
print(`PW(L,x,y): inputs a probability distibution [i,Pr(i)] where i can be positive or negative`):
print(`and the prob. of [1,i] is Pr(i)`):
print(`and positive integer x and non-negative integer y, outputs the probability of winding up`):
print(`at (x,y) without ever going in y>=0. Try:`):
print(`PW([[1,1/2],[-1,1/2]],6,0);`):


elif nops([args])=1 and op(1,[args])=SeqFromRec then
print(`SeqFromRec(Ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence Ope[1](n,N)f(n)=Ope[2](n)`):
print(`extends it to the first K values. Try:`):
print(`SeqFromRec([N-2,1],n,N,[1],10); `):



elif nops([args])=1 and op(1,[args])=SOS11Rfair then
print(`SOS11Rfair(m,M,K): the recurrence operator, followed by the right hand side of the sum of squares for the (1,-1) case with probability`):
print(`1/2, followed by exact checking up to K terms. Try:`):
print(`SOS11Rfair(m,M,30);`):

elif nops([args])=1 and op(1,[args])=SOS1kRS then
print(`SOS1kRS(k0,p,m,f): `):
print(`Inputs a positive integer k0, a SYMBOL p, a symbol m and a symbol f.`):
print(`Suppose you play a gambling game with stakes {1,-k0} with Pr(1)=p and Pr(-k0)=1-p `):
print(`where p is a SYMBOL (or number larger than k0/(k0+1) and less than 1), assumed to be larger than k0/(k0+1) and less than 1. `):
print(`Since p>k0/(k0+1), sooner or later you would reach m dollars.`):
print(`Let b(k0,n,m,p) be the probability of reaching m dollars for the first time in EXACTLY n rounds, and let`):
print(`f(m):=Sum(b(k0.n,m,p)^2,n=0..infinity). The output is a linear recurrence equation with polynomial coefficients`):
print(`satisfied by f(m). Try:`):
print(`SOS1kRS(1,p,m,f);`):
print(`SOS1kRS(2,p,m,f);`):
print(`SOS1kRS(3,4/5,m,f);`):
print(`SOS1kRS(3,p,m,f);`):


elif nops([args])=1 and op(1,[args])=SOS1k then
print(`SOS1k(k0,m0,p): the EXACT value for the`):
print(` sum of squares of the coefficient of the p.g.f. for "duration" if Pr(1)=p, Pr(-k)=1-p, using K terms`):
print(`with the goal of reaching m0. p>=k/(k+1)  Try: `):
print(`SOS1k(1,1,1/2);`):

elif nops([args])=1 and op(1,[args])=SOS11fairF then
print(`SOS11fairF(K): The value of SOS1k(1,K,1/2) using the pre-computed recurrence.  Try:`):
print(`SOS11fairF(100);`):

elif nops([args])=1 and op(1,[args])=SOS1kN then
print(`SOS1kN(k0,m0,p,K): the numerical approximation for the`):
print(` sum of squares of the coefficient of the p.g.f. for "duration" if Pr(1)=p, Pr(-k)=1-p, using K terms`):
print(`with the goal of reaching m0.  Try: `):
print(`SOS1kN(1,1,1/2,40);`):


elif nops([args])=1 and op(1,[args])=SOS1kR then
print(`SOS1kR(k0,p,m,f):  Inputs a positive integer k0, and symbols p, m, and f, outputs  a linear recurrence equation `):
print(`with polynomial coefficients for f(m), SOS1k(k0,m,p) (i.e. playing until you get m dollars if Pr(1)=p Pr(-k)=1-p. p>=k/(k+1). `):
print(`p must be >=k0/(k0+1) . `):
print(`Try: `):
print(`SOS1kR(1,1/2,m,f);`):
print(`SOS1kR(1,2/3,m,f);`):
print(`SOS1kR(1,3/4,m,f);`):
print(`SOS1kR(2,2/3,m,f);`):
print(`SOS1kR(2,3/4,m,f);`):
print(`SOS1kR(3,3/4,m,f);`):
print(`SOS1kR(3,4/5,m,f);`):

elif nops([args])=1 and op(1,[args])=SOS1kRop then
print(`SOS1kRop(k0,p,m,M):  Inputs a positive integer k0, and symbols p, m, and M, outputsA recurrence operator, in n,`):
print(` annihilating SOS1k(k0,m,p) where M is the shift operator in m. Try: `):
print(`SOS1kRop(1,1/2,m,M);`):

elif nops([args])=1 and op(1,[args])=SOS1kRopCheck then
print(`SOS1kRopCheck(k0,p,K1,K2):  checks SOS1kRop(k0,p,m,M) against the numerical approximations`):
print(`SOS1kN(k0,m0,p,K1) up to K2 terms. Try`):
print(`SOS1kRopCheck(1,1/2,100,30);`):

elif nops([args])=1 and op(1,[args])=SOS1kRopOld then
print(`SOS1kRopOld(k0,p,m,M):  Inputs a positive integer k0, and symbols p, m, and M, outputsA recurrence operator, in n,`):
print(` annihilating SOS1k(k0,m,p) where M is the shift operator in m. Try: `):
print(`SOS1kRopOld(1,1/2,m,M);`):

elif nops([args])=1 and op(1,[args])=tDkm1 then
print(`tDkm1(f,k,p,t,m): (td/dt)^m f(t), where f(t) is the prob. generating function for duration where Pr(k)=1-p Pr(-1)=p`):
print(`when it ends fhte first time the player has positive money`):
print(`and p>1/(k+1). Try:`):
print(`tDkm1(f,k,p,t,3);`):

elif nops([args])=1 and op(1,[args])=tDmfn then
print(`tDmfn(f,k,p,t,m,n): (td/dt)^m f(t)^n, where f(t) is the prob. generating function for duration until reaching 1 for the`):
print(`first time, and hence f(t)^n is the p.g.f. for reaching n for the first time`):
print(` where Pr(1)=p Pr(-k)=1-p`):
print(`and p>k/(k+1). Try:`):
print(`tDmfn(f,k,p,t,3,n);`):

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


##start from Findrec
Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:


#SeqFromRec(Ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence Ope[1](n,N)f(n)=Ope[2](n)
#extends it to the first K values
SeqFromRec:=proc(Ope,n,N,Ini,K)
local ope,yemin,gu,kha,L,n1,i:

ope:=Ope[1]:
yemin:=Ope[2]:
L:=degree(ope,N):

ope:=expand(subs(n=n-L,ope)/N^L):


yemin:=subs(n=n-L,yemin):
yemin:=yemin/coeff(ope,N,0):
ope:=expand(ope/coeff(ope,N,0)):
ope:=1-ope:

gu:=Ini:


for n1 from nops(Ini)+1 to K do

kha:=subs(n=n1,yemin)+add(subs(n=n1,coeff(ope,N,-i))*gu[nops(gu)+1-i],i=1..L):

gu:=[op(2..nops(gu),gu),kha]:
od:
gu:
end:

###End from Findrec

##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(G,n1,K): 
#Given a prob. distrubution of steps (all forward) and a positive integer n1, simulates one game with the goal of getting to >=n1 . 
#OneGame([seq([i,1/6],i=1..6)],100):
OneGame:=proc(G,n1,K) local j,co,L,i1,i:

j:=0:
co:=0:

 L:=[seq(G[i1][2],i1=1..nops(G))]:

for i from 1 to K while j<n1 do

 co:=co+1:

j:=j+G[OneRoll(L)][1]:

od:


co:

end:

##########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

#ManyGames(G,n1,K1,K2,k), plays OneGame(G,n1,K1) K2 times and returns the statistical information. 
#It also returns the ratio of games that finished in less than K1 moves, followed by 
#Try:
#ManyGames([[1,1/2],[2,1/2]],100,1000,4):

ManyGames:=proc(G,n1,K1,K2,k) local t,gu,j,lu:

gu:=0:

for j from 1 to K2 do

lu:=OneGame(G,n1,K1):

if lu<K1 then
 gu:=gu+t^lu:
fi:
od:

if gu=0 then
 RETURN(FAIL):
fi:

evalf(Alpha(gu/subs(t=1,gu),t,k)),evalf(subs(t=1,gu)/K2):


end:



Cata:=proc(x): (1-sqrt(1-4*x))/(2*x): end:

fCata:=proc(p,t): Cata(p*(1-p)*t^2)*p*t: end:

#MSnCata(p,N,n): 
#The average, variance ,.., up to the N-th moment for the r.v. time until you have n dollars for the first time if prob. of winning a dollar is p >1/2
#for symbolic p and n, p is larger than 1/2
MSnCata:=proc(p,N,n) local mu,gu,m0,i,t,lu,r,kha,av:

mu:=fCata(p,t)^n:

m0:=simplify(subs(t=1,mu),symbolic):

if m0<>1 then
 RETURN(FAIL):
fi:



gu:=[]:

for i from  1 to N do
mu:=t*diff(mu,t):
 gu:=[op(gu),simplify(subs(t=1,mu),symbolic)]:
od:

av:=gu[1]:
lu:=[av]:

for i from 2 to N do
kha:=add(binomial(i,r)*gu[i-r]*(-av)^r,r=0..i-1)+(-av)^i:
kha:=normal(kha):
lu:=[op(lu),kha]:
od:
factor(lu):

end:






#AlphaCata(p,n,N):  The average, variance ,.., followed the scaled moemnts 
#to the N-th moment for the r.v. time until you have n dollars for the first time if prob. of winning a dollar is p >1/2 for symbolic p and n
AlphaCata:=proc(p,n,N) local gu,i:
gu:=MSnCata(p,N,n):

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

end:


#ManyGamesCata(p,n1,K1,K2,m), When p>1/2,  ManyGames([[1,p],[-1,1-p]],n1,K1,K2,m) followed by the theoretical value.
#Try:
#ManyGamesCata(3/4,4,100,1000,4):
ManyGamesCata:=proc(p,n1,K1,K2,m) local G,gu:


if abs(p-1/2)<1/30 then
 print(p, `is too close to 1/2 for the simulation to be meaningful`):
  RETURN(FAIL):
fi:

if p<1/2 and n1>1 then
 print(`Not yet implemented`):
 RETURN(FAIL):
fi:


G:=[[1,p],[-1,1-p]]:
if p>1/2 then

gu:=ManyGames(G,n1,K1,K2,m):

if gu[2]<=1-10^(-Digits) then
 print(`Make `, K1, `bigger `):
 RETURN(FAIL):
fi:

RETURN(gu[1],evalf(AlphaCata(p,n1,m))):

else

gu:=ManyGames(G,n1,K1,K2,m):

RETURN(evalf([gu[1],AlphaCata(1-p,n1,m)]), evalf([gu[2],p/(1-p)]) ):
fi:

end:

#f1mkPre(f,k,p,t): the derivative of the p.g.f, in terms of f=f(t), for the duration until you have 1 Dollar with Prob(Winning a dollar)=p
#Pr(losing -k dollars)=1-p, in terms of f=f(t), k, t, and p. It is assumed that p>k/(k+1). Try:
#f1mkPre(f,k,p,t);
f1mkPre:=proc(f,k,p,t) local f1,gu:

gu:=f1-f/t-(k+1)*f1*(f-p*t)/f;
gu:=numer(gu):
factor(solve(gu,f1)):
end:



#f1mk(f,k,p,t): the derivative of the p.g.f, in terms of f=f(t), for the duration until you have 1 Dollar with Prob(Winning a dollar)=p
#Pr(losing -k dollars)=1-p, in terms of f=f(t), k, t, and p. It is assumed that p>k/(k+1). Try:
#f1mk(f,k,p,t);
f1mk:=proc(f,k,p,t):
-f^2/t/(k*f-k*p*t-p*t):
end:


#tDmfn(f,k,p,t,m,n): (td/dt)^m f(t)^n, where f(t) is the prob. generating function for duration where Pr(1)=p Pr(-k)=1-p
#and p>k/(k+1). Try:
#tDmfn(f,k,p,t,3,n);
tDmfn:=proc(f,k,p,t,m,n) local mu,gu:
option remember:

if not (type(m,integer) and m>=0) then
 print(`bad input`):
 RETURN(FAIL):
fi:


if m=0 then
 RETURN(f^n):
fi:

mu:=f1mk(f,k,p,t):

if m=1 then
 RETURN(normal(t*n*f^(n-1)*mu)):
fi:

gu:=tDmfn(f,k,p,t,m-1,n):

gu:=normal(diff(gu,t)+diff(gu,f)*mu):

gu:=normal(t*gu):

end:






#MSn1k(p,k,n,N): 
#The average, variance ,.., up to the N-th moment for the r.v. time until you have n dollars for the first time if prob. of winning a dollar is p >k/(k+1)
#and losing k dollars is 1-p,
#for symbolic p and n and k
MSn1k:=proc(p,k,n,N) local f,m,gu,r,t,av,i,lu,kha:

gu:=[seq(tDmfn(f,k,p,t,m,n),m=1..N)]:

gu:=normal(subs({t=1,f=1},gu)):

av:=gu[1]:
lu:=[av]:

for i from 2 to N do
kha:=add(binomial(i,r)*gu[i-r]*(-av)^r,r=0..i-1)+(-av)^i:
kha:=normal(kha):
lu:=[op(lu),kha]:
od:
factor(lu):

end:



#Checkf1mk(p): Checks f1mk(f,k,p,t): for k=1 against the Catalan case.
#Try:
#Checkg1mk(p);
Checkf1mk:=proc(p) local f,k,t,gu,lu,lu1,lu2,i:


gu:=f1mk(f,k,p,t):

gu:=subs(k=1,gu):

lu:=fCata(p,t):

lu1:=normal(simplify(diff(lu,t),symbolic)):

lu2:=normal(subs(f=lu,gu)):

lu:=taylor(lu1-lu2,t=0,70);

{seq(normal(coeff(lu,t,i)),i=0..60)};

end:


#Alpha1k(p,k,n,N):  The average, variance ,.., followed the scaled moemnts 
#to the N-th moment for the r.v. time until you have n dollars for the first time if prob. of winning a dollar is p >k/(k+1) 
#and losing k dollars is 1-p,  for SYMBOLIC p and n . Try:
#Alpha1k(p,k,n,N);
Alpha1k:=proc(p,k,n,N) local gu,i:

if type(p,numeric) and p<=k/(k+1) then
 print(`p must be a symbol or a number greater that`, k/(k+1)):
 RETURN(FAIL):
fi:


gu:=MSn1k(p,k,n,N):
simplify([gu[1],gu[2],seq(gu[i]/gu[2]^(i/2), i=3..N)]):

end:

#ManyGames1k(p,k,n1,K1,K2,m), When p>k/(k+1),  ManyGames([[1,p],[-k,1-p]],n1,K1,K2,m) followed by the theoretical value.
#Try:
#ManyGames1k(3/4,1,4,1000,1000,4):
ManyGames1k:=proc(p,k,n1,K1,K2,m) local G,gu:

if not p>k/(k+1) and p<=1 then
 print(p, `should be larger than 1/2 and less than 1`):
 RETURN(FAIL):
fi:

G:=[[1,p],[-k,1-p]]:

gu:=ManyGames(G,n1,K1,K2,m):

if gu[2]<=1-10^(-Digits) then
 print(`Make `, K1, `bigger `):
 RETURN(FAIL):
fi:

gu[1],evalf(Alpha1k(p,k,n1,m)):


end:



#fkm1Pre(f,k,t): the derivative of the p.g.f, in terms of f=f(t), for the duration until you have k Dollar with Prob(Winning k dollars)=1-p
#Pr(losing 1 dollar)=p, in terms of f=f(t), k, t, and p. It is assumed that p<1/(k+1).DONE DIRECTLY. Try:
#fkm1Pre(f,k,t);
fkm1Pre:=proc(f,k,t) local f1,gu:

gu:=f1-(k+1)*(f-1)/t-(k+1)*f1*(f-1)/f;
gu:=numer(gu):
factor(solve(gu,f1)):
end:

#fkm1(f,k,t): the derivative of the p.g.f, in terms of f=f(t), for the duration until you have k Dollar with Prob(Winning k dollars)=1-p
#Pr(losing 1 dollar)=p, in terms of f=f(t), k, t, and p. It is assumed that p<1/(k+1). PRE-COMPUTED. Try:
#fkm1(f,k,t);
fkm1:=proc(f,k,t):
-(k+1)*(f-1)*f/t/(k*f-k-1):
end:



#tDkm1(f,k,p,t,m): (td/dt)^m f(t), where f(t) is the prob. generating function for duration where Pr(k)=1-p Pr(-1)=p
#when it ends fhte first time the player has positive money
#and p>1/(k+1). Try:
#tDkm1(f,k,p,t,3);
tDkm1:=proc(f,k,p,t,m) local mu,gu:
option remember:

if not (type(m,integer) and m>=0) then
 print(`bad input`):
 RETURN(FAIL):
fi:


if m=0 then
 RETURN(f*p*t*(f-1)/(1-t*p*f) ):
fi:

mu:=fkm1(f,k,t):


gu:=tDkm1(f,k,p,t,m-1):

gu:=normal(diff(gu,t)+diff(gu,f)*mu):

gu:=normal(t*gu):

end:



#Av21S(p): inputs a symbol p assumed less than 2/3 and ouptuts the expected duration in terms of p
#ab initio. Because Maple is not very good at simplifying radicals, it gives something more
#complicated than necessary. Try:
#Av21S(p);
Av21S:=proc(p) local mu,mu1, t,g,g1,sol,i,lu:

mu:=(1-p)*t*g+(1-p)*p*t^2*g^2;

mu1:=subs(t=1,mu):

lu:=diff(mu,t)+diff(mu,g)*g1:
lu:=subs(t=1,lu):

sol:=EQk1(2,p,t,g)[2]:

mu:=[seq([subs(g=sol[i][1],mu1), subs({g=sol[i][1],g1=sol[i][2]},  lu ) ],i=1..nops(sol))]:

for i from 1 to nops(mu) do
  if simplify(subs(p=1/2,mu[i][1]))=1  and  evalf(subs(p=1/2,mu[i][2]))>0  then
   RETURN(simplify(mu[i][2])):
  fi:
od:
FAIL:
end:






#MyPashet(L): Our simplification of radicals. Try: MyPashet((sqrt(3)+1)/(sqrt(3)-1),3);
MyPashet:=proc(L,R) local x, gu,gu2C,i,gu1,gu2:

gu:=subs(sqrt(R)=x,L):

gu1:=expand(numer(gu)):
gu2:=expand(denom(gu)):

gu2C:=coeff(gu2,x,1)*x-coeff(gu2,x,0):

gu2:=expand(coeff(gu2,x,1)^2*R-coeff(gu2,x,0)^2):

gu1:=expand(gu1*gu2C):


gu1:=x*add(R^i*coeff(gu1,x,2*i+1),i=0..degree(gu1,x))+add(coeff(gu1,x,2*i)*R^i,i=0..degree(gu1,x)):

gu:=gu1/gu2:

gu:=factor(subs(x=sqrt(R),gu)):

if simplify(gu-L)<>0 then
L:
else
 gu:
fi:

end:


#Av21SS(p): The pre-computed, simplified expression for p<2/3, for the expected duration of a gambling game
#where the prob. of losing a dollar is p and of winning 2 dollars is 1-p, that ends as soon as the player has
#a positive amout (2 dollars or 1 dollars, as the case may be). Try:
#Av21SS(p):
Av21SS:=proc(p):
-1/2*(3*p+(-(3*p+1)*(-1+p))^(1/2)-1)/p/(-2+3*p):
end:




#EQk1(k,p,t,g): all the option for [value of root of g(1),value of g'(1)] for (-1,k) walks. Try:
#EQk1(2,p,t,g);
EQk1:=proc(k,p,t,g) local  eq,sol,lu,eq1,i:

eq:=g-1-p^k*(1-p)*t^(k+1)*g^(k+1):

lu:=[eq,normal(-diff(eq,t)/diff(eq,g))]:

eq1:=subs(t=1,eq):


sol:=[solve(eq1,g)]:

lu,[seq([sol[i],subs(g=sol[i],subs(t=1,lu[2]))], i=1..nops(sol))] :

end:

#Avk1N(k,p): inputs a pos. integer k and a number p less than k/(k+1) and ouptuts the expected duration of getting above 
#the x-axis Pr(-1)=p Pr(k)=1-p. Try:
#Avk1N(2,1/2);
Avk1N:=proc(k,p) local mu,mu1, t,g,g1,sol,i,lu:

if p>=k/(k+1) then
 print(`p must be less than`, k/(k+1)):
 RETURN(FAIL):
fi:

mu:=(1-p)*t*g*add((p*t*g)^i,i=0..k-1):

mu1:=subs(t=1,mu):

lu:=diff(mu,t)+diff(mu,g)*g1:
lu:=subs(t=1,lu):

sol:=EQk1(k,p,t,g)[2]:

mu:=[seq([subs(g=sol[i][1],mu1), subs({g=sol[i][1],g1=sol[i][2]},  lu ) ],i=1..nops(sol))]:

for i from 1 to nops(mu) do
  if simplify(mu[i][1])=1  and  evalf(mu[i][2])>0  then
   RETURN(simplify(mu[i][2])):
  fi:
od:
FAIL:
end:


#Av31S(p): inputs a symbol p assumed less than 2/3 and ouptuts the expected duration in terms of p
#ab initio. Because Maple is not very good at simplifying radicals, it gives something more
#complicated than necessary. Try:
#Av31S(p);
Av31S:=proc(p) local mu,mu1, t,g,g1,sol,i,lu:


mu:=(1-p)*t*g*add((p*t*g)^i,i=0..2):
mu1:=subs(t=1,mu):

lu:=diff(mu,t)+diff(mu,g)*g1:
lu:=subs(t=1,lu):

sol:=EQk1(3,p,t,g)[2]:

mu:=[seq([subs(g=sol[i][1],mu1), subs({g=sol[i][1],g1=sol[i][2]},  lu ) ],i=1..nops(sol))]:

for i from 1 to nops(mu) do
  if simplify(subs(p=1/2,mu[i][1]))=1  and  evalf(subs(p=1/2,mu[i][2]))>0  then
   RETURN(simplify(mu[i][2])):
  fi:
od:
FAIL:
end:



#Av31SS(p): The pre-computed expression outputted by Av31S(p).
#for p<3/4, for the expected duration of a gambling game
#where the prob. of losing a dollar is p and of winning 3 dollars is 1-p, that ends as soon as the player has
#a positive amout (3 dollars, 2 dollars or 1 dollars, as the case may be). Try:
#Av31SS(p):
Av31SS:=proc(p):
1/12*(8*2^(2/3)*p^2-16*2^(2/3)*p-20*2^(1/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)*p-4*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(2/3)+
8*2^(2/3)-7*2^(1/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)+3*2^(1/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)*3^(1/2)*(3+8*p+16*
p^2)^(1/2))*(2^(2/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(2/3)-4*2^(1/3)+8*2^(1/3)*p-4*2^(1/3)*p^2+2*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p
)^2)^(1/3)-2*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)*p)/(-80*p^4-8*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)*2^(2/3)*p^3+12*3^(1/2
)*(3+8*p+16*p^2)^(1/2)*p^3+32*p^3+2*2^(2/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)*3^(1/2)*(3+8*p+16*p^2)^(1/2)*p^2-16*2^(1/3)*((-7-20*p+3*3^(1/2
)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(2/3)*p^2-9*3^(1/2)*(3+8*p+16*p^2)^(1/2)*p^2+6*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)*2^(2/3)*p^2+141*p^2-4*2^(
2/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)*3^(1/2)*(3+8*p+16*p^2)^(1/2)*p+14*2^(1/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(2/3)*p
-58*p+2*2^(1/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(2/3)*p*3^(1/2)*(3+8*p+16*p^2)^(1/2)-18*3^(1/2)*(3+8*p+16*p^2)^(1/2)*p+12*((-7-20*p+3*3^(1/2)*(3
+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)*2^(2/3)*p-35-10*2^(2/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)+15*3^(1/2)*(3+8*p+16*p^2)^(1/2)+2*2^(1/3)*((-7
-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(2/3)+2*2^(2/3)*((-7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(1/3)*3^(1/2)*(3+8*p+16*p^2)^(1/2)-2*2^(1/3)*((-\
7-20*p+3*3^(1/2)*(3+8*p+16*p^2)^(1/2))*(-1+p)^2)^(2/3)*3^(1/2)*(3+8*p+16*p^2)^(1/2))/p:
end:



#PW(L,x,y): inputs a probability distibution [i,Pr(i)] where i can be positive or negatice
#and the prob. of [1,i] is Pr(i)
#and positive integer x and non-negative integer y, outputs the probability of winding up
#at (x,y) without ever going in y>=0. Try:
#PW([[1,1/2],[-1,1/2]],6,0);
PW:=proc(L,x,y) local i:
option remember:

if x=0 then
 if y=0 then
   RETURN(1):
 else
 RETURN(0):
 fi:
fi:

if y>0 then
  RETURN(0):
fi:

expand(add(PW(L,x-1,y-L[i][1])*L[i][2],i=1..nops(L))):
end:






#PaperCata(p0,n0,CutOff,R,N): A paper about simple random walks stating up to the R-th moment with numerical and simulation checking. Try:
#PaperCata(2/3,3,400,6,1000);
PaperCata:=proc(p0,n0,Cutoff,R,N) local p,n,t,gu,i,lu,t0:

t0:=time():

print(`Statistical Analysis  of the duration of  the number of  coin tosses until you have n dollars`):
print(`if at each round the probability of winning a dollar is p and losing a dollar is 1-p`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose  that at each round, you win a dollar with probability p and lose a dollar with probability 1-p and you quit `):
print(` as soon as you reach n dollars. If p>1/2, then, of course, sooner or later you will reach your goal, how long should it take?`):
print(``):
print(`The probability generating function is given explicity`):
print(``):
print(fCata(p,t)^n):
print(``):
print(`and in Maple format `):
print(``):
lprint(fCata(p,t)^n):
print(``):

gu:=AlphaCata(p,n,R):

print(`The expected duration is`):
print(``):
print(gu[1]):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu[1]):

print(``):

print(`The variance is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu[2]):

for i from 3 to  R do
print(``):
print(`The scaled `, i, `-th moment about the mean is `):
print(``):
print(gu[i]):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu[i]):
od:

print(``):

gu:=evalf(subs({p=p0,n=n0},gu)):

print(`Let's check it against numerical computations  for p=`, p0, `and n= `, n0):
print(`The exact value for these values is `):


print(gu):
print(``):
print(`Let's see what happened after`, Cutoff, ` rounds`):

lu:=evalf(NS([[1,p0],[-1,1-p0]],Cutoff,R,n0));

print(`The probability of ending in <=`, Cutoff , ` rounds is `, lu[2]):

lu:=lu[1]:

print(`The approximate values for the average, variance, and higher scaled-moments about the mean are`):
print(``):
print(lu):
print(``):
print(`compared to the exact values `):

print(``):
print(gu):
print(``):

print(`Now let's run a simulation with`,N,  `times (where you quit after`, Cutoff, ` rounds) `):

lu:=evalf(ManyGames([[1,p0],[-1,1-p0]],n0,Cutoff,N,R)):

print(`The fraction of games that ended in <=`, Cutoff,  ` rounds is`):
print(``):
print(lu[2]):
print(``):
print(`the statistical data for the simulation turned out to be`):
lu:=lu[1]:
print(``):
print(lu):
print(``):
print(`Of course, this changes each time, but it is not too far off.`):
print(``):
print(`recall that  the exact value is`):
print(``):
print(gu):
print(``):

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




#Paper1k(k0,p0,n0,CutOff,R,N): A verbose version of Alpha1k(p,k,n,R) (q.v.). Try:
#Paper1k(2,3/4,1,400,6,1000);
Paper1k:=proc(k0,p0,n0,Cutoff,R,N) local p,n,k,t,gu,i,lu,t0,g:

if p0<=k0/(k0+1) then
 print(p0, ` should have been more than`, k0/(k0+1)):
 RETURN(FAIL):
fi:

t0:=time():

print(`Statistical Analysis  of the  number of  coin tosses until you have n dollars`):
print(`if at each toss the probability of winning a dollar is p and of losing k dollars is `, 1-p):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose  that at each round, you win a dollar with probability p and lose`, k,  ` dollars with probability 1-p and you quit `):
print(` as soon as you reach 1 dollar. If p>k/(k+1), then, of course, sooner or later you will reach your goal. How long should it take?`):
print(``):
print(`Let g(t), be the formal power series, in t, satisfying the algebraic equation`):
print(``):
print(g(t)-1-p^k*(1-p)*t^(k+1)*g(t)^(k+1)=0):
print(``):
print(`and in Maple format `):
print(``):
lprint(g(t)-1-p^k*(1-p)*t^(k+1)*g(t)^(k+1)=0):
print(``):
print(`The probability generating function for the number of rounds until reaching n dollars is`, (p*t)^n*g(t)^n ):
print(``):
print(`By implicit differentiation, we can compute any desired derivative, and hence the expectation, variance, and higher moments. `):
print(``):
gu:=Alpha1k(p,k,n,R):

print(`The expected duration (for arbitrary p,k,n), provided p>k/(k+1), is`):
print(``):
print(gu[1]):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu[1]):

print(``):

print(`The variance is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu[2]):

for i from 3 to  R do
print(``):
print(`The scaled `, i, `-th moment about the mean is `):
print(``):
print(gu[i]):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu[i]):
od:

print(``):

gu:=evalf(subs({p=p0,n=n0,k=k0},gu)):

print(`Let's check it against numerical computations and computer simulations for p=`, p0, `and k=`, k0, `and n= `, n0):
print(`The exact value for these quantities are `):


print(gu):
print(``):
print(`Let's see what happened after`, Cutoff, ` rounds`):

lu:=evalf(NS([[1,p0],[-k0,1-p0]],Cutoff,R,n0));

print(`The probability of ending in <=`, Cutoff , ` rounds is `, lu[2]):

lu:=lu[1]:

print(`The approximate values for the average, variance, and higher scaled-moments about the mean are`):
print(``):
print(lu):
print(``):
print(`compared to the exact values `):

print(``):
print(gu):
print(``):


print(`Now let's run a simulation with`, N,  `times (where you quit after`, Cutoff, ` rounds) `):

lu:=evalf(ManyGames([[1,p0],[-k0,1-p0]],n0,Cutoff,N ,R)):

print(`The fraction of games that ended in <=`, Cutoff,  ` rounds is`):
print(``):
print(lu[2]):
print(``):
print(`the statistical data for the simulation turned out to be`):
lu:=lu[1]:
print(``):
print(lu):
print(``):
print(`Of course, this changes each time, but it is not too far off.`):
print(``):
print(`recall that  the exact value is`):
print(``):
print(gu):
print(``):

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



#NS(P,K,R,n0): inputs a probability distribution P of the form [i,p[i]] where i is an integer 
# and p[i] is the prob. of getting i dollars, and the game ends as soon as you have a n0
#computes the probability of getting there in <=K moves, followed by the list of length r whose first
#entry is the expected duration, the second the variance, followed by the scaled moments. Try:
#NS([[1,2/3],[-1,1/3]],50,4,1);
NS:=proc(P,K,R,n0) local t,gu,tot:
gu:=NSgf(P,K,t,n0):
tot:=subs(t=1,gu):

Alpha(gu/tot,t,R),tot:
end:

#NSgf(P,K,t,n0): inputs a probability distributions P  form [i,Pi] where i is an integer 
# and Pi is the prob. of getting i dollars, and the game ends as soon as you have have >=n0 dollars
#computes the first K terms of the probability generating function, in t, for the first time it gets above ground
#NSgf([[1,2/3],[-1,1/3]],50,t,1);
NSgf:=proc(P,K,t,n0) local gu,x,mu,lu,kha,i,j:
option remember:
lu:=add(P[i][2]*x^P[i][1],i=1..nops(P)):
mu:=lu:
gu:=0:
for i from 1 to K do
 kha:=add(coeff(mu,x,j)*x^j,j=n0..degree(mu,x)):
 gu:=gu+subs(x=1,kha)*t^i:
 mu:=mu-kha:
 mu:=expand(mu*lu):
od:
gu:

end:


#Alphak1N(k,p,R): inputs a pos. integer k and a number p less than k/(k+1) and ouptuts the 
#list of length R, whose first entry is the expectation, followed by the variance, followed by the scaled
#moments up to R of the r.v. duration until getting above 
#the x-axis with Pr(-1)=p Pr(k)=1-p. Try:
#Alphak1N(2,1/2,2);
Alphak1N:=proc(k,p,R) local mu,mu1, t,g,g1,sol,i,lu,tov,eqG,G1,MOM,lu1,av,kha,r:

if p>=k/(k+1) then
 print(`p must be less than`, k/(k+1)):
 RETURN(FAIL):
fi:

mu:=(1-p)*t*g*add((p*t*g)^i,i=0..k-1):
lu:=diff(mu,t)+diff(mu,g)*g1:

lu:=subs(t=1,lu):

mu1:=subs(t=1,mu):

sol:=EQk1(k,p,t,g)[2]:

mu:=[seq([subs(g=sol[i][1],mu1), subs({g=sol[i][1],g1=sol[i][2]},  lu ) ],i=1..nops(sol))]:

for i from 1 to nops(mu)   while not( simplify(mu[i][1])=1  and  evalf(mu[i][2])>0) do 
od:

tov:=i:

if i=nops(mu)+1 then
 RETURN(FAIL):
#else
# RETURN(simplify(mu[tov][2])):
fi:

lu:=(1-p)*t*g*add((p*t*g)^i,i=0..k-1):


eqG:=g-1-p^k*(1-p)*t^(k+1)*g^(k+1):

G1:=normal(-diff(eqG,t)/diff(eqG,g)):

MOM:=[]:

for i from 1 to R do
lu1:=subs({t=1, g=sol[tov][1], g1=sol[tov][2]},t*diff(lu,t)+t*diff(lu,g)*g1):
MOM:=[op(MOM),simplify(lu1)]:
lu:=normal(t*diff(lu,t)+t*diff(lu,g)*G1):
od:

av:=MOM[1]:
lu:=[av]:

for i from 2 to R do
kha:=add(binomial(i,r)*MOM[i-r]*(-av)^r,r=0..i-1)+(-av)^i:
kha:=simplify(kha):
lu:=[op(lu),kha]:
od:

lu:=[op(1..2,lu),seq(lu[i]/lu[2]^(i/2),i=3..R)]:
lu:

end:

#Alphak1(k,p,R): inputs a pos. integer k and a symbol p  assumed to be < k/(k+1) and ouptuts the 
#list of length R, whose first entry is the expectation, followed by the variance, followed by the scaled
#moments up to R of the r.v. duration until getting above for the first time (i.e. getting 1 or more)
#the x-axis with Pr(-1)=p Pr(k)=1-p. Try:
#Alphak1(2,p,2);
Alphak1:=proc(k,p,R) local mu,mu1, t,g,g1,sol,i,lu,tov,eqG,G1,MOM,lu1,av,kha,r:

if type(k,symbol) then
 print(k, `should have been a specific positive integer larger than 1`):
 RETURN(FAIL):
fi:

if not (type(k,integer) and k>1) then
 print(k, `should have been a specific positive integer larger than 1`):
 RETURN(FAIL):
fi:

if type(p,numeric) and p>=k/(k+1) then
 print(`p must be a symbol or a number less that`, k/(k+1)):
 RETURN(FAIL):
fi:

mu:=(1-p)*t*g*add((p*t*g)^i,i=0..k-1):
lu:=diff(mu,t)+diff(mu,g)*g1:

lu:=subs(t=1,lu):

mu1:=subs(t=1,mu):

sol:=EQk1(k,p,t,g)[2]:

mu:=[seq([subs(p=1/2,subs(g=sol[i][1],mu1)), subs(p=1/2,subs({g=sol[i][1],g1=sol[i][2]},  lu )) ],i=1..nops(sol))]:

for i from 1 to nops(mu)   while not( simplify(mu[i][1])=1  and  evalf(mu[i][2])>0) do 
od:

tov:=i:

if i=nops(mu)+1 then
 RETURN(FAIL):
fi:

lu:=(1-p)*t*g*add((p*t*g)^i,i=0..k-1):


eqG:=g-1-p^k*(1-p)*t^(k+1)*g^(k+1):

G1:=normal(-diff(eqG,t)/diff(eqG,g)):

MOM:=[]:

for i from 1 to R do
lu1:=subs({t=1, g=sol[tov][1], g1=sol[tov][2]},t*diff(lu,t)+t*diff(lu,g)*g1):
MOM:=[op(MOM),simplify(lu1)]:
lu:=normal(t*diff(lu,t)+t*diff(lu,g)*G1):
od:

av:=MOM[1]:
lu:=[av]:

for i from 2 to R do
kha:=add(binomial(i,r)*MOM[i-r]*(-av)^r,r=0..i-1)+(-av)^i:
kha:=simplify(kha):
lu:=[op(lu),kha]:
od:

lu:=[op(1..2,lu),seq(lu[i]/lu[2]^(i/2),i=3..R)]:


if k=2 then
 lu:=[simplify(Av21SS(p)),op(2..nops(lu),lu)]:
fi:

lu:

end:


#Paperk1(k0,p0,CutOff,R,N): A verbose version of Alphak1(k0,p,R) (q.v.) with the number of rounds cutoff CutOff and N the number of simulated games. Try:
#Paperk1(2,1/2,400,6,1000);
Paperk1:=proc(k0,p0,Cutoff,R,N) local p,t,gu,i,lu,t0,mu,g:

if p0 >= k0/(k0+1) then
 print(p0, ` should have been less than`, k0/(k0+1)):
 RETURN(FAIL):
fi:

t0:=time():

print(`Statistical Analysis  of  the number of  coin tosses until you have at least ONE  dollar`):
print(`if at each toss the probability of losing a dollar is p and of winning`, k0,  ` dollars is `, 1-p):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose  that at each round, you lose a dollar with probability p and win`, k0,  ` dollars with probability 1-p and you quit `):
print(` as soon as you reach at least 1 dollar. If p is less than `, k0/(k0+1)):
print(` then, of course, sooner or later you will reach your goal. How long should it take?`):
print(``):
print(`Let g(t), be the formal power series, in t, satisfying the algebraic equation`):
print(``):
print(g(t)-1-p^k0*(1-p)*t^(k0+1)*g(t)^(k0+1)=0):
print(``):
print(`and in Maple format `):
print(``):
lprint(g(t)-1-p^k0*(1-p)*t^(k0+1)*g(t)^(k0+1)=0):
print(``):
mu:=(1-p)*t*g*add((p*t*g)^i,i=0..k0-1):
print(`The probability generating function for the number of rounds until having a positive capital is` ):
print(``):
print(mu):
print(``):
print(` and in Maple format `):
print(``):
lprint(mu):
print(``):
print(`By implicit differentiation, we can compute any desired derivative, and hence the expectation, variance, and higher moments. `):
print(``):
gu:=Alphak1(k0,p,R):

print(`The expected duration, provided, p  is smaller than `, k0/(k0+1), `is `):
print(``):
print(gu[1]):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu[1]):

print(``):

print(`The variance is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu[2]):

for i from 3 to  R do
print(``):
print(`The scaled `, i, `-th moment about the mean is `):
print(``):
print(gu[i]):
print(``):
print(`and in Maple format `):
print(``):
lprint(gu[i]):
od:

print(``):

gu:=evalf(subs(p=p0,gu)):

print(`Let's check it against numerical computations and computer simulations for p=`, p0):
print(`The exact value for these quantities are `):


print(gu):
print(``):
print(`Let's see what happened after`, Cutoff, ` rounds`):

lu:=evalf(NS([[-1,p0],[k0,1-p0]],Cutoff,R,1));

print(`The probability of ending in <=`, Cutoff , ` rounds is `, lu[2]):

lu:=lu[1]:

print(`The approximate values for the average, variance, and higher scaled-moments about the mean are`):
print(``):
print(lu):
print(``):
print(`compared to the exact values `):

print(``):
print(gu):
print(``):


print(`Now let's run a simulation with`, N, ` times (where you quit after`, Cutoff, ` rounds) `):

lu:=evalf(ManyGames([[-1,p0],[k0,1-p0]],1,Cutoff,N,R)):

print(`The fraction of games that ended in <=`, Cutoff,  ` rounds is`):
print(``):
print(lu[2]):
print(``):
print(`the statistical data for the simulation turned out to be`):
lu:=lu[1]:
print(``):
print(lu):
print(``):
print(`Of course, this changes each time, but it is not too far off.`):
print(``):
print(`recall that  the exact value is`):
print(``):
print(gu):
print(``):

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







#GF1kDirect(k0,p,t,K,m0): the first K terms of the prob. generating function for walks with Pr(1)=p Pr(-k)=1-p until reaching m0 dollars
#done directly
#Try:
#GF1kDirect(1,p,t,10,1);
GF1kDirect:=proc(k0,p,t,K,m0):
 NSgf([[-k0,1-p],[1,p]],K,t,m0);
end:



#GF1k(k0,p,t,K,m0): the first K terms of the prob. generating function for walks with Pr(1)=p Pr(-k)=1-p until reaching m0 dollars
#Using the formula obtianed from Lagrange inversion
#Try:
#GF1k(1,p,t,10,1);
GF1k:=proc(k0,p,t,K,m0) local n:
add(
m0*((k0+1)*n+m0-1)!/(n!*(k0*n+m0)!)*(p^k0*(1-p))^n*p^m0
*t^((k0+1)*n+m0 )
,n=0..trunc((K-m0)/(k0+1))
):

end:

#SOS1kN(k0,m0,p,K): the numerical approximation for the
#sum of squares of the coefficient of the p.g.f. for "duration" if Pr(1)=p, Pr(-k)=1-p, using K terms
#with the goal of reaching m0.  Try: 
#SOS1kN(1,1,1/2,40);
SOS1kN:=proc(k0,m0,p,K) local n:
add(
(m0*((k0+1)*n+m0-1)!/(n!*(k0*n+m0)!)*(p^k0*(1-p))^n*p^m0)^2

,n=0..K):
end:

#SOS1k(k0,m0,p): the EXACT sum of squares of the coefficient of the p.g.f. for "duration" if Pr(1)=p, Pr(-k)=1-p
#with the goal of reaching m0.  Try: 
#SOS1k(1,1,1/2);
SOS1k:=proc(k0,m0,p) local n:
normal(
sum(
(m0*((k0+1)*n+m0-1)!/(n!*(k0*n+m0)!)*(p^k0*(1-p))^n*p^m0)^2
,n=0..infinity)):
end:



#SOS1kRopOld(k0,p,m,M):  Inputs a positive integer k0, and symbols p, m, and M, outputs a recurrence operator, in m,
#annihilating SOS1k(k0,m,p) where M is the shift operator in m. Try, followed by the left side. Try
#SOS1kRopOld(1,1/2,m,M);
SOS1kRopOld:=proc(k0,p,m,M) local n,gu:
gu:=SumTools[Hypergeometric][Zeilberger]((m*((k0+1)*n+m-1)!/(n!*(k0*n+m)!)*(p^k0*(1-p))^n*p^m)^2,m,n,M):
[gu[1],normal(limit(subs(p=3/4,gu[2]),m=infinity)) ]:
end:


#SOS1kR(k0,p,m,f):  Inputs a positive integer k0, and symbols p, m, and M, outputsA recurrence equation in f(n)
#annihilating SOS1k(k0,m,p) where M is the shift operator in m. Try:
#SOS1kR(1,1/2,m,f);
SOS1kR:=proc(k0,p,m,f) local n:
SumTools[Hypergeometric][ZeilbergerRecurrence]((m*((k0+1)*n+m-1)!/(n!*(k0*n+m)!)*(p^k0*(1-p))^n*p^m)^2,m,n,f,0..infinity):
end:


#SOS1kRop(k0,p,m,M):  Inputs a positive integer k0, and symbols p, m, and M, outputs a recurrence operator, in m,
#annihilating SOS1k(k0,m,p) where M is the shift operator in m. Try, followed by the left side. Try
#SOS1kRop(1,1/2,m,M);
SOS1kRop:=proc(k0,p,m,M) local ope,f,j:
ope:=SumTools[Hypergeometric][ZeilbergerRecurrence]((m*((k0+1)*n+m-1)!/(n!*(k0*n+m)!)*(p^k0*(1-p))^n*p^m)^2,m,n,f,0..infinity):
[subs({seq(f(m+j)=M^j,j=0..100)},op(1,ope)),op(2,ope)]:
end:

#SOS11Rfair(m,M,K): the recurrence operator, followed by the right hand side of the sum of squares for the (1,-1) case with probability
#1/2, followed by exact checking up to K terms. Try:
#SOS11Rfair(m,M,30);
SOS11Rfair:=proc(m,M,K) local gu,m0,mu,vu,i,Ope:
gu:=SOS1kRop(1,1/2,m,M):
mu:=[seq(SOS1k(1,m0,1/2),m0=1..K)]:

vu:=expand([seq(add(subs(m=m0,coeff(gu[1],M,i))*mu[m0+i],i=0..degree(gu[1],M))-subs(m=m0,gu[2]),m0=1..nops(mu)-degree(gu[1],M))]);

if vu<>[0$nops(vu)] then
 print(gu, `did not work out`):
 RETURN(FAIL):
else
RETURN(gu):
fi:

end:




#SOS1kRopCheck(k0,p,K1,K2):  checks SOS1kRop(k0,p,m,M) against the numerical approximations
#SOS1kN(k0,m0,p,K1) up to K2 terms. Try
#SOS1kRopCheck(1,1/2,100,30);

SOS1kRopCheck:=proc(k0,p,K1,K2) local m,M,gu,vu,mu,i,m0:

gu:=SOS1kRop(k0,p,m,M):

print(`gu is`, gu):

mu:=evalf([seq(SOS1kN(k0,m0,p,K1),m0=1..K2)]):
vu:=evalf([seq(add(subs(m=m0,coeff(gu[1],M,i))*mu[m0+i],i=0..degree(gu[1],M))-subs(m=m0,gu[2]),m0=1..nops(mu)-degree(gu[1],M))]);

vu:

end:


#PaperSOS11Rfair(m,K1): verbose form of SOS11Rfair(m,M,K). Also gives the K1-th value. Try
#PaperSOS11Rfair(m,1000);

PaperSOS11Rfair:=proc(m,K1) local gu,f,M,Ini,lu,Ope:

gu:=SOS1kR(1,1/2,m,f):
print(``):
print(`On the probability of the first player reaching m first if, starting at 0, at  each round they go one step right or left each`):
print(`with probability 1/2`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Consider  a game where players take turns tossing a fair coin and at each step either move left or right one unit`):
print(`and whoever reaches the location m first is declared the winner.`):
print(`What can you say about the probability of the first player  winning the game?`):
print(``):
print(`Theorem: The probabability of the first player winning is`):
print(``):
print( (1+f(m))/2 ):
print(``):
print(`where f(m) satisfies the  recurrence `):
print(``):
print(gu):
print(``):

print(`subject to the initial conditions `):

print(f(1)=SOS1k(1,1,1/2),f(2)=SOS1k(1,2,1/2)):

print(`and in Maple format `):
print(``):
lprint(gu):
print(``):

print(`subject to the initial conditions `):

lprint(f(1)=SOS1k(1,1,1/2),f(2)=SOS1k(1,2,1/2)):

print(``):

print(`Just for fun here is the sum-of-squares of the coefficient of the probability generating function`):
print(`for the duration until reaching`, K1, `dollars for the first time `):

Ope:=SOS1kRop(1,1/2,m,M):

Ini:=[SOS1k(1,1,1/2),SOS1k(1,2,1/2)]:

lu:=normal(SeqFromRec(Ope,m,M,Ini,K1)[-1]):

print(lu):
print(``):
print(`and in Maple format`):
print(``):
lprint(lu):

end:



#PaperSOS1kR(k0,p,m,f): verbose form of SOS1k(k0,p,m,f) Try:
#PaperSOS1kR(2,3/4,m,f);
PaperSOS1kR:=proc(k0,p,m,f) local gu:

gu:=SOS1kR(k0,p,m,f):

print(``):
print(`On the probability of the first player reaching m first if, starting at 0, at  each round they go one step to the right with probability `, p):
print(`or  `, k0 , `steps to the left with probabilty`, 1-p):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Consider  a game where players take turns tossing a fair coin and at each step either move  right one unit, with probability`, p):
print(`or move`, k0, `units left with probability `, 1-p):
print(``):
print(`and whoever reaches the location m first is declared the winner.`):
print(`What can you say about the probability of the first player of winning the game?`):
print(``):
print(`Theorem: The probabability of the first player winning is`):
print(``):
print( (1+f(m))/2 ):
print(``):
print(`where f(m) satisfies the  recurrence `):
print(``):
print(gu):
print(``):
print(`subject to the appropriate initial conditions `):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu):
print(``):
print(`subject to the appropriate initial conditions `):



end:



 



#SOS11fairF(K): The value of SOS1k(1,K,1/2) using the pre-computed recurrence.  Try:
#SOS11fairF(100);
SOS11fairF:=proc(K) local m,M:
option remember:
if K=1 then
-(-4+Pi)/Pi
elif
K=2 then
-(-16+5*Pi)/Pi
else
normal(SeqFromRec([2*m^2+3*m+M^2*(2*m^2+5*m+2)+M*(-12*m^2-24*m-10), -8/Pi],m,M, [-(-4+Pi)/Pi, -(-16+5*Pi)/Pi ],K)[-1]) :
fi:
end:







#PrAheadN(P,m0,K): inputs a probability distrubtion P with positive expectation and a positive integer m0
#computes the probability of the first player winning if two players follow P and the first to come to m0
#or above is declared the winner. Try:
#PrAheadN([[1,1/3],[4,1/3],[-2,1/3]],2,100);
PrAheadN:=proc(P,m0,K) local gu,t,i:

if not (add(P[i][1]*P[i][2],i=1..nops(P))>0) then
 print(`The expectation of`, P, `is not positive `):
 RETURN(FAIL):
fi:

gu:=NSgf(P,K,t,m0):
(1+add(coeff(gu,t,i)^2,i=0..degree(gu,t)))/2:
end:


#PaperN(P,m0,K,R): inputs P, a probability distribution with positive expecation, a positive integer m0, a large integer K
#and a positive (not too big) integer R, outputs an article about following the walk in P and
#terminating as soon as the capital is >=0 or it exceeded K rounds. It tells you
#(i) the probability of finishing in <=K rounds 
#(ii) the expectation, variance, up to the scaled R-th moment  about the mean of the random variable duration
#(iii) If there are two players who take turns walking according to this probability distribution
#the probability of the first player winning conditioned on them both terminating in <=K moves. Try
#PaperN([[2,1/2],[-1,1/2]],1,300,4);
PaperN:=proc(P,m0,K,R) local gu,lu,i,t0:

t0:=time():

if R<2 then
 print(`Make `, R, `at least 2 `):
 RETURN(FAIL):
fi:

if not (add(P[i][1]*P[i][2],i=1..nops(P))>0) then
 print(`The expectation of`, P, `is not positive `):
 RETURN(FAIL):
fi:

print(`Statistical Analysis of a certain Random Walk and its associated Game`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose that you start at the 0, and at each round, walk on the discrete line as follows`):
print(``):

for i from 1 to nops(P) do
if P[i][1]<0 then
 print(`with probability`, P[i][2], `move `, -P[i][1] , `units to the left `):
elif P[i][1]=0 then
  print(`with probability`, P[i][2], ` stay where you are `):
else
 print(`with probability`, P[i][2], `move `, P[i][1] , `units to the right `):
fi:
od:

print(`Since the expected progress in each move is positive, sooner or later you would make it to`, m0):
print(`but since life is finite, we will only do it for up to`, K, `rounds, `):
print(``):
print(`You do that until you reached, for the first time `, m0, `or to its right and then it ends`):
print(`or already moved`, K, `rounds, at which case you give up `):

   gu:=evalf(NS(P,K,R,m0)):
 print(`The probability that you achieved your goal is`, gu[2]):
 gu:=gu[1]:
 print(``):
 print(`conditioned on you finishing in <=`, K, `rounds `):
 print(``):
 print(` the expected duration until you reach`, m0, `or beyond for the first time is`):
 print(``):
 print(gu[1]):
 print(``):
 print(` the variance of the duration is `):
 print(``):
 print(gu[2]):
 print(``):
 print(`hence the standard deviation is`, sqrt(gu[2]) ):
  for i from 3 to R do
print(``):
   print(`The scaled `, i, `-th moment about the mean is`, gu[i]):
 od:

  lu:=evalf(PrAheadN(P,m0,K)):

 print(`Now let's turn it into a game, where two 1D random walkers, follow the above mentioned random walk, and take turns, independent of each other`):
 print(`The player to first reach the goal of ` , m0, `or above is declared the winner.`):
 print(`conditioned on both players finishing in <= `,K, `rounds, the probability of the first-to-move player to win that stupid game is`):
 print(``):
 print(lu):
 print(``):
 print(`--------------------------------------------`):
  print(``):
print(`This ends this article that took`, time()-t0, `seconds to generate. `):
print(``):
end:



#PGF1k(k,p,t,N): The begining of the formal power series, in t, that is the probability generating function of the duration of 
#a "gambler's ruin game with unlimited credit"
#where P(1)=p, P(-k)=1-p where you win a dollar with prob. p and lose k dollars with probability 1-p and quit as soon as you have one dollar
#in your pocket, through term  term (k+1)*N+1. Try:
#PGF1k(3,p,t,10);
PGF1k:=proc(k,p,t,N) local n:
add(1/(k*n+1)*binomial((k+1)*n,n)*(p^k*(1-p))^n*p*t^((k+1)*n+1),n=0..N):
end:

#CheckPGF1k(k,N): Checks PGF1k(k,p,t,N). Try:
#CheckGF1k(3,10);
CheckPGF1k:=proc(k,N) local p,t:
evalb(expand(PGF1k(k,p,t,N)-NSgf([[1,p],[-k,1-p]],(k+1)*N+1,t,1))=0):
end:



#PGFk1(k,p,t,N): The begining of the formal power series, in t, that is the probability generating function of the duration of 
#a "gambler's ruin game with unlimited credit"
#where P(-1)=p, P(k)=1-p where you win a dollar with prob. p and lose k dollars with probability 1-p and quit as soon as you have positive capital
#in your pocket, through term  term (k+1)*N+k. Try:
#PGFk1(3,p,t,10);
PGFk1:=proc(k,p,t,N) local n,r,gu,gu1:

gu:=0:
for n from 0 to N do

 for r from 1 to k do
  gu1:=r/((k+1)*n+r)
*binomial((k+1)*n+r,n)
*(p^k*(1-p))^n*p^(r-1)*(1-p)*t^((k+1)*n+r):
  gu:=gu+gu1:
 od:
od:
gu:
end:

#CheckPGFk1(k,N): Checks PGFk1(k,p,t,N). Try:
#CheckGFk1(3,10);
CheckPGFk1:=proc(k,N) local p,t:
evalb(expand(PGFk1(k,p,t,N)-NSgf([[-1,p],[k,1-p]],(k+1)*N+k,t,1))=0):
end:



#SOS1kRS(k0,p,m,f): 
#Inputs a positive integer k0, a SYMBOL p, a symbol m and a symbol f.
#Suppose you play a gambling game with stakes {1,-k0} with Pr(1)=p and Pr(-k0)=1-p 
#where p is a SYMBOL (or number larger than k0/(k0+1) and less than 1), assumed to be larger than k0/(k0+1) and less than 1. 
#Since p>k0/(k0+1), sooner or later you would reach m dollars.
#Let b(k0,n,m,p) be the probability of reaching m dollars for the first time in EXACTLY n rounds, and let
#f(m):=Sum(b(k0.n,m,p)^2,n=0..infinity). The output is a linear recurrence equation with polynomial coefficients
#satisfied by f(m). Try:
#SOS1kRS(1,p,m,f);
#SOS1kRS(2,p,m,f);
#SOS1kRS(3,4/5,m,f);
#SOS1kRS(3,p,m,f);

SOS1kRS:=proc(k0,p,m,f) local ope,M,ORD1,i:

if not (type(k0,integer) and k0>=1 and  type(m,symbol) and type(f,symbol)) then
 print(`bad input`):
 RETURN(FAIL):
fi:

if  type(p,numeric) and not (p>k0/(k0+1) and p<1) then
 print(`bad input`):
 RETURN(FAIL):
fi:

ope:= SumTools[Hypergeometric][Zeilberger](m^2*((k0+1)*n+m-1)!^2/n!^2/(k0*n+m)!^2*((p^k0*(1-p))^n)^2*(p^m)^2,m,n,M)[1]:

ORD1:=degree(ope,M):

ope:=expand(subs(m=m-ORD1,ope)/M^ORD1):

ope:=add(factor(coeff(ope,M,-i))*M^(-i),i=0..ORD1):

add(coeff(ope,M,-i)*f(m-i),i=0..ORD1)=0:
end:




#PaperSOS1kRS(k0,p,m,f): verbose form of SOS1kR(k0,p,m,f) (q.v.). Try:
#PaperSOS1kRS(1,p,m,f);
PaperSOS1kRS:=proc(k0,p,m,f) local gu:

gu:=SOS1kRS(k0,p,m,f):

print(``):
print(`On the probability of the first player reaching m first if, starting at 0, at  each round they go one step to the right with probability `, p):
print(`or  `, k0 , `steps to the left with probabilty`, 1-p):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Consider  a game where players take turns tossing a fair coin and at each step either move  right one unit, with probability`, p):
print(`or move`, k0, `units left with probability `, 1-p):
print(``):
print(`and whoever reaches the location m first is declared the winner.`):
print(`What can you say about the probability of the first player of winning the game?`):
print(``):
print(`Theorem: The probabability of the first player winning is`):
print(``):
print( (1+f(m))/2 ):
print(``):
print(`where f(m) satisfies the  recurrence `):
print(``):
print(gu):
print(``):
print(`subject to the appropriate initial conditions `):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu):
print(``):
print(`subject to the appropriate initial conditions `):



end:


#ExF2m1Old(n): The exact formula (obtained by semi-rigorous experimental mathematics) for the expected number of
#moves in a "one step forward two steps backwards random walk" until reaching a location >=n for  the first time.
#In terms of sqrt(5)
#Try:
#ExF2m1Old(10);
ExF2m1Old:=proc(n)
(fibonacci(n+2)-1)*sqrt(5)+(2*n+3-3*fibonacci(n+1)-fibonacci(n)):
end:

#ExF2m1(n): The exact formula (obtained by semi-rigorous experimental mathematics) for the expected number of
#moves in a "one step forward two steps backwards random walk" until reaching a location >=n for  the first time.
#In terms of the Golden ratio
#Try:
#ExF2m1(10);
ExF2m1:=proc(n) local phi:
phi:=(1+sqrt(5))/2:
2*(fibonacci(n+2)*phi-fibonacci(n+3) + 2 -phi +n):
end: