######################################################################
## GR3moms.txt Save this file as  GR3moms.txt  to use it,            #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read  GR3moms.txt <Enter>                         #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(`First Written: Sept. 2023: tested for Maple 2020 `):
print(`Version  : Sept. 2023  `):
print():
print(`This is GR3moms.txt, A Maple package to find explicit expressions for many moments of the Duration of a 3-player Gambler's Ruin.`):
print(`accompanying Shalsoh B. Ekhad and Doron Zeilberger's  article: `):
print(`"Explicit Formulas for the  Moments of the Duration of a (fair) 3-Players Gambler's Ruin" `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/GR3moms.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():

print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();".`):
print(` For specific help type "ezra(procedure_name);" `):

print(`----------------------`):

print(`For a list of the supporting functions type: ezra1();`):
print():
print(` For specific help type "ezra(procedure_name);" `):
 print(`----------------------`):


print(`For a list of the  functions  regarding 2-players Gambler's Ruin type: ezra2();`):
print():
print(` For specific help type "ezra(procedure_name);" `):
 print(`----------------------`):




ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: f3iabcT2, f3iabcTie, FL `):
print(``):

else
ezra(args):

fi:


end:


ezra2:=proc()
if args=NULL then

print(`The  procedures regarding 2-players Gambler's Ruin are:`):
print(` AllMoms2, AllMoms2T, AllMoms2am, AllMoms2amT, AllMoms2amS,  AllMoms2amST ,AllMoms2L, AllMoms2Le, f2iab, M2iab, M2iabT `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` GR3moms.txt: A Maple package for  computing moments of the duration of a 3-player (and also 2-player) Gambler's Ruin `):

print(`The MAIN procedures are: AllMoms3, AllMoms3am, AllMoms3amS, AllMom3L, AllMoms3Le, f3iabc,  M3iabc `):
print(``):

elif nargs=1 and args[1]=AllMoms2 then
print(`AllMoms2(K,a,b): The list of the first K (straignt) moments for the duration of a 2-player Gambler's Ruin with initial capitals a,b. Try:`):
print(`AllMoms2(6,a,b);`):


elif nargs=1 and args[1]=AllMoms2T then
print(`AllMoms2T(K,a,b): The list of the first K (straignt) moments for the duration of a 2-player Gambler's Ruin with TIES, with initial capitals a,b. Try:`):
print(`AllMoms2T(6,a,b);`):

elif nargs=1 and args[1]=AllMoms2am then
print(`AllMoms2am(K,a,b): The list of the first K  moments ABOUT THE MEAN, for the duration of a 2-player Gambler's Ruin with initial capitals a,b`):
print(`AllMoms2am(6,a,b);`):
print(`The first and second entries are the same as those of AllMoms2(K,a,b), in other words, the expectation and variance `):


elif nargs=1 and args[1]=AllMoms2amT then
print(`AllMoms2amT(K,a,b): The list of the first K  moments ABOUT THE MEAN, for the duration of a 2-player Gambler's Ruin with TIES, with initial capitals a,b`):
print(`AllMoms2amT(6,a,b);`):
print(`The first and second entries are the same as those of AllMoms2T(K,a,b,c), in other words, the expectation and variance `):

elif nargs=1 and args[1]=AllMoms2amS then
print(`AllMoms2amS(K,a,b): The list of the first K  SCALED moments ABOUT THE MEAN, for the duration of a 2-player Gambler's Ruin with initial capitals a,b`):
print(`The first and second entries are the same as those of AllMoms2am(K,a,b), in other words, the expectation and variance. Try: `):
print(`AllMoms2amS(6,a,b);`):


elif nargs=1 and args[1]=AllMoms2amST then
print(`AllMoms2amST(K,a,b): The list of the first K  SCALED moments ABOUT THE MEAN, for the duration of a 2-player Gambler's Ruin with TIES, with initial capitals a,b`):
print(`The first and second entries are the same as those of AllMoms2am(K,a,b), in other words, the expectation and variance. Try: `):
print(`AllMoms2amST(6,a,b);`):



elif nargs=1 and args[1]=AllMoms2L then
print(`AllMoms2L(K,x1,x2): The list of length K such that if the initial capitals are x1*N, x2*N, in a 2-player gambler's ruin,`):
print(` the first enrty is the limit, as N goes to infinity, of the expectaion divided by N^2`):
print(` the second enrty is the limit, as N goes to infinity, of the variance divided by N^4`):
print(`the remaining entries are the limits of the i-th scaled moments for i from 3 to K, as N goes to infinity. Try:`):
print(`AllMoms2L(6,x1,x2);`):



elif nargs=1 and args[1]=AllMoms2Le then
print(`AllMoms2Le(K): Inputs a positive integer K>=2, and outputs a list of NUMBERS of length K such that,`):
print(`in a fair 2-player gambler's ruin, where both players have initial capital a dollars`):
print(` the first entry is the expected duration divided by a^2`):
print(`the first entry is the limit of the variance of the duration divided by a^4, as a goes to infinity`):
print(`The third through the K-th entry are the limits of the scaled moments, as a goes to infinity. Try:`):
print(`AllMoms2Le(10);`):

elif nargs=1 and args[1]=AllMoms3 then
print(`AllMoms3(K,a,b,c): The list of the first K SCALED (straignt) moments for the duration of a 3-player Gambler's Ruin with initial capitals a,b,c. Try:`):
print(`AllMoms3(6,a,b,c);`):

elif nargs=1 and args[1]=AllMoms3am then
print(`AllMoms3am(K,a,b,c): The list of the first K moments ABOUT THE MEAN,for the duration of a 3-player Gambler's Ruin with initial capitals a,b,c. Try:`):
print(`AllMoms3am(6,a,b,c);`):


elif nargs=1 and args[1]=AllMoms3amS then
print(`AllMoms3amS(K,a,b,c): The list of the first SCALED K moments ABOUT THE MEAN,for the duration of a 3-player Gambler's Ruin with initial capitals a,b,c`):
print(`The first and second entries are the same as those of AllMoms3am(K,a,b,c), in other words, the expectation and variance. Try: `):
print(`AllMoms3amS(6,a,b,c);`):


elif nargs=1 and args[1]=AllMoms3L then
print(`AllMoms3L(K,x1,x2,x3): The list of length K such that if the initial capitals are x1*N, x2*N, x3*N, in a 3-player gambler's ruin,`):
print(` the first enrty is the limit, as N goes to infinity, of the expectaion divided by N^2`):
print(` the second enrty is the limit, as N goes to infinity, of the variance divided by N^4`):
print(`the remaining entries are the limits of the i-th scaled moments for i from 3 to K, as N goes to infinity. Try:`):
print(`AllMoms3L(6,x1,x2.x3);`):


elif nargs=1 and args[1]=AllMoms3Le then
print(`AllMoms3Le(K): Inputs a positive integer K>=2, and outputs a list of NUMBERS of length K such that,`):
print(`in a fair 3-player gambler's ruin, where all three players have initial capital a dollars`):
print(` the first entry is the expected duration divided by a^2`):
print(`the first entry is the limit of the variance of the duration divided by a^4, as a goes to infinity`):
print(`The third through the K-th entry are the limits of the scaled moments, as a goes to infinity. Try:`):
print(`AllMoms3Le(10);`):

elif nargs=1 and args[1]=FL then
print(`FL(f,L): Applies the fair Gambler's ruin operator with nops(L) players  to the expression f in the variables of L. Try:`):
print(`FL(a*b*c,[a,b,c]);`):

elif nargs=1 and args[1]=f2iab then
print(`f2iab(i,a,b): The explicit expression for the i-th binomial moment of the duration of a 2-player Gambler's ruin with initial capitals a,b. Try:`):
print(`f2iab(2,a,b);`):


elif nargs=1 and args[1]=f3iabc then
print(`f3iabc(i,a,b,c): The explicit expression for the i-th binomial moment of the duration of a 3-player Gambler's ruin with initial capitals a,b,c until one of the goes broke. Try:`):
print(`f3iabc(4,a,b,c);`):


elif nargs=1 and args[1]=f3iabcTie then
print(`f3iabcTie(i,a,b,c): The explicit expression (if it exists) for the i-th binomial moment of the duration of a 3-player Gambler's ruin with  TIES (see p. 42 of Jiri Andel's book) with initial capitals a,b,c where ties are allowed).`):
print(`Note that if FAILS for i>1, but works for i=1. Try:`):
print(`f3iabcTie(1,a,b,c);`):

elif nargs=1 and args[1]=f3iabcT2 then
print(`f3iabcT2(i,a,b,c): The explicit expression for the i-th binomial moment of the duration of a 3-player Gambler's ruin with initial capitals a,b,c, until two of them are broke, or FAIL.`):
print(`Note that there does not seem to be a nice answer if i>1, hence you would get FAIL, but it works for the expectation. Try: `):
print(`f3iabcT2(1,a,b);`):

elif nargs=1 and args[1]=M2iab then
print(`M2iab(i,a,b): The i-th moment of the duration of a 2-player gambler's ruin with starting capitals a,b. Try:`):
print(` M2iab(4,a,b); `):


elif nargs=1 and args[1]=M3iabc then
print(`M3iabc(i,a,b,c): The i-th moment of the duration of a 3-player gambler's ruin with starting capitals a,b,c,  until one of the goes broke. Try:`):
print(` M3iab(4,a,b,c): `):

elif nargs=1 and args[1]=Paper2 then
print(`Paper2(R,a1,a2): outputs a paper about the first R moments of the duration of 2-player gambler's ruin, with initial capitals a1,a2, until one of them goes broke. Try:`):
print(`Paper2(6,A,B);`):


elif nargs=1 and args[1]=Paper2terse then
print(`Paper2terse(R,a1,a2): outputs a paper about the first R moments of the duration of 2-player gambler's ruin, with initial capitals a1,a2, until one of them goes broke, but tersely. Try:`):
print(`Paper2terse(6,A,B);`):


elif nargs=1 and args[1]=Paper2terseT then
print(`Paper2terseT(R,a1,a2): outputs a paper about the first R moments of the duration of 2-player gambler's ruin with TIES, with initial capitals a1,a2, until one of them goes broke, but tersely. Try:`):
print(`Paper2terseT(6,A,B);`):


elif nargs=1 and args[1]=Paper3 then
print(`Paper3(R,a1,a2,a3): outputs a paper about the first R moments of the duration of 3-player gambler's ruin with initial capitals, a1,a2,a3, until one of them goes broke. Try:`):
print(`Paper3(6,A,B,C);`):


elif nargs=1 and args[1]=Paper3terse then
print(`Paper3terse(R,a1,a2,a3): outputs a TERSE paper about the first R moments of the duration of 3-player gambler's ruin with initial capitals, a1,a2,a3, until one of them goes broke. Try:`):
print(`Paper3terse(6,A,B,C);`):


else
 print(`There is no such thing as`, args):

fi:


end:

ez2:=proc(): print(`  Paper2(R,a1,a2)  `): end:

with(combinat):

#f2iab(i,a,b): The explicit expression for the i-th binomial moment of the duration of a 2-player Gambler's ruin with initial capitals a,b. Try:
#f2iab(2,a,b);
f2iab:=proc(i,a,b) local eq,var,P,c,P1,P0,i1,i2,P11:
option remember:
if i=0 then 
RETURN(1):
fi:
P:=a*b*add(add(c[i1,i2]*a^i1*b^i2,i2=0..2*i-2-i1),i1=0..2*i-2):
 var:={seq(seq(c[i1,i2],i2=0..2*i-2-i1),i1=0..2*i-2)}:


P0:=f2iab(i-1,a,b):

P1:=expand(FL(P,[a,b])+FL(P0,[a,b])-P0):


eq:={}:


for i1 from 0 to degree(P1,a) do
P11:=coeff(P1,a,i1):
for i2 from 0 to degree(P11,b) do
eq:=eq union {coeff(P11,b,i2)}:
od:
od:
eq:=eq minus {0}:

var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:


P:=subs(var,P):
factor(P):

end:



#M2iab(i,a,b): The i-th moment of the duration of a 2-player gambler's ruin with starting capitals a,b. Try:
#M2iab(4,a,b):
M2iab:=proc(i,a,b) local i1: option remember:factor(expand(add(stirling2(i,i1)*f2iab(i1,a,b)*i1!,i1=1..i))):end:


#M2iabT(i,a,b): The i-th moment of the duration of a 2-player gambler's ruin with TIES, with starting capitals a,b. Try:
#M2iabT(4,a,b):
M2iabT:=proc(i,a,b) local i1: option remember:factor(expand(add(stirling2(i,i1)*f2iabT(i1,a,b)*i1!,i1=1..i))):end:


#AllMoms2(K,a,b): The list of the first K (straignt) moments
AllMoms2:=proc(K,a,b) local i: option remember: [seq(M2iab(i,a,b),i=1..K)]:end:


#AllMoms2T(K,a,b): The list of the first K (straignt) moments of gambler's ruin with TIES
AllMoms2T:=proc(K,a,b) local i: option remember: [seq(M2iabT(i,a,b),i=1..K)]:end:


#AllMoms2am(K,a,b): The list of first K moments about the meam of the r.v. duration in a fair gambler's ruin with starting capitals a,b
#AllMoms2am(4,a,b);
AllMoms2am:=proc(K,a,b) local lu,i,mu,gu,gu1,j:

lu:=AllMoms2(K,a,b):

mu:=lu[1]:

gu:=[mu]:

for i from 2 to K do
 gu1:=factor(add(lu[i-j]*(-mu)^j*binomial(i,j),j=0..i-1)+(-mu)^i):
 gu:=[op(gu),gu1]:
od:

gu:
end:


#AllMoms2amT(K,a,b): The list of first K moments about the meam of the r.v. duration in a fair gambler's ruin wtih TIES, with starting capitals a,b
#AllMoms2amT(4,a,b);
AllMoms2amT:=proc(K,a,b) local lu,i,mu,gu,gu1,j:

lu:=AllMoms2T(K,a,b):

mu:=lu[1]:

gu:=[mu]:

for i from 2 to K do
 gu1:=factor(add(lu[i-j]*(-mu)^j*binomial(i,j),j=0..i-1)+(-mu)^i):
 gu:=[op(gu),gu1]:
od:

gu:
end:


#AllMoms2amS(K,a,b): The list of first K SCALED moments about the meam of the r.v. duration in a fair gambler's ruin with starting capitals a,b
AllMoms2amS:=proc(K,a,b) local gu,i:
gu:=AllMoms2am(K,a,b):
factor([gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..nops(gu))]):
end:


#AllMoms2amST(K,a,b): The list of first K SCALED moments about the meam of the r.v. duration in a fair gambler's ruin with TIES with starting capitals a,b
AllMoms2amST:=proc(K,a,b) local gu,i:
gu:=AllMoms2amT(K,a,b):
factor([gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..nops(gu))]):
end:


#AllMoms2L(K,x1,x2): The list whose first entry is the exact expectation of the duration of a 2-player fair gambler's ruin with starting capital N*x1, N*x2 divided by N^2
#(where of course x1+x2=1)
#the second entry is the limit of the variance divided by N^4, and the remaining ones are the limits of the scaled moments as N goes to infinity for the 3rd up to the K-th.
#Try:
#AllMoms2L(4,x1,x2);
AllMoms2L:=proc(K,x1,x2) local gu,a,b,N,i:
gu:=AllMoms2amS(K,a,b):
[normal(subs({a=x1*N,b=x2*N},gu[1])/N^2), limit(normal(subs({a=x1*N,b=x2*N},gu[2])/N^4),N=infinity),seq(limit(subs({a=x1*N,b=x2*N},gu[i]),N=infinity),i=3..nops(gu))]:
end:




#Paper2(R,a1,a2): outputs a paper about the first R moments of the duration of 2-player gambler's ruin until one of them goes broke. Try:
#Paper2(6,A,B);
Paper2:=proc(R,a1,a2)  local gu,guAM, guS,guL,x1,x2,i,t0,gv:

t0:=time():

gu:=AllMoms2(R,a1,a2):

guAM:=AllMoms2am(R,a1,a2):

guS:=AllMoms2amS(R,a1,a2):
guL:=AllMoms2L(R,x1,x2):


print(`The First`, R, `moments of the Duration of a 2-Player Gambler's Ruin until one of them goes Broke`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose that 2 gamblers start with initial capitals of`, a1,a2, `dollars. At each round one of the players, picked uniformly at random, has to give one dollar`):
print(`to the other player. The game continues until one of the players goes broke. Let X be duration of this game`):
print(``):
print(`In this paper we will rederive, from scracth, the classical  expressions for expectation, and  variance`):
print(``):
print(`and state rigorously proved explicit expressions for the third all the way to the`, R, `-th moment. In particular we will give expressions for the skewness and kurtosis.`):
print(``):
print(``):
print(`The Expectated Duration is:`):
print(``):
print(gu[1]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu[1]):
print(``):
print(`Note that if N is large and the starting capitals are x1*N,x2*N  (where, of course, x1+x2=1),, then the expectation is N^2 times `):
print(``):
print(guL[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guL[1]):
print(``):

print(``):
print(`The Second Moment is:`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu[2]):
print(``):
print(`hence the variance (aka second moment about the mean), is: `):

print(``):
print(guAM[2]):
print(``):

print(`and in Maple notation`):
print(``):
lprint(guAM[2]):
print(``):

print(`Note that if N is large and the starting capitals are x1*N,x2*N,x3*N, (where, of course, x1+x2+x3=1), then the limit of the variance divided by N^4, as N goes to infinity, is: `):
print(``):
print(guL[2]):
print(``):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guL[2]):
print(``):


for i from 3 to R do
print(``):
print(` The `, i, ` -th (straight) moment of the Duration is: `):
print(``):
print(gu[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[i]):
print(``):
print(`Hence the`, i, ` -th moment about the mean is `):

print(``):
print(guAM[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guAM[i]):
print(``):

print(`Hence the scaled`, i, `-th moment about the mean is`):


if i=3 then
 print(`aka  skewness`):
fi:

if i=4 then
 print(`aka kurtosis`):
fi:

print(``):
print(guS[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guS[i]):
print(``):


print(`Note that if N is large and the starting capitals are x1*N,x2*N,x3*N, (where, of course, x1+x2+x3=1), then the limit of the`, i, `-th  scaled moment about the mean  tends to `):
print(``):
print(guL[i]):
print(``):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guL[i]):
print(``):
od:


guS:=subs(a2=a1,guS):


print(`Finally, let's look what happens when both players start out with the same capital`, a1, `dollars each`):
print(``):
print(`The expected duration is`):

print(``):
print(guS[1]):
print(``):

print(`and in Maple notation`):
print(``):
lprint(guS[1]):
print(``):


print(``):
print(`The variance is`):

print(``):
print(guS[2]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(guS[2]):
print(``):

if R>=3 then

print(``):
print(`The skewness is`):

print(``):
print(guS[3]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(guS[3]):
print(``):
print(`Note that, as`, a1, `goes to infinity it tends to`):
gv:=limit(guS[3],a1=infinity):
lprint(gv):
print(``):
lprint(`that, in decimals is`, evalf(gv)):
fi:




if R>=4 then

print(``):
print(`The kurtosis is`):

print(``):
print(guS[4]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(guS[4]):
print(``):
print(`Note that, as`, a1, `goes to infinity it tends to`):
gv:=limit(guS[4],a1=infinity):
lprint(gv):
print(``):
lprint(`that, in decimals is`, evalf(gv)):
fi:

for i from 5 to R do


print(``):
print(`The scaled `, i, `-th moment about the mean is`):

print(``):
print(guS[i]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(guS[i]):
print(``):
print(`Note that, as`, a1, `goes to infinity it tends to`):
gv:=limit(guS[i],a1=infinity):
lprint(gv):
print(``):
lprint(`that, in decimals is`, evalf(gv)):

od:


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

end:

##start three-players fair Gambler's Ruin
ez3:=proc(): print(` FL(f,L), Paper3(R,a1,a2,a3) `): end:



#f3iabcOld(i,a,b,c): The explicit expression for the i-th binomial moment of the duration of a 3-player Gambler's ruin with initial capitals a,b,c. Try:
#f3iabcOld(2,a,b,c);
f3iabcOld:=proc(i,a,b,c) local eq,var,P,C,P1,P11,P111,P0,i1,i2,i3:
option remember:
if i=0 then 
RETURN(1):
fi:
P:=a*b*c*add(add(add(C[i1,i2,i3]*a^i1*b^i2*c^i3,i3=0..2*i-2-i1-i2), i2=0..2*i-2-i1),i1=0..2*i-2)/(a+b+c):
var:={seq(seq(seq(C[i1,i2,i3],i3=0..2*i-2-i1-i2), i2=0..2*i-2-i1),i1=0..2*i-2)}:

P0:=f3iabc(i-1,a,b,c):

P1:=normal(
P
-1/6*subs({a=a+1,b=b-1},P)-1/6*subs({a=a-1,b=b+1},P)-1/6*subs({a=a+1,c=c-1},P)-1/6*subs({a=a-1,c=c+1},P)-1/6*subs({b=b+1,c=c-1},P)-1/6*subs({b=b-1,c=c+1},P)
-1/6*subs({a=a+1,b=b-1},P0)-1/6*subs({a=a-1,b=b+1},P0)-1/6*subs({a=a+1,c=c-1},P0)-1/6*subs({a=a-1,c=c+1},P0)-1/6*subs({b=b+1,c=c-1},P0)-1/6*subs({b=b-1,c=c+1},P0)
):

#P1:=expand(FL(P,[a,b])+FL(P0,[a,b])-P0):

P1:=expand(numer(P1)):

eq:={}:


for i1 from 0 to degree(P1,a) do
P11:=coeff(P1,a,i1):
for i2 from 0 to degree(P11,b) do
P111:=coeff(P11,b,i2):
for i3 from 0 to degree(P111,c) do
eq:=eq union {coeff(P111,c,i3)}:
od:
od:
od:

eq:=eq minus {0}:

var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:


P:=subs(var,P):
factor(P):

end:



#f3iabc(i,a,b,c): The explicit expression for the i-th binomial moment of the duration of a 3-player Gambler's ruin with initial capitals a,b,c. Try:
#f3iabc(2,a,b,c);
f3iabc:=proc(i,a,b,c) local eq,var,P,C,P1,P11,P111,P0,i1,i2,i3:
option remember:
if i=0 then 
RETURN(1):
fi:
P:=a*b*c*add(add(add(C[i1,i2,i3]*a^i1*b^i2*c^i3,i3=0..2*i-2-i1-i2), i2=0..2*i-2-i1),i1=0..2*i-2)/(a+b+c):
var:={seq(seq(seq(C[i1,i2,i3],i3=0..2*i-2-i1-i2), i2=0..2*i-2-i1),i1=0..2*i-2)}:

P0:=f3iabc(i-1,a,b,c):
P1:=expand(FL(P,[a,b,c])+FL(P0,[a,b,c])-P0):

P1:=expand(numer(P1)):

eq:={}:


for i1 from 0 to degree(P1,a) do
P11:=coeff(P1,a,i1):
for i2 from 0 to degree(P11,b) do
P111:=coeff(P11,b,i2):
for i3 from 0 to degree(P111,c) do
eq:=eq union {coeff(P111,c,i3)}:
od:
od:
od:

eq:=eq minus {0}:

var:=solve(eq,var):

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


P:=subs(var,P):
factor(P):

end:



#M3iabc(i,a,b,c): The i-th moment of the duration of a 3-player gambler's ruin with starting capitals a,b,c. Try:
#M3iabc(4,a,b,c):
M3iabc:=proc(i,a,b,c) local i1: option remember:factor(expand(add(stirling2(i,i1)*f3iabc(i1,a,b,c)*i1!,i1=1..i))):end:


#AllMoms3(K,a,b,c): The list of the first K (straignt) moments
AllMoms3:=proc(K,a,b,c) local i: option remember: [seq(M3iabc(i,a,b,c),i=1..K)]:end:



#AllMoms3am(K,a,b,c): The list of first K moments about the meam of the r.v. duration in a fair 3-player gambler's ruin with starting capitals a,b,c. Try:
#AllMoms3am(4,a,b,c);
AllMoms3am:=proc(K,a,b,c) local lu,i,mu,gu,gu1,j:
option remember:
lu:=AllMoms3(K,a,b,c):

mu:=lu[1]:

gu:=[mu]:

for i from 2 to K do
 gu1:=factor(add(lu[i-j]*(-mu)^j*binomial(i,j),j=0..i-1)+(-mu)^i):
 gu:=[op(gu),gu1]:
od:

gu:
end:


#FL(f,L): Applies the fair Gambler's ruin operator with nops(L) players  to the expression f in the variables of L. Try:
#FL(a*b*c,[a,b,c]);
FL:=proc(f,L) local gu,i,j,k:
k:=nops(L): 
gu:=f:
for i from 1 to k do
 for j from i+1 to k do
  gu:=gu-1/(k*(k-1))*subs({L[i]=L[i]+1, L[j]=L[j]-1},f):
  gu:=gu-1/(k*(k-1))*subs({L[i]=L[i]-1, L[j]=L[j]+1},f):
 od:
od:

expand(gu):
end:

#AllMoms3amS(K,a,b,c): The list of first K SCALED moments about the meam of the r.v. duration in a fair 3-player gambler's ruin with starting capitals a,b,c. Try:
AllMoms3amS:=proc(K,a,b,c) local gu,i:
gu:=AllMoms3am(K,a,b,c):
factor([gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..nops(gu))]):
end:


#AllMoms3L(K,x1,x2,x3): The list whose first entry is the exact expectation of the duration of a 3-player fair gambler's ruin with starting capital N*x1, N*x2 N*x3, divided by N^2
#(where of course x1+x2+x3=1)
#the second entry is the limit of the variance divided by N^4, and the remaining ones are the limits of the scaled moments as N goes to infinity for the 3rd up to the K-th.
#Try:
#AllMoms3L(4,x1,x2,x3);
AllMoms3L:=proc(K,x1,x2,x3) local gu,a,b,c,N,i:
gu:=AllMoms3amS(K,a,b,c):
[normal(subs({a=x1*N,b=x2*N,c=x3*N},gu[1])/N^2), limit(normal(subs({a=x1*N,b=x2*N,c=x3*N},gu[2])/N^4),N=infinity),seq(limit(subs({a=x1*N,b=x2*N,c=x3*N},gu[i]),N=infinity),i=3..nops(gu))]:
end:





#f3iabcT2(i,a,b,c): The explicit expression for the i-th binomial moment of the duration of a 3-player Gambler's ruin with initial capitals a,b,c. Try:
#f3iabcT2(2,a,b,c);
f3iabcT2:=proc(i,a,b,c) local eq,var,P,C,P1,P11,P111,P0,i1,i2,i3,ka1:
option remember:
if i=0 then 
RETURN(1):
fi:
P:=add(add(add(C[i1,i2,i3]*a^i1*b^i2*c^i3,i3=0..2*i-i1-i2), i2=0..2*i-i1),i1=0..2*i):


var:={seq(seq(seq(C[i1,i2,i3],i3=0..2*i-i1-i2), i2=0..2*i-i1),i1=0..2*i)}:


P0:=f3iabcT2(i-1,a,b,c):
P1:=expand(FL(P,[a,b,c])+FL(P0,[a,b,c])-P0):

P1:=expand(P1):

eq:={}:


for i1 from 0 to degree(P1,a) do
P11:=coeff(P1,a,i1):
for i2 from 0 to degree(P11,b) do
P111:=coeff(P11,b,i2):
for i3 from 0 to degree(P111,c) do
eq:=eq union {coeff(P111,c,i3)}:
od:
od:
od:

eq:=eq minus {0}:


ka1:=subs(a=0,P)-f2iab(i,b,c):

for i1 from 0 to degree(ka1,b) do
P11:=coeff(ka1,b,i1):
for i2 from 0 to degree(P11,c) do
P111:=coeff(P11,c,i2):
eq:=eq union {P111}:
od:
od:



ka1:=subs(b=0,P)-f2iab(i,a,c):

for i1 from 0 to degree(ka1,a) do
P11:=coeff(ka1,a,i1):
for i2 from 0 to degree(P11,c) do
P111:=coeff(P11,c,i2):
eq:=eq union {P111}:
od:
od:



ka1:=subs(c=0,P)-f2iab(i,a,b):

for i1 from 0 to degree(ka1,a) do
P11:=coeff(ka1,a,i1):
for i2 from 0 to degree(P11,b) do
P111:=coeff(P11,b,i2):
eq:=eq union {P111}:
od:
od:




var:=solve(eq,var):

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


P:=subs(var,P):

factor(P):


end:

ezraSt:=proc()
if args=NULL then

print(`The Story procedures are: Paper2, Paper2terse, Paper2terseT, Paper3, Paper3terse `):

else
ezra(args):

fi:

end:


#Paper3(R,a1,a2,a3): outputs a paper about the first R moments of the duration of a fair 3-player gambler's ruin, until one of them goes broke. Try:
#Paper3(6,A,B,C);
Paper3:=proc(R,a1,a2,a3)  local gu,guAM, guS,guL,x1,x2,x3,i,t0,gv:

t0:=time():

gu:=AllMoms3(R,a1,a2,a3):

guAM:=AllMoms3am(R,a1,a2,a3):

guS:=AllMoms3amS(R,a1,a2,a3):
guL:=AllMoms3L(R,x1,x2,x3):


print(`The First`, R, `moments of the Duration of a 3-Player Gambler's Ruin until one of them goes Broke`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Suppose that 3 gamblers start with initial capitals of`, a1,a2,a3, `dollars. At each round one of the players, picked uniformly at random, has to give one dollar`):
print(`to one of the other players (also uniformly at random). The game continues until one of the players goes broke. Let X be duration of this game`):
print(``):
print(`In this paper we will rederive, from scracth, Engel's expression for the expectation, and Bruss et. al.'s and Strizaker expression for the variance`):
print(``):
print(`and state rigorously proved explicit expressions for the third all the way to the`, R, `-th moment. In particular we will give expressions for the skewness and kurtosis.`):
print(``):
print(``):
print(`The Expectated Duration is:`):
print(``):
print(gu[1]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu[1]):
print(``):
print(`Note that if N is large and the starting capitals are x1*N,x2*N,x3*N  (where, of course, x1+x2+x3=1),, then the expectation is N^2 times `):
print(``):
print(guL[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guL[1]):
print(``):

print(``):
print(`The Second Moment is:`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple notation `):
print(``):
lprint(gu[2]):
print(``):
print(`hence the variance (aka second moment about the mean), is: `):

print(``):
print(guAM[2]):
print(``):

print(`and in Maple notation`):
print(``):
lprint(guAM[2]):
print(``):

print(`Note that if N is large and the starting capitals are x1*N,x2*N,x3*N, (where, of course, x1+x2+x3=1), then the limit of the variance divided by N^4, as N goes to infinity, is: `):
print(``):
print(guL[2]):
print(``):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guL[2]):
print(``):


for i from 3 to R do
print(``):
print(` The `, i, ` -th (straight) moment of the Duration is: `):
print(``):
print(gu[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[i]):
print(``):
print(`Hence the`, i, ` -th moment about the mean is `):

print(``):
print(guAM[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guAM[i]):
print(``):

print(`Hence the scaled`, i, `-th moment about the mean is`):


if i=3 then
 print(`aka  skewness`):
fi:

if i=4 then
 print(`aka kurtosis`):
fi:

print(``):
print(guS[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guS[i]):
print(``):


print(`Note that if N is large and the starting capitals are x1*N,x2*N,x3*N, (where, of course, x1+x2+x3=1), then the limit of the`, i, `-th  scaled moment about the mean  tends to `):
print(``):
print(guL[i]):
print(``):
print(``):
print(`and in Maple notation`):
print(``):
lprint(guL[i]):
print(``):
od:

guS:=subs({a2=a1,a3=a1},guS):


print(`Finally, let's look what happens when all three  players start out with the same capital`, a1, `dollars each`):
print(``):
print(`The expected duration is`):

print(``):
print(guS[1]):
print(``):

print(`and in Maple notation`):
print(``):
lprint(guS[1]):
print(``):


print(``):
print(`The variance is`):

print(``):
print(guS[2]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(guS[2]):
print(``):

if R>=3 then

print(``):
print(`The skewness is`):

print(``):
print(guS[3]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(guS[3]):
print(``):
print(`Note that, as`, a1, `goes to infinity it tends to`):
gv:=limit(guS[3],a1=infinity):
lprint(gv):
print(``):
lprint(`that, in decimals is`, evalf(gv)):
fi:




if R>=4 then

print(``):
print(`The kurtosis is`):

print(``):
print(guS[4]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(guS[4]):
print(``):
print(`Note that, as`, a1, `goes to infinity it tends to`):
gv:=limit(guS[4],a1=infinity):
lprint(gv):
print(``):
lprint(`that, in decimals is`, evalf(gv)):
fi:

for i from 5 to R do


print(``):
print(`The scaled `, i, `-th moment about the mean is`):

print(``):
print(guS[i]):
print(``):


print(`and in Maple notation`):
print(``):
lprint(guS[i]):
print(``):
print(`Note that, as`, a1, `goes to infinity it tends to`):
gv:=limit(guS[i],a1=infinity):
lprint(gv):
print(``):
lprint(`that, in decimals is`, evalf(gv)):

od:



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

end:



#AllMoms2Le(K): Inputs a positive integer K>=2, and outputs a list of NUMBERS of length K such that,
#in a fair 2-player gambler's ruin, where both players have initial capital a dollars
# the first entry is the expected duration divided by a^2
#the first entry is the limit of the variance of the duration divided by a^4, as a goes to infinity
#The third through the K-th entry are the limits of the scaled moments, as a goes to infinity. Try:
#AllMoms2Le(10);
AllMoms2Le:=proc(K) local a, b,gu,i:
gu:=AllMoms2amS(K,a,b):
gu:=normal(subs(b=a,gu)):
[gu[1]/a^2, limit(gu[2]/a^4,a=infinity), seq(limit(gu[i],a=infinity),i=3..K)]:
end:



#AllMoms3Le(K): Inputs a positive integer K>=2, and outputs a list of NUMBERS of length K such that,
#in a fair 3-player gambler's ruin, where all three players have initial capital a dollars
# the first entry is the expected duration divided by a^2
#the first entry is the limit of the variance of the duration divided by a^4, as a goes to infinity
#The third through the K-th entry are the limits of the scaled moments, as a goes to infinity. Try:
#AllMoms3Le(10);
AllMoms3Le:=proc(K) local a, b,c,gu,i:
gu:=AllMoms3amS(K,a,b,c):
gu:=normal(subs({b=a,c=a},gu)):
[gu[1]/a^2, limit(gu[2]/a^4,a=infinity), seq(limit(gu[i],a=infinity),i=3..K)]:
end:





####start with ties




#f2iabT(i,a,b): The explicit expression for the i-th binomial moment of the duration of a 2-player Gambler's ruin with initial capitals a,b, where ties are allowed. Try:
#f2iabT(2,a,b);
f2iabT:=proc(i,a,b) local eq,var,P,c,P1,P0,i1,i2,P11:
option remember:
if i=0 then 
RETURN(1):
fi:
P:=a*b*add(add(c[i1,i2]*a^i1*b^i2,i2=0..2*i-2-i1),i1=0..2*i-2):
 var:={seq(seq(c[i1,i2],i2=0..2*i-2-i1),i1=0..2*i-2)}:


P0:=f2iabT(i-1,a,b):

P1:=expand(1/2*P-1/4*subs({a=a-1,b=b+1},P)-1/4*subs({a=a+1,b=b-1},P)    +1/2*P0-1/4*subs({a=a-1,b=b+1},P0)-1/4*subs({a=a+1,b=b-1},P0)-P0):

eq:={}:


for i1 from 0 to degree(P1,a) do
P11:=coeff(P1,a,i1):
for i2 from 0 to degree(P11,b) do
eq:=eq union {coeff(P11,b,i2)}:
od:
od:
eq:=eq minus {0}:

var:=solve(eq,var):
if var=NULL then
 RETURN(FAIL):
fi:


P:=subs(var,P):
factor(P):

end:





#Paper2terse(R,a1,a2): outputs a terse paper about the first R moments (binomial, then straight, then about-the-mean, then sacled) of the duration of 2-player gambler's ruin until one of them goes broke. Try:
#Paper2terse(6,A,B);
Paper2terse:=proc(R,a1,a2)  local guB, gu,guAM, guS,i,a:


guB:=[seq(f2iab(i,a1,a2),i=1..R)]:

gu:=AllMoms2(R,a1,a2):

guAM:=AllMoms2am(R,a1,a2):

guS:=AllMoms2amS(R,a1,a2):


print(`The first`, R, `binomial  moments of the duration of a fair gambler's ruin game with 2 players with initial capitals`, a1,a2, `are `):
print(``):
lprint(guB):
print(``):


print(`The first`, R, `straight  moments, i.e. E[D^i], of the duration of a fair gambler's ruin game with 2 players with initial capital`, a1,a2, `are `):
print(``):
lprint(gu):
print(``):


print(`The first`, R, ` moments-about the mean, i.e. E[(D-mu)^i], where mu=E[D], of the duration of a fair gambler's ruin game with 2 players with initial capital`, a1,a2, `are `):
print(``):
lprint(guAM):
print(``):


print(`The first`, R, ` Scaled moments-about the mean, i.e. E[(D-mu)^i]/Var(D)^(i/2), where mu=E[D], of the duration of a fair gambler's ruin game with 2 players with initial capital`, a1,a2, `are `):
print(``):
lprint(guS):
print(``):
guS:=subs({a1=a,a2=a},guS):
guS:=[guS[1]/a^2,limit(guS[2]/a^4,a=infinity), seq(limit(guS[i],a=infinity),i=3..nops(gu))]:
print(``):
print(`If the initial capitals are equal, both a, then the first moment divided by a^2, the limit of the variance divided by a^4, as a goes to infinity, and the limits of the scaled moments as a goes to infinity up to the`, R, `are `):
print(``):
lprint(guS):
print(``):
print(`and in floating points`):
print(``):
lprint(evalf(guS)):
print(``):

print(`----------------------------------------`):
print(``):
print(`This took `, time(), `seconds . `):
print(``):
end:


#Paper2terseT(R,a1,a2): outputs a terse paper about the first R moments (binomial, then straight, then about-the-mean, then sacled) of the duration of 2-player gambler's ruin with TIES, until one of them goes broke. Try:
#Paper2terseT(6,A,B);
Paper2terseT:=proc(R,a1,a2)  local guB, gu,guAM, guS,i,a:


guB:=[seq(f2iabT(i,a1,a2),i=1..R)]:

gu:=AllMoms2T(R,a1,a2):

guAM:=AllMoms2amT(R,a1,a2):

guS:=AllMoms2amST(R,a1,a2):


print(`The first`, R, `binomial  moments of the duration of a fair gambler's ruin game with 2 players with initial capitals`, a1,a2, `are `):
print(``):
lprint(guB):
print(``):


print(`The first`, R, `straight  moments, i.e. E[D^i], of the duration of a fair gambler's ruin game with 2 players with initial capital`, a1,a2, `are `):
print(``):
lprint(gu):
print(``):


print(`The first`, R, ` moments-about the mean, i.e. E[(D-mu)^i], where mu=E[D], of the duration of a fair gambler's ruin game with 2 players with initial capital`, a1,a2, `are `):
print(``):
lprint(guAM):
print(``):


print(`The first`, R, ` Scaled moments-about the mean, i.e. E[(D-mu)^i]/Var(D)^(i/2), where mu=E[D], of the duration of a fair gambler's ruin game with 2 players with initial capital`, a1,a2, `are `):
print(``):
lprint(guS):
print(``):
guS:=subs({a1=a,a2=a},guS):
guS:=[guS[1]/a^2,limit(guS[2]/a^4,a=infinity), seq(limit(guS[i],a=infinity),i=3..nops(gu))]:
print(``):
print(`If the initial capitals are equal, both a, then the first moment divided by a^2, the limit of the variance divided by a^4, as a goes to infinity, and the limits of the scaled moments as a goes to infinity up to the`, R, `are `):
print(``):
lprint(guS):
print(``):
print(`and in floating points`):
print(``):
lprint(evalf(guS)):
print(``):

print(`----------------------------------------`):
print(``):
print(`This took `, time(), `seconds . `):
print(``):
end:






#Paper3terse(R,a1,a2,a3): outputs a terse paper about the first R moments (binomial, then straight, then about-the-mean, then sacled) of the duration of 2-player gambler's ruin until one of them goes broke. Try:
#Paper3terse(6,A,B,C);
Paper3terse:=proc(R,a1,a2,a3)  local guB, gu,guAM, guS,i,a:


guB:=[seq(f3iabc(i,a1,a2,a3),i=1..R)]:

gu:=AllMoms3(R,a1,a2,a3):

guAM:=AllMoms3am(R,a1,a2,a3):

guS:=AllMoms3amS(R,a1,a2,a3):


print(`The first`, R, `binomial  moments of the duration of a fair gambler's ruin game with 3 players with initial capitals`, a1,a2,a3,  `are `):
print(``):
lprint(guB):
print(``):


print(`The first`, R, `straight  moments, i.e. E[D^i], of the duration of a fair gambler's ruin game with 3 players with initial capitals`, a1,a2, a3, `are `):
print(``):
lprint(gu):
print(``):


print(`The first`, R, ` moments-about the mean, i.e. E[(D-mu)^i], where mu=E[D], of the duration of a fair gambler's ruin game with 3 players with initial capital`, a1,a2, a3, `are `):
print(``):
lprint(guAM):
print(``):


print(`The first`, R, ` Scaled moments-about the mean, i.e. E[(D-mu)^i]/Var(D)^(i/2), where mu=E[D], of the duration of a fair gambler's ruin game with 2 players with initial capitals`, a1,a2, a3, `are `):
print(``):
lprint(guS):
print(``):
guS:=subs({a1=a,a2=a,a3=a},guS):
guS:=[guS[1]/a^2,limit(guS[2]/a^4,a=infinity), seq(limit(guS[i],a=infinity),i=3..nops(gu))]:
print(``):
print(`If the initial capitals are equal, both a, then the first moment divided by a^2, the limit of the variance divided by a^4, as a goes to infinity, and the limits of the scaled moments as a goes to infinity up to the`, R, `are `):
print(``):
lprint(guS):
print(``):
print(`and in floating points`):
print(``):
lprint(evalf(guS)):
print(``):

print(`----------------------------------------`):
print(``):
print(`This took `, time(), `seconds . `):
print(``):
end:







#f3iabcTie(i,a,b,c): The explicit expression for the i-th binomial moment of the duration of a 3-player Gambler's ruin with  TIES (see p. 42 of Jiri Andel's book) with initial capitals a,b,c where ties are allowed. Try:
#f3iabcTie(2,a,b,c);
f3iabcTie:=proc(i,a,b,c) local eq,var,P,C,P1,P0,i1,i2,i3,P11,P111:
option remember:
if i=0 then 
RETURN(1):
fi:
P:=a*b*c*add(add(add(C[i1,i2,i3]*a^i1*b^i2*c^i3,i3=0..2*i-2-i1-i2),i2=0..2*i-2-i1),i1=0..2*i-2)/(a+b+c-2):
 var:={seq(seq(seq(C[i1,i2,i3],i3=0..2*i-2-i1-i2),i2=0..2*i-2-i1),i1=0..2*i-2)}:


P0:=f3iabcTie(i-1,a,b,c):

P1:=3/4*P-1/4*subs({a=a+2,b=b-1,c=c-1},P)-1/4*subs({b=b+2,a=a-1,c=c-1},P)-1/4*subs({c=c+2,a=a-1,b=b-1},P)+3/4*P0-1/4*subs({a=a+2,b=b-1,c=c-1},P0)-1/4*subs({b=b+2,a=a-1,c=c-1},P0)-1/4*subs({c=c+2,a=a-1,b=b-1},P0)-P0:
P1:=expand(numer(P1)):

eq:={}:


for i1 from 0 to degree(P1,a) do
P11:=coeff(P1,a,i1):
for i2 from 0 to degree(P11,b) do
P111:=coeff(P11,b,i2):
for i3 from 0 to degree(P111,c) do
eq:=eq union {coeff(P111,c,i3)}:
od:
od:
od:

eq:=eq minus {0}:


var:=solve(eq,var):
lprint(`var is`, var):
if var=NULL then
 RETURN(FAIL):
fi:


P:=subs(var,P):
factor(P):

end: