# restart; read `c:/Users/Asus/Google Drive/Aek # /Pile Game/Piles.txt`; # read `c:/Users/Asus/Google Drive/Guess.txt`; # First version: Jan 18,2019 with(LinearAlgebra): Digits :=500: ####################################### # Maple program accompanied the paper # Probabilistic Two Piles Game by # Ho Hon Leung and Thotsaporn Aek # Thanatipanonda # ####################################### # Section 0: Path counting # Section 1: Introduction with case 0 0 then S := S union {seq([op(P),a], P in RecPath(p+a,s-1,t,n,a,b))}; fi: if p+b < n and t >0 then S := S union {seq([op(P),b], P in RecPath(p+b,s,t-1,n,a,b))}; fi: S; end: # The number of paths which arises from # take a or b chips with prob 1/2 # with total of k turns, # never get to n chips, at the end # of k turns has exactly v chips. # Try: # NumPath(1,-1,2,18,0); # nops(Paths(1,-1,2,18,0)); NumPath := proc(n,a,b,k,v) option remember; local s,t,M; M := solve({s+t=k, a*s+b*t=v}); s := subs(M,s); t := subs(M,t); if type(s,nonnegint) = false or type(t,nonnegint) = false then return(0); fi: CountPath(0,s,t,n,a,b); end: # At position p # The number of paths # which arises from # taking a chips s times # and b chips t times and # never get equal or more # than n chips # Try: # CountPath(0,3,3,1,-1,1); # nops(RecPath(0,3,3,1,-1,1)); CountPath := proc(p,s,t,n,a,b) option remember; local P,S; if s=0 and t=0 then return(1); fi: S := 0; if p+a < n and s > 0 then S := S + CountPath(p+a,s-1,t,n,a,b); fi: if p+b < n and t >0 then S := S + CountPath(p+b,s,t-1,n,a,b); fi: S; end: ################################# # Section 1: Case 0 a or a >= b then ERROR("a must non-neg and b must be greater than a"); elif k <0 then return(0); fi: m := (n-a*k)/(b-a); if type(m,nonnegint) = false then return(0); else return(binomial(k,m)); fi: end: # The probability of winning at turn k, # R1(n,k) for 0 < a < b; # Try: # seq(print([seq(R1(n,k,1,2),k=1..n)]),n=13..15); # d := (n,k) ->1/2^k*(binomial(k,n-k)+binomial(k-1,n-k)): # seq(print([seq(d(n,k),k=1..n)]),n=13..15); R1 := proc(n,k,a,b) option remember; add(ForPath(n-i,k-1,a,b)/2^(k-1),i=1..a) +add(ForPath(n-i,k-1,a,b)/2^k,i=a+1..b); end: # Generate Q1 by equation (2) # Try: # n:=39: # [seq( [add( binomial(k,i),i=0..n-k-1)/2^k, # Q1(n,k,1,2)], k=0..32)]; # n:=39: a:=2: b:=3: # [seq( [add( binomial(k,i),i=0..(n-a*k-1)/(b-a))/2^k, # Q1(n,k,a,b)], k=0..32)]; Q1 := proc(n,k,a,b) option remember; local j; if n<0 then ERROR(NoCase); elif n=0 then return(0); fi: 1-add(R1(n,j,a,b),j=1..k); end: # Probability of winning # for second player. Check # the formula with definition # Try: # Prob1(230,1,2); Prob1 := proc(n,a,b) option remember; local k,L,R; L := 1/2-1/2*add(R1(n,k,a,b)^2,k=1..n); R := add( Q1(n,k,a,b)*R1(n,k,a,b) ,k=1..n); [L, simplify(L-R)]; end: # Test sum of r^2(n,k) when a=1, b=2 # Try: A:=[seq(R1Sq12(n),n=3..40)]; # Guess1D(A,5,1,3,n,N); R1Sq12 := proc(n) option remember; local k,r; r := k-> (binomial(k,n-k)+binomial(k-1,n-k))/2^k; add( r(k)^2 ,k=floor(n/2)..n); end: ################################# # Section 2: Case a=-1, b=1 # # Section 2.1, 2.2: Winning Prob. # Functions: # Cat(n,k), R2(n,k), # Q2(n,k), SumR2Sq(n), # TRec1(n), Prob2All(n,K), # ################################# # Generalize Catalan numbers # Number of ways to play k-turns # with -1 and 1 move and ended up # exactly at n-1 (never get to n). # Try: # n:=5:[seq(NumPath(n,-1,1,k,n-1),k=0..14)]; # [seq(Cat(n,k),k=0..14)]; Cat := proc(n,k) option remember; local m,s; if n mod 2 = k mod 2 then return(0); elif n mod 2 =1 then m := k/2; s := (n-1)/2; return((2*s+1)*binomial(2*m,m-s)/(m+s+1)); elif n mod 2 =0 then m := (k+1)/2; s := n/2; return(s*binomial(2*m,m-s)/m); fi: end: # R(n,k) for (a,b)= (-1,1) # Try: # [seq(R2(3,k),k=0..14)]; # [seq(R2(2,k),k=1..14)]; R2 := proc(n,k) option remember; Cat(n,k-1)/2^k; end: # Q(n,k) for (a,b)= (-1,1) # Try: # n:=3: [seq( [add(NumPath(n,-1,1,k,i),i=-k..n-1)/2^k, # Q2(n,k)], k=0..32)]; Q2 := proc(n,k) option remember; local j; if n<0 then ERROR(NoCase); elif n=0 then return(0); fi; 1-add(R2(n,j),j=1..k); end: # sum(r(n,k)^2,k=1..infinity); # Try: # n:=4: A:=SumR2Sq(n); evalf(A); # add(evalf(R2(n,k)^2),k=1..5000); # A:=[seq(SumR2Sq(n),n=1..20)]: # Guess1D(A,3,1,1,n,N); # B:=[seq(SumR2Sq(2*n),n=0..35)]: # Guess1D(B,3,4,0,n,N); SumR2Sq := proc(n) option remember; local m,s,SS; if n mod 2 =1 then s := (n-1)/2; SS := ((2*s+1)/2)^2 *sum( (binomial(2*m,m-s)/(m+s+1)/4^m)^2 ,m=0..infinity); elif n mod 2 =0 then s := n/2; SS := s^2*sum( (binomial(2*m,m-s)/m/4^m)^2 ,m=1..infinity); fi: return(SS); end: # Compute sum(r(n,k)^2,k=1..infinity); # from the recurrence # Try: # TRec1(7); TRec1 := proc(n) option remember; local A; if n =1 then return((4-Pi)/Pi); elif n =2 then return((16-5*Pi)/Pi); elif n =3 then return((236-75*Pi)/(3*Pi)); fi: A := ((7*n-5)*TRec1(n-1) -(7*n-16)*TRec1(n-2) +(n-3)*TRec1(n-3))/n; simplify(A); end: # Check the formula of # the winning prob. with # definition # Try: Prob2All(3,3000); Prob2All := proc(n,K) option remember; local L; L := 1/2-1/2*SumR2Sq(n); [simplify(L),evalf(L-Prob2K(n,K))]; end: ################################# # Extra Section to prove formally # the guessing recurrences on # section 2.2 # # Functions: # MyTEven(N), MapleTEven(N), # HomTEven(n,T), # MyTOdd(N), # ################################# # {a,b}={-1,1} # My test on T_(2n) from # my Guessed Recurrence # Try: MyTEven(10); MyTEven := proc(N) option remember; local n; [seq( simplify(n*(2*n+1)*(12*n^2+48*n+43)*SumR2Sq(2*n) +(-687-3939*n-6718*n^2-4236*n^3-840*n^4)*SumR2Sq(2*n+2) +(3000+12717*n+13954*n^2+5844*n^3+840*n^4)*SumR2Sq(2*n+4) -(2*n+5)*(n+3)*(12*n^2+24*n+7)*SumR2Sq(2*n+6)), n=1..N)]; end: # {a,b}={-1,1} # My test on T_(2n) from # Maple's guessed recurrence # Try: MapleTEven(10); MapleTEven := proc(N) option remember; local n; [seq( simplify( (48*n^4+168*n^3+178*n^2+53*n)*SumR2Sq(2*n) +(-1632*n^4-6528*n^3-8884*n^2-4712*n-774)*SumR2Sq(2*n+2) +(48*n^4+216*n^3+322*n^2+179*n+30)*SumR2Sq(2*n+4) -( -2/Pi^2/n/(n+1)^2/(2+n) *(192*n^6*Pi+1152*n^5*Pi+2608*n^4*Pi +192*n^4-192*n^4*cos(Pi*(2+n))^2+768*n^3 -768*n^3*cos(Pi*(2+n))^2+2752*n^3*Pi+1144*n^2 -1144*n^2*cos(Pi*(2+n))^2+1328*n^2*Pi+752*n -752*n*cos(Pi*(2+n))^2+224*n*Pi+159 -159*cos(Pi*(2+n))^2) ) ),n=1..N)]; end: # Try to find a simple homogenous # relation from MapleTEven for T(2n) # Note: This is the same now with my guess # MyTEven(10); # Try: HomTEven(n); HomTEven := proc(n,T) option remember; local i,S,C,A,Z,XX; S := (48*n^4+168*n^3+178*n^2+53*n)*T(2*n) +(-1632*n^4-6528*n^3-8884*n^2-4712*n-774)*T(2*n+2) +(48*n^4+216*n^3+322*n^2+179*n+30)*T(2*n+4) -( -2/Pi^2/n/(n+1)^2/(2+n) *(192*n^6*Pi+1152*n^5*Pi+2608*n^4*Pi +192*n^4-192*n^4*cos(Pi*(2+n))^2+768*n^3 -768*n^3*cos(Pi*(2+n))^2+2752*n^3*Pi+1144*n^2 -1144*n^2*cos(Pi*(2+n))^2+1328*n^2*Pi+752*n -752*n*cos(Pi*(2+n))^2+224*n*Pi+159 -159*cos(Pi*(2+n))^2) ): C:=( -2/Pi^2/n/(n+1)^2/(2+n) *(192*n^6*Pi+1152*n^5*Pi+2608*n^4*Pi +192*n^4-192*n^4*cos(Pi*(2+n))^2+768*n^3 -768*n^3*cos(Pi*(2+n))^2+2752*n^3*Pi+1144*n^2 -1144*n^2*cos(Pi*(2+n))^2+1328*n^2*Pi+752*n -752*n*cos(Pi*(2+n))^2+224*n*Pi+159 -159*cos(Pi*(2+n))^2) ) ; A :=simplify(C*subs(n=n+1,S)-subs(n=n+1,C)*S); XX := 0; for i from 0 to 3 do Z := coeff(A,T(2*(n+i)),1); XX := XX+factor(subs( {cos(Pi*(2+n))^2=1, cos(Pi*(3+n))^2=1},Z))*T(2*(n+i)); od: XX := factor(XX*Pi/(24*n^2+72*n+53)/32); add(factor(coeff(XX,T(2*(n+i)),1))*T(2*(n+i)),i=0..3); end: # {a,b}={-1,1} # My test on T_(2n+1) # from my guessed recurrence # Note: It could be done # automatically from Maple command # as well. # Try: MyTOdd(10); MyTOdd := proc(N) option remember; local n; [seq( simplify( (2*n+1)*(n+1)*(6*n^2+30*n+35)*SumR2Sq(2*n+1) +(-2459-7127*n-7166*n^2-2958*n^3-420*n^4)*SumR2Sq(2*n+3) +(6815+15737*n+11990*n^2+3762*n^3+420*n^4)*SumR2Sq(2*n+5) -(2*n+7)*(3+n)*(6*n^2+18*n+11)*SumR2Sq(2*n+7) ),n=0..N)]; end: ############################# # Section 2.3: Winning Prob. # within K turns # # Functions: # Prob2K(n,K), For1Prob2K(k), # ############################# # Probability of second player # to win within K turn # Try: [seq(Prob2K(1,K),K=1..20)]; # [seq(ForProb2K(K),K=1..20)]; Prob2K := proc(n,K) option remember; local j,k,A,B; A:=add( Q2(n,k)*R2(n,k) ,k=1..K); return(A); if n mod 2 =1 then B := add( Q2(n,2*j-1)*R2(n,2*j-1) ,j=1..floor((K+1)/2)); else B := add( Q2(n,2*j)*R2(n,2*j) ,j=1..floor(K/2)); fi: return(B); end: # Closed form formula for Prob2K # when n=1, player play k turns # Try: [seq(For1Prob2K(k),k=1..20)]; # [seq(Prob2K(1,k),k=1..20)]; For1Prob2K := proc(k) option remember; local L; L := floor((k+1)/2); 1- (2*L+1)/16^L*binomial(2*L,L)^2; end: ######################### # Section 2.4: Average # duration of Game Play # # Functions: # AvgGame2(n), # ######################### # Fix L as in the paper # Maple evaluate the infinite # sum for n=1,2,3,4 only. # Try: AvgGame2(n); AvgGame2 := proc(n) option remember; local s,L,q; if n mod 2 =1 then s := (n-1)/2; q := 1-(2*s+1)/2*sum( binomial(2*m,m-s)/(m+s+1)/4^m ,m=0..L-1); elif n mod 2 =0 then s := n/2; q := 1-s*sum( binomial(2*m,m-s)/m/4^m ,m=1..L); fi: q := simplify(q); print(q); sum(q^2,L=0..infinity); end: ################################# # Section 3: Case a=-1, b=2 # # Functions: # RecD(n,k), R3(n,k), # Q3(n,k), SumR3Sq(n), # Prob3(n,K), # ################################# # Generalize 3-Raney numbers: # The number of ways to play k-turns # with -1 and 2 move and ended up # exactly at n-1 (never get to n). # Try: # n:=15:[seq(NumPath(n,-1,2,k,n-1) # +NumPath(n,-1,2,k,n-2),k=0..34)]; # [seq( RecD(n,k),k=0..34)]; RecD := proc(n,k) option remember; local m; if n <= 0 then return(0); elif n = 1 and k mod 3 = 0 then m := k/3; return( binomial(3*m,m)/(2*m+1)); elif n = 1 and k mod 3 = 1 then m := (k-1)/3; return( binomial(3*m+1,m+1)/(2*m+1)); fi: RecD(n-1,k+1)-RecD(n-3,k); end: # R(n,k) for (a,b)= (-1,2) # Try: # [seq(R3(2,k),k=1..20)]; # [seq(R3(3,k),k=1..20)]; R3 := proc(n,k) option remember; RecD(n,k-1)/2^k; end: # Q(n,k) for (a,b)= (-1,2) # Try: # n:=3: [seq( [add(NumPath(n,-1,2,k,i),i=-k..n-1)/2^k, # Q3(n,k)], k=0..32)]; Q3 := proc(n,k) option remember; local j; if n<0 then ERROR(NoCase); elif n=0 then return(0); fi; 1-add(R3(n,j),j=1..k); end: # sum(R(n,k)^2,k=1..infinity); # No closed form of these sums # The numerical value seems correct. # Try: # n:=1: A:=SumR3Sq(n); evalf(A); # add(evalf(R3(n,k))^2,k=1..7000); SumR3Sq := proc(n) option remember; local m,SS; if n =1 then SS := sum( (binomial(3*m,2*m+1)/m/8^m/2)^2 +((3*m+1)!/(2*m+1)!/(m+1)!/8^m/4)^2 ,m=0..infinity); elif n =2 then SS := sum( (binomial(3*m,2*m+1)/m/8^m)^2 ,m=1..infinity) +sum( ((3*m+1)!/(2*m+1)!/(m+1)!/8^m/2)^2 ,m=0..infinity); elif n =6 then SS := sum( ((3*m)!/(2*m+3)!/m!*3 *(m-1)*(5*m+4)/8^m*2)^2 +((3*m+1)!/(2*m+3)!/(m+2)! *3*(m+1)*(5*m^2+4*m-4)/8^m)^2 ,m=1..infinity); else ERROR("NoCase"); fi: return(SS); end: # Probability of winning # for second player. Check # the formula with definition # Try: # Prob3(3,5000); Prob3 := proc(n,K) option remember; local k,L,R; L := 1/2-1/2*add(R3(n,k)^2,k=1..K); R := add( Q3(n,k)*R3(n,k) ,k=1..K); [evalf(L),evalf(R),evalf(L-R)]; end: ################################# # Section 4: # Winning with n1 and n2 chips # # Functions: # DefWin(n1,n2), Win1(n1,n2), # Prop7(n1,n2,K), OldWin2(n1,n2), # Win2(n1,n2), # ################################# # Winning probability for the # second player (player 1,2 wins # if reaches n1,n2 chips # respectively) # Test with case a=1,b=2! # Try: [seq(DefWin(13,n),n=6..27)]; DefWin := proc(n1,n2) option remember; local k; add( Q1(n1,k,1,2)*R1(n2,k,1,2),k=1..n1); end: # Section 5: page 11 of Wong & Xu # Recursive way to calculate # prob. of Bob wins # (Case a=1, b=2) # Try: # [seq(Win1(n,n),n=1..20)]; # [seq(Prob1(n,1,2)[1],n=1..20)]; # [seq(Win1(13,n),n=6..27)]; # [seq(DefWin(13,n),n=6..27)]; Win1 := proc(n1,n2) option remember; if n1>=1 and n2 <= 0 then return(1); fi: 1-1/2*(Win1(n2,n1-1)+Win1(n2,n1-2)); end: ###################### # Recurrence for case # a=-1, b=1 ###################### # Proposition 7 page 11 # (Generalize of theorem 1) # for case a=-1, b=1 # Try: # Prop7(1,2,5000); Prop7 := proc(n1,n2,K) option remember; local k,A; A := add( R2(n1,k)*R2(n2,k),k=1..K); evalf(Win2(n1,n2)+Win2(n2,n1)+A); end: # Recursive calculation for # case a=-1, b=1. # Prob. that Bob wins # given that Alice starts # From # OldWin2(a,b) := 1-1/2*(OldWin2(b,a+1)+OldWin2(b,a-1)); # Try: # Matrix([seq([seq(OldWin2(a,b),b=1..5)],a=1..5)]); OldWin2 := proc(a,b) option remember; local P; if a=b then P := 1/2-1/2*SumR2Sq(a); elif b = 0 then P := 1; elif a < b then P := 2-2*OldWin2(b-1,a)-OldWin2(a,b-2); elif a > b and a+b mod 2 = 1 then P := 1-OldWin2(b,a); elif a > b and a+b mod 2 = 0 then P := 2-2*OldWin2(b+1,a)-OldWin2(a,b+2); fi: simplify(P); end: # Win2 from definition in section 4 # Try: # Matrix([seq([seq(Win2(n1,n2),n2=1..5)],n1=1..5)]); # [seq([seq(Win2(n1,n2)-OldWin2(n1,n2),n2=1..15)],n1=1..15)]; Win2 := proc(n1,n2) option remember; local j,m,s1,s2,a,Q,R1,R2; if n1 mod 2 =1 and n2 mod 2 =1 then s1 := (n1-1)/2; R1 := m -> (2*s1+1)/(m+s1)/2/4^(m-1) *binomial(2*m-2,m-1-s1); s2 := (n2-1)/2; R2 := m -> (2*s2+1)/(m+s2)/2/4^(m-1) *binomial(2*m-2,m-1-s2); a := 1: elif n1 mod 2 =1 and n2 mod 2 = 0 then s1 := (n1-1)/2; R1 := m -> (2*s1+1)/(m+s1)/2/4^(m-1) *binomial(2*m-2,m-1-s1); s2 := n2/2; R2 := m-> s2/m/4^m*binomial(2*m,m-s2); a := 1: elif n1 mod 2 =0 and n2 mod 2 = 1 then s1 := n1/2; R1 := m-> s1/m/4^m*binomial(2*m,m-s1); s2 := (n2-1)/2; R2 := m -> (2*s2+1)/(m+s2+1)/2/4^(m) *binomial(2*m,m-s2); a := 0: elif n1 mod 2 =0 and n2 mod 2 = 0 then s1 := n1/2; R1 := m-> s1/m/4^m*binomial(2*m,m-s1); s2 := n2/2; R2 := m-> s2/m/4^m*binomial(2*m,m-s2); a := 1: fi: Q := m-> factor(1-sum(R1(j),j=1..m)); sum( Q(m)*R2(m),m=a..infinity); end: