######################################################################
##RUIN: Save this file as RUIN   To use it, stay in the              #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read RUIN : <Enter>                                #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers  dot edu.                              # 
#######################################################################
 
#Created: Aug.-Sept. 2006
 
print(`Created: Aug.-Sept. 2006`):
print(`This is RUIN`):
print(`to study the Gambler's Ruin Problem.`):
print(`It accompanies the paper `):
print(`Symbol-Crunching with the Gambler's Ruin Problem`):
print(`by  Doron Zeilberger`):
print(`available from the author's website`):
print():
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print():
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg`):
 print(`For a list of the procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print():

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedurs are: BVh, CheckGRlife, CheckGRmomentFair`):
 print(` CheckGRmomentLoaded, CheckGRmoment `):
 print( ` GuessPol, GRlifeBmoment, GRlifeMoment1, LoadedDie, PowerToBinom`):
 print(` RmomentSlow `):
else
ezra(args):
fi:

end:



ezra:=proc()

if args=NULL then
 print(`The main procedures are: BmomentR, BmomentRp, BmomentRw `):
 print(` ExpLife, ExpLifeW, `):
 print(` Game, GRlife, GRlifeMoment, GRlifeW, GRwin, GRwinS`):
 print(`  PrAtnF,PrAtnFad, PrAtnW, PrAtnWad, Rmoment, Rmomentp`):
 print(`RmomentW, Story, StoryW, Stat, Tourn `):
 


elif nops([args])=1 and op(1,[args])=BmomentR then
print(`BmomentR(r,n,A): the r^th Binomial moment of the r.v. duration`):
print(`of game in a (fair-coin)`):
print(`Gambler's Ruin problem in [0,A] starting at $n$`):
print(`For example, try: BmomentR(2,n,A);`):

elif nops([args])=1 and op(1,[args])=BmomentRp then
print(`BmomentRp(r,p,n,A): the r^th Binomial moment of the r.v. duration`):
print(`of game in a (loaded coin with prob. p of winning a dollar)`):
print(`Gambler's Ruin problem in [0,A] starting at $n$`):
print(`For example, try: BmomentRp(2,p,n,A);`):

elif nops([args])=1 and op(1,[args])=BmomentRw then
print(`BmomentRw(r,n,A): the r^th Binomial moment of the r.v. duration`):
print(`of game in a WINNING (fair coin) `):
print(`Gambler's Ruin problem in [0,A] starting at $n$`):
print(`For example, try: BmomentR(2,n,A);`):

elif nops([args])=1 and op(1,[args])=BVh then

print(`BVh(p,x,n,A,a,L,Left1,Right1):`):
print(`given a Laurent polynomial p(x) where the multipliicity`):
print(`of its roots are given by the list L`):
print(`Solves (symbolically) in terms of the roots`):
print(`corresponding to L labelled by the indexed variable`):
print(`a, and two lists, Left1, Right1, of length`):
print(`-ldegree(p,x) and degree(p,x) resp.`):
print(`the homog. boundary value problem`):
print(`p(N)f(n)=0 (where N is the shift in n: Ng(n):=g(n+1))`):
print(`subject to the boundary conditions`):
print(`f(-i)=Left1[i+1], i=0..nops(Left1)-1`):
print(`f(A+i)=Right1[A+i], i=0..nops(Right1)-1`):
print(`For example, try:`):
print(`BVh(1-(x+1/x)/2,x,n,A,a,[2],[0],[1]);`):


elif nops([args])=1 and op(1,[args])=CheckGRlife then
print(`CheckGRlife(P,L,i,A,K): Checks the exact value`):
print(`by simulation,  the exact value for the expectation`):
print(`of the r.v. duration of a Gambler's Ruin game`):
print(`with loaded-die P,L in [0,A] starting at i. For example, try:`):
print(`CheckGRlife([1/2,1/2],[-1,1],2,4,100);`):

elif nops([args])=1 and op(1,[args])=CheckGRmomentFair then
print(`CheckGRmomentFair(i,A,K,r): Checks `):
print(`by simulation,  the exact value for the r^th moment.`):
print(`ABOUT THE MEAN, of the r.v. duration of a Gambler's Ruin game`):
print(`with a fair coin, [0,A] starting at i. For example, try:`):
print(`CheckGRmomentFair(2,4,100,3);`):

elif nops([args])=1 and op(1,[args])=CheckGRmomentLoaded then
print(`CheckGRmomentLoaded(i,A,p, K,r): Checks `):
print(`by simulation,  the exact value for the r^th moment.`):
print(`ABOUT THE MEAN, of the r.v. duration of a Gambler's Ruin game`):
print(`with a loaded coin (Prob. of winning being p), `):
print(` in [0,A] starting at i. For example, try:`):
print(`CheckGRmomentLoaded(2,4,1/3,100,3);`):

elif nops([args])=1 and op(1,[args])=CheckGRmoment then
print(`CheckGRmoment(P,L,i,A,K,r): Checks, `):
print(`by simulation,  the exact value for the r^th moment.`):
print(`ABOUT THE MEAN, of the r.v. duration of a Gambler's Ruin game`):
print(`with loaded-die P,L in [0,A] starting at i. For example, try:`):
print(`CheckGRmoment([1/2,1/2],[-1,1],2,4,100,3);`):

elif nops([args])=1 and op(1,[args])=Game then
print(`Game(P,L,i,A): Given a list of non-negative rational`):
print(`numbers P and a list of outcomes L`):
print(`such that the prob. that it will lend at L[i] is P[i] simulates`):
print(`a game where one exists as soon as one has <=0 dollars`):
print(`or >=A dollars and starts at i`):
print(`it returns the list of partial outcomes`):
print(`For example, try: Game([1/2,1/2],[-1,1],2,4);`):

elif nops([args])=1 and op(1,[args])=GRlife then
print(`GRlife(P,L,A): Inputs a loaded die P,L and a positive integer`):
print(`A and outputs the list whose i^th item is the`):
print(`expected duration in a Gambler's Ruin game in [0,A] that`):
print(`starts at i`):
print(`For example, try: GRlife([1/2,1/2],[-1,1],4) ;`):

elif nops([args])=1 and op(1,[args])=GRlifeBmoment then
print(`GRlifeBmoment(P,L,A,r): Inputs a loaded die P,L and a`):
print(`positive integer A and outputs the list whose i^th item is the`):
print(`r-th binomial moment of the r.v. duration of game if it `):
print(`starts at i`):
print(`For example, try: GRlifeBmoment([1/2,1/2],[-1,1],4,2) ;`):

elif nops([args])=1 and op(1,[args])=GRlifeMoment then
print(`GRlifeMoment(P,L,A,r): Inputs a loaded die P,L and a`):
print(`positive integer A and outputs the list whose i^th item is the`):
print(`r-th moment ABOUT THE MEAN of the r.v. duration of game if it `):
print(`starts at i`):
print(`For example, try: GRlifeMoment([1/2,1/2],[-1,1],4,2) ;`):

elif nops([args])=1 and op(1,[args])=GRlifeMoment1 then
print(`GRlifeMoment1(P,L,A,r): Inputs a loaded die P,L and a`):
print(`positive integer A and outputs the list whose i^th item is the`):
print(`r-th moment of the r.v. duration of game if it `):
print(`starts at i`):
print(`For example, try: GRlifeMoment1([1/2,1/2],[-1,1],4,2) ;`):


elif nops([args])=1 and op(1,[args])=GRlifeW then
print(`GRlifeW(P,L,A): Inputs a loaded die P,L and a positive integer`):
print(`A and outputs the list whose i^th item is the`):
print(`(relative) expected duration in a Gambler's Ruin game in [0,A] that`):
print(`starts at i, assuming that he won`):
print(`For example, try: GRlifeW([1/2,1/2],[-1,1],4) ;`):

elif nops([args])=1 and op(1,[args])=GRwin then
print(`GRwin(P,L,A): Inputs a loaded die P,L and a positive integer`):
print(`A and outputs the list whose i^th item is the`):
print(`prob. of winning in a Gambler's Ruin game in [0,A] that`):
print(`starts at i`):
print(`For example, try: GRwin([1/2,1/2],[-1,1],4) ;`):

elif nops([args])=1 and op(1,[args])=GRwinS then
print(`GRwinS(P,L,n,A): Inputs a loaded die P,L and SYMBOLS n and`):
print(`A and outputs the expression in terms of the`):
print(`roots of a certain polynomial for the explicit expression`):
print(`for the prob. of winning in a Gambler's Ruin game in [0,A] that`):
print(`starts at i`):
print(`For example, try: GRwinS([1/2,1/2],[-1,1],n,A) ;`):

elif nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(L,n,s0): guesses a polynomial of degree d in n for`):
print(` the list L, such that L[i]=P[i] for i>=s0 for example, try: `):
print(`GuessPol([(seq(i,i=1..10),n,1);`):

elif nops([args])=1 and op(1,[args])=ExpLife then
print(`ExpLife(n,A): Given symbols n and A,`):
print(`guesses an explicit expression`):
print(`(polynomial in n and A)`):
print(`for expectation of the r.v.`):
print(`duration of a fair-coin gambler's ruin game in [0,A]`):
print(`if you start at n. For example, try:`):
print(`ExpLife(n,A):`):


elif nops([args])=1 and op(1,[args])=ExpLifeW then
print(`ExpLifeW(n,A): Given symbols n and A,`):
print(`guesses an explicit expression`):
print(`(polynomial in n, Laurent polynomial in A)`):
print(`for expectation of the r.v.`):
print(`duration of a WINNING fair-coin gambler's ruin game in [0,A]`):
print(`if you start at n. For example, try:`):
print(`ExpLifeW(n,A):`):

elif nops([args])=1 and op(1,[args])=PrAtnF then
print(`PrAtnF(P,L,n,A,t): Inputs a loaded die P,L and a positive integer`):
print(`A and  a positive integer r, and an integer n between 0 and A,`):
print(`outputs the probability that the game will end exactly after`):
print(`t tosses.`):
print(`For example, try: PrAtnF([1/2,1/2],[-1,1],2,4,2) ;`):

elif nops([args])=1 and op(1,[args])=PrAtnFad then
print(`PrAtnFad(P,L,n,A,t): Inputs a loaded die P,L and a positive integer`):
print(`A and  a positive integer r, and an integer n between 0 and A,`):
print(`outputs the probability that the game will after <=`):
print(`t tosses.`):
print(`For example, try: PrAtnFad([1/2,1/2],[-1,1],2,4,2) ;`):

elif nops([args])=1 and op(1,[args])=PrAtnW then
print(`PrAtnW(P,L,n,A,t): Inputs a loaded die P,L and a positive integer`):
print(`A and  a positive integer r, and an integer n between 0 and A,`):
print(`outputs the probability that the game will end exactly after`):
print(`t tosses.`):
print(`For example, try: PrAtnW([1/2,1/2],[-1,1],2,4,2) ;`):

elif nops([args])=1 and op(1,[args])=PrAtnWad then
print(`PrAtnWad(P,L,n,A,t): Inputs a loaded die P,L and a positive integer`):
print(`A and  a positive integer r, and an integer n between 0 and A,`):
print(`outputs the probability that the game will end exactly after`):
print(`<=t tosses.`):
print(`For example, try: PrAtnWad([1/2,1/2],[-1,1],2,4,2) ;`):

elif nops([args])=1 and op(1,[args])=Rmoment then
print(`Rmoment(n,A,r): Given symbols n and A and a positive`):
print(`integer r>=2 finds an explicit expression`):
print(`(polynomial in n, Laurent polynomial in A)`):
print(`for the r^th moment about the mean of the r.v.`):
print(`duration of a fair-coin gambler's ruin game in [0,A]`):
print(`if you start at n. For example, try:`):
print(`Rmoment(n,A,2):`):

elif nops([args])=1 and op(1,[args])=Rmomentp then
print(`Rmomentp(n,A,p,r): Given symbols n and A and a positive`):
print(`integer r>=2 finds an explicit expression`):
print(`(polynomial in n, Laurent polynomial in A)`):
print(`for the r^th moment about the mean of the r.v.`):
print(`duration of `):
print(`gambler's ruin game (with prob. p of getting a dollar), in [0,A]`):
print(`if you start at n. `):
print(` NOTE:p could be symbolic or numeric, if  it is symbolic it is`):
print(`understood that p is NOT 1/2`):
print(`For example, try:`):
print(`Rmomentp(n,A,p,2):`):

elif nops([args])=1 and op(1,[args])=RmomentW then
print(`RmomentW(n,A,r): Given symbols n and A and a positive`):
print(`integer r>=2 finds an explicit expression`):
print(`(polynomial in n and A)`):
print(`for the r^th moment about the mean of the r.v.`):
print(`duration of a WINNING fair-coin gambler's ruin game in [0,A]`):
print(`if you start at n. For example, try:`):
print(`RmomentW(n,A,2):`):

elif nops([args])=1 and op(1,[args])=RmomentSlow then
print(`RmomentSlow(n,A,r): Given symbols n and A and a positive`):
print(`integer r>=2 guesses an explicit expression`):
print(`(polynomial in n, Laurent polynomial in A)`):
print(`for the r^th moment about the mean of the r.v.`):
print(`duration of a fair-coin gambler's ruin game in [0,A]`):
print(`if you start at n. For example, try:`):
print(`RmomentSlow(n,A,2):`):

elif nops([args])=1 and op(1,[args])=LoadedDie then
print(`LoadedDie(P,L): Given a list of non-negative rational`):
print(`numbers P and a list of outcomes L`):
print(`such that the prob. that it will lend at L[i] is P[i] simulates`):
print(`a die with that prob. dist.`):
print(`For example, try: LoadedDie([1/2,1/2],[-1,1]);`):

elif nops([args])=1 and op(1,[args])=PowerToBinom then
print(`PowerToBinom(r): Given r, finds the vector of length`):
print(`r+1 such that x^r=a[0]+a[1]*binomial(x,1)+...+a[r]*binomial(x,r)`):
print(`For example, try:`):
print(`PowerToBinom(2);`):

elif nops([args])=1 and op(1,[args])=Story then
print(`Story(n,A,x,R): inputs symbols n, A, x and a positive integer`):
print(`R, outputs information about the r.v. expected life`):
print(`of a gambler's ruin game in [0,A], starting at n, where at each step`):
print(`with prob. 1/2 one goes one step to the left or one step to the `):
print(` right. It computes (symbolically and rigorously!) everything up to`):
print(`the R^th moment. It also finds the asymptotics for A large and`):
print(`n=Ax and if R>=3 the asymptotic skewness and if R>=4 the`):
print(`asymptotic Kurtosis`):
print(`For example, try:`):
print(`Story(n,A,x,4);`):

elif nops([args])=1 and op(1,[args])=StoryW then
print(`StoryW(n,A,x,R): inputs symbols n, A, x and a positive integer`):
print(`R, outputs information about the r.v. expected life`):
print(`of a WINNING`):
print(`gambler's ruin game in [0,A], starting at n, where at each step`):
print(`with prob. 1/2 one goes one step to the left or one step to the `):
print(` right. It computes (symbolically and rigorously!) everything up to`):
print(`the R^th moment. It also finds the asymptotics for A large and`):
print(`n=Ax and if R>=3 the asymptotic skewness and if R>=4 the`):
print(`asymptotic Kurtosis`):
print(`For example, try:`):
print(`StoryW(n,A,x,4);`):

elif nops([args])=1 and op(1,[args])=Stat then
print(`Stat(P,L,i,A,K,r): The stats of Playing Game(P,L,i,A) K times`):
print(`with info up to the r^th moment.`):
print(`It returns 3 lists:`):
print(`the first one has [prob. of losing, expected time it takes to lose,`):
print(`j^th moment about the mean of time] for j=2..r`):
print(`the second list has the same about winning`):
print(`the third list is the same about finishing`):
print(`For example, try:`):
print(`Stat([1/2,1/2],[-1,1],2,4,100,3);`):

elif nops([args])=1 and op(1,[args])=Tourn then
print(`Tourn(P,L,i,A,K): Playing Game(P,L,i,A) K times`):
print(`returns the list of durations when it lost`):
print(`followed by the list of durations when it won`):
print(`For example, try:`):
print(`Tourn([1/2,1/2],[-1,1],2,4,100);`):

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

#GuessPol1(L,d,n,s0): guesses a polynomial of degree d in n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#GuessPol1([(seq(i,i=1..10),1,n,1);
GuessPol1:=proc(L,d,n,s0) local P,i,a,eq,var:
if d>=nops(L)-s0-2 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=s0..s0+d+3)}:

var:=solve(eq,var):

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

subs(var,P):

end:



#GuessPol(L,n,s0): guesses a polynomial of degree d in n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#GuessPol([(seq(i,i=1..10),n,1);
GuessPol:=proc(L,n,s0) local d,gu:

for d from 0 to nops(L)-s0-3 do
 gu:=GuessPol1(L,d,n,s0):
 if gu<>FAIL then
    RETURN(gu):
 fi:
od:

FAIL:

end:

#LoadedDie(P,L): Given a list of non-negative rational
#numbers P and a list of outcomes L
#such that the prob. that it will lend at L[i] is P[i] simulates
#a die with that prob. dist.
LoadedDie:=proc(P,L) local i,M,P1,ra,r:
if nops(P)<>nops(L) or min(op(P))<0 then
 ERROR(`Bad input`):
fi:

M:=lcm(seq(denom(P[i]),i=1..nops(P))):
P1:=[seq(P[i]*M,i=1..nops(P))]:

ra:=rand(1..M):
r:=ra():
if r<=P1[1] then
 RETURN(L[1]):
fi:
for i from 2 to nops(P1) do
if r>P1[i-1] and r<=convert([op(1..i,P1)],`+`) then
 RETURN(L[i]):
fi:
od:
FAIL:
end:


#Game(P,L,i,A): Given a list of non-negative rational
#numbers P and a list of outcomes L
#such that the prob. that it will lend at L[i] is P[i] simulates
#a game where one exists as soon as one has <=0 dollars
#or >=A dollars and starts at i
#it returns the list of partial outcomes
#For example, try: Game([1/2,1/2],[-1,1],2,4);

Game:=proc(P,L,i,A) local i1,M,P1,ra,r,G,x:
if nops(P)<>nops(L) or min(op(P))<0 then
 ERROR(`Bad input`):
fi:

if convert(P,`+`)<>1 then
 ERROR(`Bad input`):
fi:

M:=lcm(seq(denom(P[i1]),i1=1..nops(P))):
P1:=[seq(P[i1]*M,i1=1..nops(P))]:

ra:=rand(1..M):
G:=[i]:
x:=i:
while x>0 and x<A do
r:=ra():
if r<=P1[1] then
  x:=x+L[1]:
  G:=[op(G),x]:
else
for i1 from 2 to nops(P1) do
if r>P1[i1-1] and r<=convert([op(1..i1,P1)],`+`) then
 x:=x+L[i1]:
 G:=[op(G),x]:
 break:
fi:
od:
fi:
od:
G:
end:





#Tourn(P,L,i,A,K): Playing Game(P,L,i,A) K times
#returns the list of durations when it lost
#followed by the sorted list of durations when it won
#followed by the sorted list of everything
#For example, try:
#Tourn([1/2,1/2],[-1,1],2,4,100);
Tourn:=proc(P,L,i,A,K) local LO,WI,g,i1:
LO:=[]:
WI:=[]:

for i1 from 1 to K do
g:=Game(P,L,i,A):
if g[nops(g)]<=0 then
 LO:=[op(LO),nops(g)-1]:
else
  WI:=[op(WI),nops(g)-1]:
fi:
od:
sort(LO),sort(WI),sort([op(LO),op(WI)]):
end:




#Stat(P,L,i,A,K,r): The stats of Playing Game(P,L,i,A) K times
#with info up to the r^th moment.
#It returns 3 lists:
#the first one has [prob. of losing, expected time it takes to lose,
#j^th moment about the mean of time] for j=2..r
#the second list has the same about winning
#the third list is the same about finishing
#For example, try:
#Stat([1/2,1/2],[-1,1],2,4,100,3);
Stat:=proc(P,L,i,A,K,r) local gu,LO,WI,ALL,probL,probW,probA,avL,avW,avA,i1, j:
if convert(P,`+`)<>1 then
 ERROR(`Bad input`):
fi:
 gu:=Tourn(P,L,i,A,K):
 LO:=gu[1]: WI:=gu[2]: ALL:=[op(gu[1]),op(gu[2])]:
if LO=[] or WI=[] then RETURN(FAIL): fi:

 probL:=nops(LO)/K:
 probW:=nops(WI)/K:
 probA:=nops(ALL)/K:
 avL:=convert(LO,`+`)/nops(LO):
 avW:=convert(WI,`+`)/nops(WI):
 avA:=convert(ALL,`+`)/nops(ALL):
[ [probL,avL,seq(add((LO[i1]-avL)^j,i1=1..nops(LO))/nops(LO),j=2..r)],
 [probW,avW,seq(add((WI[i1]-avW)^j,i1=1..nops(WI))/nops(WI),j=2..r)],
 [probA,avA,seq(add((ALL[i1]-avA)^j,i1=1..nops(ALL))/nops(ALL),j=2..r)]]:
end:
 


#GRwin(P,L,A): Inputs a loaded die P,L and a positive integer
#A and outputs the list whose i^th item is the
#prob. of winning in a Gambler's Ruin game in [0,A] that
#starts at i
#For example, try: GRwin([1/2,1/2],[-1,1],4) ;
GRwin:=proc(P,L,A) local eq,var,mi1,ma1,p,i1,var1,j:
if nops(P)<>nops(L) or not type(A,integer) or A<0 then
 ERROR(`Bad input`):
fi:
mi1:=min(op(L))+1:
ma1:=A+max(op(L))-1:
var:={seq(p[i1],i1=mi1..ma1)}:
eq:={seq(p[i1]=0,i1=mi1..0),seq(p[i1]=1,i1=A..ma1),
seq(p[i1]=add(P[j]*p[i1+L[j]],j=1..nops(P)),i1=1..A-1)}:
var1:=solve(eq,var):
if var1=NULL then
 RETURN(FAIL):
fi:
[seq(subs(var1,p[i1]),i1=1..A)]:
end:


#GRlife(P,L,A): Inputs a loaded die P,L and a positive integer
#A and outputs the list whose i^th item is the
#expected duration of a Gambler's Ruin game in [0,A] that
#starts at i
#For example, try: GRlife([1/2,1/2],[-1,1],4) ;
GRlife:=proc(P,L,A) local eq,var,mi1,ma1,p,i1,var1,j:
option remember:
if nops(P)<>nops(L) or not type(A,integer) or A<0 then
 ERROR(`Bad input`):
fi:
mi1:=min(op(L))+1:
ma1:=A+max(op(L))-1:
var:={seq(p[i1],i1=mi1..ma1)}:
eq:={seq(p[i1]=0,i1=mi1..0),seq(p[i1]=0,i1=A..ma1),
seq(p[i1]=1+add(P[j]*p[i1+L[j]],j=1..nops(P)),i1=1..A-1)}:
var1:=solve(eq,var):
if var1=NULL then
 RETURN(FAIL):
fi:
[seq(subs(var1,p[i1]),i1=1..A)]:
end:


#GRlifeW(P,L,A): Inputs a loaded die P,L and a positive integer
#A and outputs the list whose i^th item is the
#(relative) expected duration of a Winning Gambler's Ruin 
#game in [0,A] that starts at i
#For example, try: GRlifeW([1/2,1/2],[-1,1],4) ;
GRlifeW:=proc(P,L,A) local eq,var,mi1,ma1,p,i1,var1,j,grw,gu:
if nops(P)<>nops(L) or not type(A,integer) or A<0 then
 ERROR(`Bad input`):
fi:
grw:=GRwin(P,L,A):
mi1:=min(op(L))+1:
ma1:=A+max(op(L))-1:
var:={seq(p[i1],i1=mi1..ma1)}:
eq:={seq(p[i1]=0,i1=mi1..0),seq(p[i1]=0,i1=A..ma1),
seq(p[i1]=grw[i1]+add(P[j]*p[i1+L[j]],j=1..nops(P)),i1=1..A-1)}:
var1:=solve(eq,var):
if var1=NULL then
 RETURN(FAIL):
fi:
gu:=[seq(subs(var1,p[i1])/grw[i1],i1=1..A)]:
gu:
end:


#GRlifeBmoment(P,L,A,r): Inputs a loaded die P,L and a positive integer
#A and outputs the list whose i^th item is the
#r-th binomial moment of the r.v. duration of game if it 
#starts at i
#For example, try: GRlifeBmoment([1/2,1/2],[-1,1],4,2) ;
GRlifeBmoment:=proc(P,L,A,r) local eq,var,mi1,ma1,p,i1,var1,j,grw,T:

if nops(P)<>nops(L) or not type(A,integer) or A<0 then
 ERROR(`Bad input`):
fi:
if r=0 then
 RETURN([1$A]):
fi:
grw:=GRlifeBmoment(P,L,A,r-1):
mi1:=min(op(L))+1:
ma1:=A+max(op(L))-1:

for i1 from mi1-1 to 0 do
 T[i1]:=0:
od:

if r=1 then
 T[0]:=1:
fi:

for i1 from A+1 to ma1+1 do
 T[i1]:=0:
od:

for i1 from 1 to A do
 T[i1]:=grw[i1]:
od:

var:={seq(p[i1],i1=mi1..ma1)}:
eq:={seq(p[i1]=0,i1=mi1..0),seq(p[i1]=0,i1=A..ma1),
seq(p[i1]=add(P[j]*(T[i1+L[j]]+p[i1+L[j]]),j=1..nops(P)),i1=1..A-1)}:
var1:=solve(eq,var):
if var1=NULL then
 RETURN(FAIL):
fi:
[seq(subs(var1,p[i1]),i1=1..A)]:
end:

#PowerToBinom(x,r): Given r, finds the vector of length
#r+1 such that x^r=a[0]+a[1]*binomial(x,1)+...+a[r]*binomial(x,r)
#For example, try:
#PowerToBinom(2);
PowerToBinom:=proc(r) local a,i,eq,var,gu,x:
var:={seq(a[i],i=0..r)}:
gu:=expand(x^r-expand(add(a[i]*binomial(x,i),i=0..r))):
eq:={seq(coeff(gu,x,i),i=0..r)}:
var:=solve(eq,var):
[seq(subs(var,a[i]),i=0..r)]:
end:



#GRlifeMoment1(P,L,A,r): Inputs a loaded die P,L and a positive integer
#A and outputs the list whose i^th item is the
#r-th moment of the r.v. duration of game if it 
#starts at i
#For example, try: GRlifeMoment1([1/2,1/2],[-1,1],4,2) ;
GRlifeMoment1:=proc(P,L,A,r) local i1,B,v,j:

v:=PowerToBinom(r):
B:=[seq(GRlifeBmoment(P,L,A,i1),i1=0..r)]:

[seq(add(v[i1+1]*B[i1+1][j],i1=0..r),j=1..A)]:
end:



#GRlifeMoment(P,L,A,r): Inputs a loaded die P,L and a positive integer
#A and outputs the list whose i^th item is the
#r-th moment (about the mean) of the r.v. duration of game if it 
#starts at i
#For example, try: GRlifeMoment([1/2,1/2],[-1,1],4,2) ;
GRlifeMoment:=proc(P,L,A,r) local v,B1,i,k,j:
v:=GRlifeMoment1(P,L,A,1):
B1:=[seq(GRlifeMoment1(P,L,A,i),i=0..r)]:
[seq(add((-1)^(r-k)*binomial(r,k)*v[j]^(r-k)*B1[k+1][j],k=0..r),j=1..A)]:
end:




#CheckGRmoment(P,L,i,A,K,r): Checks the exact value
#by simulation,  the exact value for the r^th moment.
#of the r.v.  duration of a Gambler's Ruin game
#with loaded-die P,L in [0,A] starting at i. For example, try:
#CheckGRmoment([1/2,1/2],[-1,1],2,4,100,3);
CheckGRmoment:=proc(P,L,i,A,K,r) local gu1,gu:
if r<2 then
 ERROR(`r must be >=2 `):
fi:
gu1:=GRlifeMoment(P,L,A,r)[i]:
gu:=evalf(Stat(P,L,i,A,K,r)[3][r+1]):
print(`This checks, via simulation, the exact value`):
print(``, r , `-th moment (about the mean)`):
print(`of the r.v. Duration in Gambler's Ruin in`, [0,A], ` starting at`, i):
print(`with loaded-die`, P,L ):
print(`The exact value is: `, gu1):
print(`and in floating-point:` , evalf(gu1)):
print():
print(`The  empirical value from THIS run of`, K, `games is`):
print(gu):
end:


#CheckGRlife(P,L,i,A,K): Checks the exact value
#by simulation,  the exactation
#of the r.v. expected duration of a Gambler's Ruin game
#with loaded-die P,L in [0,A] starting at i. For example, try:
#CheckGRlife([1/2,1/2],[-1,1],2,4,100);
CheckGRlife:=proc(P,L,i,A,K) local gu1,gu:

gu1:=GRlife(P,L,A)[i]:
gu:=evalf(Stat(P,L,i,A,K,1)[3][2]):
print(`This checks, via simulation, the exact value for the expecation`):
print(`of the r.v. Duration in Gambler's Ruin in`, [0,A], ` starting at`, i):
print(`with loaded-die`, P,L ):
print(`The exact value is: `, gu1):
print(`and in floating-point:` , evalf(gu1)):
print():
print(`The  empirical value from THIS run of`, K, `games is`):
print(gu):
end:

#ExpLife(n,A): Given symbols n and A 
#guesses an explicit expression
#(polynomial in n, Laurent polynomial in A)
#for the expected
#duration of a fair-coin gambler's ruin game in [0,A]
#if you start at n. For example, try:
#ExpLife(n,A);
ExpLife:=proc(n,A) local lu1,lu,A1:

lu1:=
[seq(normal(GuessPol(GRlife([1/2,1/2],[-1,1],A1),n,1)),
A1=7..12)]:

lu:=subs(A=A-6,GuessPol(lu1,A,1)):

if lu=FAIL then
 RETURN(FAIL):
fi:

factor(lu):
end:

#RmomentSlow(n,A,r): Given symbols n and A and a positive
#integer r>=2 guesses an explicit expression
#(polynomial in n, Laurent polynomial in A)
#for the r^th moment about the mean of the r.v.
#duration of a fair-coin gambler's ruin game in [0,A]
#if you start at n. For example, try:
#RmomentSlow(n,A,2);
RmomentSlow:=proc(n,A,r) local lu1,lu,A1:
option remember:

lu1:=
[seq(normal(GuessPol(GRlifeMoment([1/2,1/2],[-1,1],A1,r),n,1)/n/(A1-n)*A1^r),
A1=5+2*r..5+6*r)]:

lu:=subs(A=A-(2*r+4),GuessPol(lu1,A,1)):

if lu=FAIL then
 RETURN(FAIL):
fi:

factor(lu*n*(A-n)/A^r):
end:



#Story(n,A,x,R): inputs symbols n, A, x and a positive integer
#R, outputs information about the r.v. expected life
#of a gambler's ruin game in [0,A], starting at n, where at each step
#with prob. 1/2 one goes one step to the left or one step to the right.
#It computes (symbolically and rigorously!) everything up to
#the R^th moment. It also finds the asymptotics for A large and
#n=Ax and if R>=3 the asymptotic skewness and if R>=4 the
#asymptotic Kurtosis
#For example, try:
#Story(n,A,x,4);
Story:=proc(n,A,x,R) local gu,r,mu,lu,pu,t0:

t0:=time():
if R<4 then
 ERROR(`Last argument must be >=4`):
fi:

print(` Studies in`):
print(` the Classical  Gambler's Ruin Problem `):
print():
print(` by Shalosh B. Ekhad `):
print():
print(`Consider A Player that starts at x=`, n, `and with prob. `):
print(`1/2 walks, at each day, one step to the left or one step to`):
print(` the right, until it either reaches x=0 or x=`, A):
print(`Let D be the random variable: duration of the game`):
print(`The Expected value of D is`):
print(BmomentR(1,n,A)):

for r from 2 to R do
 gu[r]:=Rmoment(n,A,r);
 mu[r]:=normal(subs(n=A*x,gu[r])):
 lu[r]:=coeff(mu[r],A,degree(mu[r],A))*A^degree(mu[r],A):
 print(`The `, r, `-th moment (about the mean) of D is`):
 print(factor(gu[r])):
 print(`and setting n=`, x*A , `it is `):
 print(factor(mu[r])):
 print(`and asymptotically it is, as`, A, `goes to infinity`):
 print(lu[r]):
 print(`If it starts at the middle, i.e.`, x=1/2, `it is`):
 print(subs(x=1/2,lu[r])):
 print(`which is roughly`):
   print(evalf(subs(x=1/2,lu[r]))):
od:

print(`The asymptotic Skewness is`):
pu:=normal(lu[3]^2/lu[2]^3);
pu:=sqrt(pu):
print(pu):
print(`and when`, x=1/2 , `it is`):
print(subs(x=1/2,pu)):
print(`which is roughly`, evalf(subs(x=1/2,pu))):

print(`The asymptotic Kurtosis is`):
pu:=normal(lu[4]/lu[2]^2);
print(pu):
print(`and when`, x=1/2 , `it is`):
print(subs(x=1/2,pu)):
print(`which is roughly`, evalf(subs(x=1/2,pu))):

print(`The sequence of r^th-root of r^th moment for x=1/2,`):
print(` starting at the second is, divided by`, A^2, ` is :`):
print([seq(evalf(simplify((subs(x=1/2,normal
(lu[r]/A^(2*r))))^(1/r))),r=2..R)] ):

print(`The whole thing took`, time()-t0, `seconds  of CPU time`):
end:





####Begin from BVP 


#BVh(p,x,n,A,a,L,Left1,Right1):
#given a Laurent polynomial p(x) where the multipliicity
#of its roots are given by the list L
#Solves (symbolically) in terms of the roots
#corresponding to L labelled by the indexed variable
#a, and two lists, Left1, Left2, of length
#-ldegree(p,x) and degree(p,x) resp.
#the homog. boundary value problem
#p(N)f(n)=0 (where N is the shift in n: Ng(n):=g(n+1))
#subject to the boundary conditions
#f(-i)=Left1[i+1], i=0..nops(Left1)-1
#f(A+i)=Right1[A+i], i=0..nops(Right1)-1
#For example, try:
#BVh(1-(x+1/x)/2,x,n,A,a,[2],[0],[1]);
BVh:=proc(p,x,n,A,a,L,Left1,Right1) local J1,J2,f,b,i,j,i1,eq,var:
J1:=-ldegree(p,x):J2:=degree(p,x):

if nops(Left1)<>J1 or nops(Right1)<>J2 then
 ERROR(`Bad input`):
fi:

if convert(L,`+`)<>J1+J2 then
 ERROR(`Bad input`):
fi:

f:=add(add(b[i,j]*a[i]^n*n^j,j=0..L[i]-1),i=1..nops(L)):
var:={seq(seq(b[i,j],j=0..L[i]-1),i=1..nops(L))}:
eq:={seq(subs(n=-i1,f)=Left1[i1+1],i1=0..J1-1),
  seq(subs(n=A+i1,f)=Right1[i1+1],i1=0..J2-1)}:
var:=solve(eq,var):

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

subs(var,f):

end:

##########From BVP##################




#GRwinS(P,L,n,A): Inputs a loaded die P,L and SYMBOLS n and
#A and outputs the expression in terms of the
#roots of a certain polynomial for the explicit expression
#for the prob. of winning in a Gambler's Ruin game in [0,A] that
#starts at i
#For example, try: GRwinS([1/2,1/2],[-1,1],n,A) ;
GRwinS:=proc(P,L,n,A) local p,x,mu,L1,i1,j,j1,Ones,gu,a,Left1,Right1,A1:
if nops(P)<>nops(L) then
 ERROR(`Bad input`):
fi:

if convert(P,`+`)<>1 then
 ERROR(`Bad input`):
fi:

p:=1-add(P[j]*x^L[j],j=1..nops(P)):

mu:=[solve(p,x)]:

Ones:=[]:
for i1 from 1 to nops(mu) do
 if mu[i1]=1 then
  Ones:=[op(Ones),i1]:
 fi:
od:


mu:=[1$nops(Ones),op(1..Ones[1]-1,mu),
seq(op(Ones[j1]+1..Ones[j1+1]-1,mu),j1=1..nops(Ones)-1),
op(Ones[nops(Ones)]+1..nops(mu),mu)]:
if nops([mu])-nops({op(mu)}) >=2 then
  print(`Not implemented`):
  RETURN(FAIL):
fi:



if nops(Ones)=1 then
 L1:=[1$nops(mu)]:
elif nops(Ones)=2 then
 L1:=[2,1$(nops(mu)-2)]:
 mu:=[1,op(3..nops(mu),mu)]:
else
 RETURN(FAIL):
fi:


Left1:=[0$(-ldegree(p,x))]:
Right1:=[1$degree(p,x)]:

gu:=BVh(p,x,n,A,a,L1,Left1,Right1):
gu:=subs({seq(a[j]=mu[j],j=1..nops(mu))},gu):

gu:=simplify(gu):

if {seq(evalb(GRwin(P,L,A1)=[seq(simplify(
subs({n=i1,A=A1},gu)),i1=1..A1)]), A1=1..20)}<>{true} then
 print(gu, `is wrong`):
 RETURN(FAIL):
fi:
gu:

end:
 

#ExpLifeW(n,A): Given symbols n and A 
#guesses an explicit expression
#(polynomial in n, Laurent polynomial in A)
#for the expected
#duration of a WINNING  fair-coin gambler's ruin game in [0,A]
#if you start at n. For example, try:
#ExpLifeW(n,A);
ExpLifeW:=proc(n,A) local lu1,lu,A1:

lu1:=
[seq(normal(GuessPol(GRlifeW([1/2,1/2],[-1,1],A1),n,1)),
A1=7..12)]:

lu:=subs(A=A-6,GuessPol(lu1,A,1)):

if lu=FAIL then
 RETURN(FAIL):
fi:

factor(lu):
end:


#BmomentR(r,n,A): the r^th Binomial moment of the r.v. duration
#of game in a Gambler's Ruin problem in [0,A] starting at $n$
#For example, try: BmomentR(2,n,A);
BmomentR:=proc(r,n,A)
local gu1,c1,c2,a,d,P,eq,var,P1,i,lu,var1:
option remember:

if r=0 then
 RETURN(1):
fi:

gu1:=BmomentR(r-1,n,A):
d:=degree(gu1,n)+2:

P:=add(a[i]*n^i,i=2..d):
var:={seq(a[i],i=2..d)}:
lu:=P-(1/2)* subs(n=n-1,P)-(1/2)* subs(n=n+1,P)-
(1/2)* subs(n=n-1,gu1)-(1/2)* subs(n=n+1,gu1):
lu:=expand(lu):
eq:={seq(coeff(lu,n,i),i=0..degree(lu,n))}:
var:=solve(eq,var):

P:=subs(var,P):

P:=subs({seq(a[i]=1,i=2..d)},P):

P1:=P+c1+c2*n:
var1:=solve({subs(n=0,P1),subs(n=A,P1)},{c1,c2}):
P1:=subs(var1,P1):
factor(P1):
end:


#Rmoment(n,A,r): Given symbols n and A and a positive
#integer r>=2 finds an explicit expression
#(polynomial in n, Laurent polynomial in A)
#for the r^th moment about the mean of the r.v.
#duration of a fair-coin gambler's ruin game in [0,A]
#if you start at n. For example, try:
#Rmoment(n,A,2);
Rmoment:=proc(n,A,r) local a,k:
option remember:
a:=Rmoment1(n,A,1):
factor(normal(
add((-1)^(r-k)*binomial(r,k)*a^(r-k)*Rmoment1(n,A,k),k=0..r))):
end:


#Rmoment1(n,A,r): Given symbols n and A and a positive
#integer r>=2 finds an explicit expression
#(polynomial in n, Laurent polynomial in A)
#for the r^th moment (NOT about the mean) of the r.v.
#duration of a fair-coin gambler's ruin game in [0,A]
#if you start at n. For example, try:
#Rmoment1(n,A,2);
Rmoment1:=proc(n,A,r) local v,i1:
v:=PowerToBinom(r):
normal(add(v[i1+1]*BmomentR(i1,n,A),i1=0..r));
end:


#CheckGRmomentFair(i,A,K,r): Checks the exact value
#by simulation,  the exact value for the r^th moment.
#of the r.v.  duration of a Gambler's Ruin game
#with a fair coin, in [0,A], starting at i. For example, try:
#CheckGRmomentFair(2,4,100,3);
CheckGRmomentFair:=proc(i,A,K,r) local gu1,gu,n,A1,P,L:
if r<2 then
 ERROR(`r must be >=2 `):
fi:
P:=[1/2,1/2]:
L:=[-1,1]:
gu1:=subs({n=i,A1=A},Rmoment(n,A1,r)):
gu:=evalf(Stat(P,L,i,A,K,r)[3][r+1]):
print(`This prpgram checks, via simulation, the exact value of the `):
print(``, r , `-th moment (about the mean)`):
print(`of the r.v. Duration in Gambler's Ruin in`, [0,A], ` starting at`, i):
print(`with fair coin.`):
print(`The exact value, from the formula is: `, gu1):
print(`and in floating-point:` , evalf(gu1)):
print():
print(`The  empirical value from THIS run of`, K, `games is`):
print(gu):
end:








#BmomentRw1(r,n,A): the pre-r^th Binomial moment of the r.v. duration
#of WINNING game in a Gambler's Ruin problem in [0,A] starting at $n$
#For example, try: BmomentRw(2,n,A);
BmomentRw1:=proc(r,n,A)
local gu1,c1,c2,a,d,P,eq,var,P1,i,lu,var1:
option remember:

if r=0 then
 RETURN(n/A):
fi:

gu1:=BmomentRw1(r-1,n,A):
d:=degree(gu1,n)+2:

P:=add(a[i]*n^i,i=2..d):
var:={seq(a[i],i=2..d)}:
lu:=P-(1/2)* subs(n=n-1,P)-(1/2)* subs(n=n+1,P)-
(1/2)* subs(n=n-1,gu1)-(1/2)* subs(n=n+1,gu1):
lu:=expand(lu):
eq:={seq(coeff(lu,n,i),i=0..degree(lu,n))}:
var:=solve(eq,var):

P:=subs(var,P):

P:=subs({seq(a[i]=1,i=2..d)},P):

P1:=P+c1+c2*n:
var1:=solve({subs(n=0,P1),subs(n=A,P1)},{c1,c2}):
P1:=subs(var1,P1):
factor(P1):
end:


#BmomentRw(r,n,A): the r^th Binomial moment of the r.v. duration
#of WINNING game in a Gambler's Ruin problem in [0,A] starting at $n$
#For example, try: BmomentRw(2,n,A);

BmomentRw:=proc(r,n,A)
normal(BmomentRw1(r,n,A)/(n/A)):
end:




#Rmoment1w(n,A,r): Given symbols n and A and a positive
#integer r>=2 finds an explicit expression
#(polynomial in n, Laurent polynomial in A)
#for the r^th moment (NOT about the mean) of the r.v.
#duration of a WINNING fair-coin gambler's ruin game in [0,A]
#if you start at n. For example, try:
#Rmoment1w(n,A,2);
Rmoment1w:=proc(n,A,r) local v,i1:
v:=PowerToBinom(r):
normal(add(v[i1+1]*BmomentRw(i1,n,A),i1=0..r));
end:



#RmomentW(n,A,r): Given symbols n and A and a positive
#integer r>=2 finds an explicit expression
#(polynomial in n, Laurent polynomial in A)
#for the r^th moment about the mean of the r.v.
#duration of a WINNING fair-coin gambler's ruin game in [0,A]
#if you start at n. For example, try:
#RmomentW(n,A,2);
RmomentW:=proc(n,A,r) local a,k:
option remember:
a:=Rmoment1w(n,A,1):
factor(normal(
add((-1)^(r-k)*binomial(r,k)*a^(r-k)*Rmoment1w(n,A,k),k=0..r))):
end:

#StoryW(n,A,x,R): inputs symbols n, A, x and a positive integer
#R, outputs information about the r.v. expected life
#of a WINNING
#gambler's ruin game in [0,A], starting at n, where at each step
#with prob. 1/2 one goes one step to the left or one step to the right.
#It computes (symbolically and rigorously!) everything up to
#the R^th moment. It also finds the asymptotics for A large and
#n=Ax and if R>=3 the asymptotic skewness and if R>=4 the
#asymptotic Kurtosis
#For example, try:
#StoryW(n,A,x,4);
StoryW:=proc(n,A,x,R) local gu,r,mu,lu,pu,t0:

t0:=time():
if R<4 then
 ERROR(`Last argument must be >=4`):
fi:

print(` Studies in`):
print(` the Classical  Gambler's Ruin Problem `):
print():
print(` by Shalosh B. Ekhad `):
print():
print(`Consider A Player that starts at x=`, n, `and with prob. `):
print(`1/2 walks, at each day, one step to the left or one step to`):
print(` the right, until it either reaches x=0 or x=`, A):
print(`Let D be the random variable: duration of a WINNING game`):
print(`The Expected value of D is`):
print(BmomentRw(1,n,A)):

for r from 2 to R do
 gu[r]:=RmomentW(n,A,r);
 mu[r]:=normal(subs(n=A*x,gu[r])):
 lu[r]:=coeff(mu[r],A,degree(mu[r],A))*A^degree(mu[r],A):
 print(`The `, r, `-th moment (about the mean) of D is`):
 print(factor(gu[r])):
 print(`and setting n=`, x*A , `it is `):
 print(factor(mu[r])):
 print(`and asymptotically it is, as`, A, `goes to infinity`):
 print(lu[r]):
 print(`If it starts at the middle, i.e.`, x=1/2, `it is`):
 print(subs(x=1/2,lu[r])):
 print(`which is roughly`):
   print(evalf(subs(x=1/2,lu[r]))):
od:

print(`The asymptotic Skewness is`):
pu:=normal(lu[3]^2/lu[2]^3);
pu:=sqrt(pu):
print(pu):
print(`and when`, x=1/2 , `it is`):
print(subs(x=1/2,pu)):
print(`which is roughly`, evalf(subs(x=1/2,pu))):

print(`The asymptotic Kurtosis is`):
pu:=normal(lu[4]/lu[2]^2);
print(pu):
print(`and when`, x=1/2 , `it is`):
print(subs(x=1/2,pu)):
print(`which is roughly`, evalf(subs(x=1/2,pu))):

print(`The sequence of r^th-root of r^th moment for x=1/2,`):
print(` starting at the second is, divided by`, A^2, ` is :`):
print([seq(evalf(simplify((subs(x=1/2,normal
(lu[r]/A^(2*r))))^(1/r))),r=2..R)] ):

print(`The whole thing took`, time()-t0, `seconds  of CPU time`):
end:



#BmomentRp(r,p,n,A): the r^th Binomial moment of the r.v. duration
#of game in a Gambler's Ruin problem in [0,A] starting at $n$
#and with prob. of winning a dollar p 
#For example, try: BmomentRp(2,p,n,A);
BmomentRp:=proc(r,p,n,A) local gu:
gu:=BmomentRp1(r,p,n,A) :
gu[1]+gu[2]*((1-p)/p)^n:
end:

#BmomentRp1(r,p,n,A): the r^th Binomial moment of the r.v. duration
#of game in a Gambler's Ruin problem in [0,A] starting at $n$
#and with prob. of winning a dollar p 
#For example, try: BmomentRp1(2,p,n,A);
BmomentRp1:=proc(r,p,n,A)
local gu1,gu2,c1,c2,a,eq,var,P1,i,lu1,var1,P2, d1, d2, gu, lu2:
option remember:
if p=1/2 then
 RETURN([BmomentR(r,n,A),0]):
fi:

if r=0 then
 RETURN([1,0]):
fi:

gu:=BmomentRp1(r-1,p,n,A):
gu1:=gu[1]:
gu2:=gu[2]:

d1:=degree(gu1,n)+1:
P1:=add(a[i]*n^i,i=1..d1):
var:={seq(a[i],i=1..d1)}:
lu1:=P1-(1-p)* subs(n=n-1,P1)-(p)* subs(n=n+1,P1)-
(1-p)* subs(n=n-1,gu1)-(p)* subs(n=n+1,gu1):
lu1:=expand(lu1):
eq:={seq(coeff(lu1,n,i),i=0..degree(lu1,n))}:
var:=solve(eq,var):
P1:=subs(var,P1):

if gu2=0 then
 d2:=1:
else
d2:=degree(gu2,n)+1:
fi:

P2:=add(a[i]*n^i,i=1..d2):
var:={seq(a[i],i=1..d2)}:
lu2:=P2-p* subs(n=n-1,P2)-(1-p)* subs(n=n+1,P2)-
(p)* subs(n=n-1,gu2)-(1-p)* subs(n=n+1,gu2):
lu2:=expand(lu2):
eq:={seq(coeff(lu2,n,i),i=0..degree(lu2,n))}:
var:=solve(eq,var):

P2:=subs(var,P2):

#Q:=P1+P2*((1-p)/p)^n+c1+c2*((1-p)/p)^n:
eq:={subs(n=0, P1+P2)+c1+c2,
subs(n=A, P1)+subs(n=A,P2)*((1-p)/p)^A+c1+c2*((1-p)/p)^A }:

var1:=solve(eq,{c1,c2}):

[subs(var1,P1)+subs(var1,c1),subs(var1,P2)+subs(var1,c2)]:


end:



#Rmoment1p(n,A,p,r): Given symbols n and A and a positive
#integer r>=2 finds an explicit expression
#(polynomial in n, Laurent polynomial in A)
#for the r^th moment (NOT about the mean) of the r.v.
#duration of a loaded coin (with prob. p of winning a dollar)
#gambler's ruin game in [0,A]
#if you start at n. For example, try:
#Rmoment1p(n,A,p,2);
Rmoment1p:=proc(n,A,p,r) local v,i1:
v:=PowerToBinom(r):
normal(add(v[i1+1]*BmomentRp(i1,p,n,A),i1=0..r));
end:

#Rmomentp(n,A,p,r): Given symbols n and A and a positive
#integer r>=2 finds an explicit expression
#(polynomial in n, Laurent polynomial in A)
#for the r^th moment about the mean of the r.v.
#duration of a loaded coin (with prob. of winning a dollar
#p) gambler's ruin game in [0,A]
#if you start at n. For example, try:
#Rmomentp(n,A,p,2);
Rmomentp:=proc(n,A,p,r) local a,k:
option remember:
a:=Rmoment1p(n,A,p,1):
factor(normal(
add((-1)^(r-k)*binomial(r,k)*a^(r-k)*Rmoment1p(n,A,p,k),k=0..r))):
end:


#CheckGRmomentLoaded(i,A,p,K,r): Checks the exact value
#by simulation,  the exact value for the r^th moment.
#of the r.v.  duration of a Gambler's Ruin game
#with loaded coin with prob. of winning
#a dollar being p, in [0,A] starting at i. For example, try:
#CheckGRmomentLoaded(2,4,1/3,100,3);
CheckGRmomentLoaded:=proc(i,A,p,K,r) local gu1,gu,n,A1,P,L,ku:
if r<2 then
 ERROR(`r must be >=2 `):
fi:
P:=[1-p,p]:
L:=[-1,1]:
gu1:=subs({n=i,A1=A},Rmomentp(n,A1,p,r)):
ku:=Stat(P,L,i,A,K,r):
if ku=FAIL then
 print(`One of them is empty`):
RETURN(FAIL):
fi:
gu:=evalf(ku[3][r+1]):
print(`This prpgram checks, via simulation, the exact value of the `):
print(``, r , `-th moment (about the mean)`):
print(`of the r.v. Duration in Gambler's Ruin in`, [0,A], ` starting at`, i):
print(`with fair coin.`):
print(`The exact value, from the formula is: `, gu1):
print(`and in floating-point:` , evalf(gu1)):
print():
print(`The  empirical value from THIS run of`, K, `games is`):
print(gu):
end:



#PrAtnF(P,L,n,A,t): Inputs a loaded die P,L and a positive integer
#A and  a positive integer r, and an integer n between 0 and A,
#outputs the probability that the game will FINISH exactly after
#t tosses.
#For example, try: PrAtnF([1/2,1/2],[-1,1],2,4,2) ;
PrAtnF:=proc(P,L,n,A,t) local j:
option remember:
if nops(P)<>nops(L) or not type(A,integer) or A<0  or t<0 then

 ERROR(`Bad input`):
fi:

if t=0 then
 if n<=0 or n>=A then
   RETURN(1):
  else
    RETURN(0):
 fi:
fi:
if n<=0 or n>=A then
 if t=0 then
   RETURN(1):
  else
   RETURN(0):
  fi:
fi:

add(P[j]*PrAtnF(P,L,n+L[j],A,t-1),j=1..nops(P)):
end:



#PrAtnW(P,L,n,A,t): Inputs a loaded die P,L and a positive integer
#A and  a positive integer r, and an integer n between 0 and A,
#outputs the probability that the game first player
#will WIN exactly after
#t tosses.
#For example, try: PrAtnW([1/2,1/2],[-1,1],2,4,2) ;
PrAtnW:=proc(P,L,n,A,t) local j:
option remember:
if nops(P)<>nops(L) or not type(A,integer) or A<0  or t<0 then

 ERROR(`Bad input`):
fi:

if t=0 then
 if  n>=A then
   RETURN(1):
  else
    RETURN(0):
 fi:
fi:

if  n>=A then
 if t=0 then
   RETURN(1):
  else
   RETURN(0):
  fi:
fi:

if n<=0 then
 RETURN(0):
fi:

add(P[j]*PrAtnW(P,L,n+L[j],A,t-1),j=1..nops(P)):
end:



#PrAtnFad(P,L,n,A,t): Probab. of finishing within <=t tosses
PrAtnFad:=proc(P,L,n,A,t) local t1:
add(PrAtnF(P,L,n,A,t1),t1=0..t): 
end:

#PrAtnWad(P,L,n,A,t): Probab. of winning within <=t tosses
PrAtnWad:=proc(P,L,n,A,t) local t1:
add(PrAtnW(P,L,n,A,t1),t1=0..t): 
end:

#PrAtnFadBM(P,L,n,A,t,r): Binomial Moment of finishing within <=t tosses
PrAtnFadBM:=proc(P,L,n,A,t,r) local t1:
add(binomial(t1,r)*PrAtnF(P,L,n,A,t1),t1=0..t): 
end:


#StoryTeX(n,A,x,R): Like Story(n,A,x,R) but with latex
StoryTeX:=proc(n,A,x,R) local gu,r,mu,lu,pu,t0:

t0:=time():
if R<4 then
 ERROR(`Last argument must be >=4`):
fi:

print(` Studies in`):
print(` the Classical  Gambler's Ruin Problem `):
print():
print(` by Shalosh B. Ekhad `):
print():
print(`Consider A Player that starts at x=`, n, `and with prob. `):
print(`1/2 walks, at each day, one step to the left or one step to`):
print(` the right, until it either reaches x=0 or x=`, A):
print(`Let D be the random variable: duration of the game`):
print(`The Expected value of D is`):
lprint(latex(BmomentR(1,n,A))):
for r from 2 to R do
 gu[r]:=Rmoment(n,A,r);
 mu[r]:=normal(subs(n=A*x,gu[r])):
 lu[r]:=coeff(mu[r],A,degree(mu[r],A))*A^degree(mu[r],A):
 print(`The `, r, `-th moment (about the mean) of D is`):
 lprint(latex(factor(gu[r]))):
 print(`and setting n=`, x*A , `it is `):
lprint(latex(factor(mu[r]))):
 print(`and asymptotically it is, as`, A, `goes to infinity`):
 lprint(latex(lu[r])):
 print(`If it starts at the middle, i.e.`, x=1/2, `it is`):
 lprint(latex(subs(x=1/2,lu[r]))):
 print(`which is roughly`):
   lprint(latex(evalf(subs(x=1/2,lu[r])))):
od:

print(`The asymptotic Skewness is`):
pu:=normal(lu[3]^2/lu[2]^3);
pu:=sqrt(pu):
lprint(latex(pu)):
print(`and when`, x=1/2 , `it is`):
lprint(latex(subs(x=1/2,pu))):
print(`which is roughly`):
lprint(evalf(subs(x=1/2,pu))):

print(`The asymptotic Kurtosis is`):
pu:=normal(lu[4]/lu[2]^2);
lprint(latex(pu)):
print(`and when`, x=1/2 , `it is`):
lprint(latex(subs(x=1/2,pu))):
print(`which is roughly`, evalf(subs(x=1/2,pu))):

print(`The sequence of r^th-root of r^th moment for x=1/2,`):
print(` starting at the second is, divided by`, A^2, ` is :`):
print([seq(evalf(simplify((subs(x=1/2,normal
(lu[r]/A^(2*r))))^(1/r))),r=2..R)] ):

print(`The whole thing took`, time()-t0, `seconds  of CPU time`):
end: