######################################################################
##AliceBob: Save this file as  AliceBob                              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read AliceBob<Enter>                               #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Dec. 15, 2009
 
print(`Created:  Dec. 15, 2009`):
print(` This is AliceBob `):
print(`to study a Coin Tossing game`):
print(`It accompanies the article `):
print(` Optimizing a Losing Game`):
print(`by Stan Wagon, Herbert S. Wilf and Doron Zeilberger`):

print(`and also available from Zeilberger'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 procedures are: `):
 print(` fSeqDirect `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Bestn, BestnSeq,f, FindMamzers,fSeq, Ope `):
 print(` `):


elif nops([args])=1 and op(1,[args])=Bestn then
print(`Bestn(p,q): the best n for Alice in the coin-tossing game`):
print(`with Alice's coin having prob. q of landing Heads`):
print(`and Bob's prob. p of landing Heads (p>q) `):
print(`by using the third-order recurrence outputted by`):
print(`the Apagodu-Zeilberger mutli-variable extensin of the`):
print(`Almkvist-Zeilberger algorithm. `):
print(`Try: `):
print(` Bestn(50/100,49/100); `):

elif nops([args])=1 and op(1,[args])=BestnSeq then
print(`BestnSeq(p,r,K): The sequence [seq(Bestn(p,p-r*i)[1],i=1..K)].`):
print(`For example, try:`):
print(`BestnSeq(1/2,1/100,10);`):

elif nops([args])=1 and op(1,[args])=FindMamzers then
print(`FindMamzers(P,R): finds the pairs (p,q)`):
print(`for p between 1/P and (P-1)/P q=p-1/(2*i) for i between 1 and R`):
print(`such that`):
print(`it is not true that Bestn(p,q)=1/(2*(p-q)). For example, try:`):
print(`FindMamzers(10,10);`):

elif nops([args])=1 and op(1,[args])=f then
print(` f(n,p,q) : the probability that Alice will win a game`):
print(`of throwing loaded coins n times, where Alice's coin has`):
print(`probabilty q of landing Heads, and Bob's coin has probability`):
print(` p of landing Head, computed directly from the formula `):
print(`for example, try:`):
print(` f(10,p,q); `):
print(` f(20,1/2,1/3) `);

elif nops([args])=1 and op(1,[args])=fSeq then
print(`fSeq(K,p,q): the sequence {f(n,p,q)}`):
print(`for n from 1 to K, computed `):
print(`by using the third-order recurrence outputted by`):
print(`the Apagodu-Zeilberger mutli-variable extensin of the`):
print(`Almkvist-Zeilberger algorithm. `):
print(` For example, try: `):
print(` fSeq(10,p,q); `):
print(` fSeq(100,1/2,1/3); `):

elif nops([args])=1 and op(1,[args])=fSeqDirect then
print(`fSeqDirect(K,p,q): the sequence {f(n,p,q)}`):
print(`for n from 1 to K, computed directly, not as fast`):
print(`as fSeq(K,p,q) that uses the recurrence outputted by`):
print(`the Apagodu-Zeilberger mutli-variable extensin of the`):
print(`Almkvist-Zeilberger algorithm. `):
print(`Warning: this is slower than fSeq(K,p,q). Only use for`):
print(`checking purposes`):
print(` For example, try: `):
print(` fSeqDirect(10,p,q); `):
print(` fSeqDirect(100,1/2,1/3); `):

elif nops([args])=1 and op(1,[args])=Ope then
print(`Ope(p,q,n,N): the third-order recurrence operator`):
print(`annihilating f(n,p,q) `):
print(`Do:`):
print(`Ope(p,q,n,N); `):
else
print(`There is no ezra for`,args):
fi:
 
end:


ez:=proc(): 
print(`f(n,p,q), f1(n,p,q), fSeqDirect(K,p,q), f1Seq(K,p,q) `): 
print(` fSeq(K,p,q), BestnD(p,q), Bestn(p,q) `): 
print(`BestnSeq(p,r,K)`):
end:

f:=proc(n,p,q) local r,s:
expand(
add(add(binomial(n,r)*binomial(n,s)*p^r*(1-p)^(n-r)*q^s*(1-q)^(n-s),
s=r+1..n),r=0..n)):
end:


f1:=proc(n,p,q) local x,y,gu,i,j,mu,mu1,mu2:
gu:=(p/x+1-p)^n*(q/y+1-q)^n:
gu:=expand(gu):
mu:=y/(1-x*y)/(1-y):
mu:=taylor(mu,x=0,n+1):
mu1:=0:
for i from 0 to n do
mu2:=coeff(mu,x,i):
mu2:=taylor(mu2,y=0,n+1):
mu2:=add(coeff(mu2,y,j)*y^j,j=0..n):
mu1:=expand(mu1+mu2*x^i):
od:

gu:=expand(gu*mu1):
coeff(coeff(gu,x,0),y,0):
end:


#read MultiAlmkvistZeilberger:
#ope:=MAZ(1,1/x/(1-x*y)/(1-y),(p/x+1-p)*(q/y+1-q),[x,y],n,N,{p,q})[1]:




#SeqFromRec(ope1,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
#For example, try:
#SeqFromRec(N-2,n,N,[1],10);
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
L:=degree(ope,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

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

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), 
expand(normal(-add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..L)))]:
od:

gu:

end:
 
fSeqDirect:=proc(K,p,q) local i:[seq(f(i,p,q),i=1..K)]: end:

f1Seq:=proc(K,p,q) local i:[seq(f1(i,p,q),i=1..K)]: end:

fSeq:=proc(K,p,q) local ope,n,N,Ini,i:
Ini:=[seq(f(i,p,q),i=1..3)]:

ope:=
-(p+q-1)^2*(n+1)*(2*n*p-2*n*q+5*p-5*q-1)/(n+3)/(-2*n*q+2*n*p+3*p-1-3*q)+(-5+6*n
^2*p-13*q+27*p-3*n-6*n*q*p-10*q*p+2*n^2*p^3+8*n^2*q^2+7*n*p^3+29*n*q^2-31*n*p^2
-8*n^2*p^2+35*q*p^2-19*n*q+27*n*p-27*p^2+23*q^2+39*n*q*p^2+10*n^2*q*p^2-6*n^2*q
-7*n*q^3-10*n^2*q^2*p-39*n*q^2*p-35*p*q^2-2*n^2*q^3+5*p^3-5*q^3)/(n+3)/(-2*n*q+
2*n*p+3*p-1-3*q)*N-(4*n^2*q^2-6*n^2*q+6*n^2*p-4*n^2*p^2-8*n^2*q^2*p+8*n^2*q*p^2
-3*n+16*n*q^2+32*n*q*p^2-23*n*q-32*n*q^2*p-16*n*p^2-4*n*q*p+27*n*p-7+29*p-16*p^
2+14*q^2-30*p*q^2-19*q+30*q*p^2-8*q*p)/(n+3)/(-2*n*q+2*n*p+3*p-1-3*q)*N^2+N^3
:
SeqFromRec(ope,n,N,Ini,K):
end:


Ope:=proc(p,q,n,N)
-(p+q-1)^2*(n+1)*(2*n*p-2*n*q+5*p-5*q-1)/(n+3)/(-2*n*q+2*n*p+3*p-1-3*q)+(-5+6*n
^2*p-13*q+27*p-3*n-6*n*q*p-10*q*p+2*n^2*p^3+8*n^2*q^2+7*n*p^3+29*n*q^2-31*n*p^2
-8*n^2*p^2+35*q*p^2-19*n*q+27*n*p-27*p^2+23*q^2+39*n*q*p^2+10*n^2*q*p^2-6*n^2*q
-7*n*q^3-10*n^2*q^2*p-39*n*q^2*p-35*p*q^2-2*n^2*q^3+5*p^3-5*q^3)/(n+3)/(-2*n*q+
2*n*p+3*p-1-3*q)*N-(4*n^2*q^2-6*n^2*q+6*n^2*p-4*n^2*p^2-8*n^2*q^2*p+8*n^2*q*p^2
-3*n+16*n*q^2+32*n*q*p^2-23*n*q-32*n*q^2*p-16*n*p^2-4*n*q*p+27*n*p-7+29*p-16*p^
2+14*q^2-30*p*q^2-19*q+30*q*p^2-8*q*p)/(n+3)/(-2*n*q+2*n*p+3*p-1-3*q)*N^2+N^3
:
end:

#BestnD(p,q): the best n for Alice 
BestnD:=proc(p,q) local n0:

if p<=q then
 RETURN(FAIL):
fi:

for n0 from 1 while f(n0,p,q)<f(n0+1,p,q) do od:
n0,f(n0,p,q):
end:



#Bestn(p,q): the best n for Alice using the recurrence
Bestn:=proc(p,q)
local ope1,gu,n1,n,N,j1,mamzer:

if p<=q or q<=0 or q>=1 or p<=0 or p>=1 then
 RETURN(FAIL):
fi:

mamzer:=(1-3*(p-q))/2/(p-q):

if type(mamzer,integer) and mamzer>0 then
 RETURN(BestnD(p,q)):
fi:

ope1:=
-(p+q-1)^2*(n+1)*(2*n*p-2*n*q+5*p-5*q-1)/(n+3)/(-2*n*q+2*n*p+3*p-1-3*q)+(-5+6*n
^2*p-13*q+27*p-3*n-6*n*q*p-10*q*p+2*n^2*p^3+8*n^2*q^2+7*n*p^3+29*n*q^2-31*n*p^2
-8*n^2*p^2+35*q*p^2-19*n*q+27*n*p-27*p^2+23*q^2+39*n*q*p^2+10*n^2*q*p^2-6*n^2*q
-7*n*q^3-10*n^2*q^2*p-39*n*q^2*p-35*p*q^2-2*n^2*q^3+5*p^3-5*q^3)/(n+3)/(-2*n*q+
2*n*p+3*p-1-3*q)*N-(4*n^2*q^2-6*n^2*q+6*n^2*p-4*n^2*p^2-8*n^2*q^2*p+8*n^2*q*p^2
-3*n+16*n*q^2+32*n*q*p^2-23*n*q-32*n*q^2*p-16*n*p^2-4*n*q*p+27*n*p-7+29*p-16*p^
2+14*q^2-30*p*q^2-19*q+30*q*p^2-8*q*p)/(n+3)/(-2*n*q+2*n*p+3*p-1-3*q)*N^2+N^3
:

ope1:=expand(subs(n=n-3,ope1)/N^3):

gu:=[f(1,p,q),f(2,p,q),f(3,p,q)]:

if f(1,p,q)>f(2,p,q) then
  RETURN(1,f(1,p,1)):
fi:

if f(2,p,q)>f(3,p,q) then
  RETURN(2,f(2,p,2)):
fi:

for n1 from 4 do
gu:=[op(2..3,gu), 
expand(normal(-add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),
j1=1..3)))]:
if gu[3]<gu[2] then
  RETURN(n1-1,gu[2]):
fi:
od:

end:


#BestnSeq(p,r,K): The sequence [seq(Bestn(p,p-r*i)[1],i=1..K)]
#For example, try:
#BestnSeq(1/2,1/100,10);
BestnSeq:=proc(p,r,K) local i:
[seq(Bestn(p,p-r*i)[1],i=1..K)]:
end:

#FindMamzers(P,R): finds the pairs (p,q)
#for p between 1/P and (P-1)/P q=p-1/(2*i) for i between 1 and R
#such that
#it is not true that Bestn(p,q)=1/(2*(p-q)). For example, try:
#FindMamzers(10,10);
FindMamzers:=proc(P,R) local j,i,p,gu:
gu:={}:
for j from 1 to P-1 do
p:=j/P:
for i from trunc(1/2/p)+1 to R do
 if Bestn(p,p-1/(2*i))[1]<>i then
  gu:=gu union {[p,p-1/(2*i)]}:
 fi:
od:
od:

gu:
end: