######################################################################
##PARRONDO: Save this file as PARRONDO. To use it, stay in the       #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read PARRONDO : <Enter>                            #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Temple University ,                   #
#zeilberg@math.temple.edu.                                           # 
#######################################################################
 
#Created: Oct. 12, 2000
#This version: Oct. 12, 2000
#PARRONDO: A Maple package to study the PARRONDO PARADOX
#It accompanies Ekhad and Zeilberger's article
#"Remarks on the Parrondo Paradox"
#Please report bugs to zeilberg@math.temple.edu
 
 
print(`PARRONDO: A Maple package to study the PARRONDO PARADOX`):
print(`Please report bugs to zeilberg@math.temple.edu`):
print(`It accompanies Ekhad and Zeilberger's article`):
print(`"Remarks on the Parrondo Paradox"`):
print(`Created: Oct. 12, 2000`):
print(`This version: Oct. 12, 2000`):
lprint(``):
print(`Written by Doron Zeilberger, zeilberg@math.temple.edu`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.temple.edu/~zeilberg/`):
print(`For a list of all procedures, type ezra1(); `):
 print(`For a list of the IMPORTANT procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 print(``):
 
ezra1:=proc()
if args=NULL then
 print(`Contains the following procedures: AABBGame, `):
 print(`AABBGame1,ABGameR, ABGameR1, AGame, AGame1,Bestp, BestWord, BGame, `):
 print(` BGame1, EWperM,G2Game, G2Game1, GBestWord, GF,GFE, GFV, `):
 print(` GWordGame, GWordGame1, Juan, Juan1, WordGame`):
 print(` WordGame0, WordGame1 `):
 print(`For help with any of these  type ezra(proc_name); `):
fi:
end:
 
ezra:=proc()
if args=NULL then
 print(` Bestp, GBestWord, GWordGame, GWordGame1, GBestWord`):
fi:
 
if nops([args])=1 and op(1,[args])=AABBGame then
print(`AABBGame(s,e): the generating function, in s,`):
print(`of the expected win of game AABB`):
fi:
 
 
if nops([args])=1 and op(1,[args])=AABBGame1 then
print(`AABBGame1(e): The expected win per move of game AABB`):
print(`normalized by dividing by 4 (since any single move of `):
print(` of AABB coreesponds to 4 moves of A or B `):
fi:
 
if nops([args])=1 and op(1,[args])=ABGameR then
print(`ABGameR(s,p,e): the generating function, in s,`):
print(`of the expected win of random combination of`):
print(`A's and B's with bias p towards A`):
fi:
 
if nops([args])=1 and op(1,[args])=ABGameR1 then
print(`ABGameR1(p,e): the asymptotic win per move`):
print(`of the expected win of random combination of`):
print(`A's and B's with bias p towards A`):
fi:
 
if nops([args])=1 and op(1,[args])=AGame then
print(`AGame(s,e): the generating function, in s,`):
print(`of the expected win of game A`):
fi:
 
if nops([args])=1 and op(1,[args])=AGame1 then
print(`AGame1(e): The expected win per move of game A`):
fi:
 
if nops([args])=1 and op(1,[args])=Bestp then
print(`Bestp(e): for a given bias, e,`):
print(`finds the  p such that in  a random`):
print(`A-B game with probability p of picking`):
print(`A and prob. 1-p of of picking B,`):
print(`then the axpected pay-off per move is`):
print(`maximized, followed by that value`):
print(`For example: try: Bestp(0); and Bestp(1/100);`):
fi:
 
if nops([args])=1 and op(1,[args])=BestWord then
print(`BestWord(n,A,B,e): best word of length n with its`):
print(`accompanying value, for bias e`):
fi:
 
if nops([args])=1 and op(1,[args])=BGame then
print(`BGame(s,e): the generating function, in s,`):
print(`of the expected win of game B`):
fi:
 
if nops([args])=1 and op(1,[args])=BGame1 then
print(`BGame1(e): The expected win per move of game B`):
fi:
 
if nops([args])=1 and op(1,[args])=EWperM then
print(`EWperM(Pmat,t): the asymptotic expected win per move`):
print(`in a game of chance governed by Pmat (in terms of t)`):
fi:
 
 
if nops([args])=1 and op(1,[args])=G2Game then
print(`G2Game(p,q,s): the g.f. (in s) for expected win per move`):
print(`in a game of chance governed by a two-state game of chance using `):
print(`letters p and q`):
fi:
 
if nops([args])=1 and op(1,[args])=G2Game1 then
print(`G2Game1(p,q): the asymptotic expected win per move`):
print(`in a game of chance governed by a two-state game of chance using `):
print(`letters p and q`):
fi:
 
 
if nops([args])=1 and op(1,[args])=GBestWord then
print(`GBestWord(n,A,B,M,pr2,e): best word of length n with its`):
print(`accompanying value, for bias e (must be a rational number)`):
print(`for generalized Parrondo with M (the usual is 3) and`):
print(` p_2=pr2-e (usual pr3=3/4)`):
print(`For example, try: GBestWord(3,A,B,3,3/4,1/100)`):
fi:
 
 
if nops([args])=1 and op(1,[args])=GF then
print(`GFk(Pmat,s): Given the transition matrix in a Markovian`):
print(`game of chance, Pmat (as a list of lists of polynomials`):
print(`in some variable), finds the rational function in s,`):
print(`the coeff. of s^n of `):
print(`whose Taylor expansion (w.r.t. s=0) gives, `):
print(`the generating polynomial (in t) of the probability`):
print(`distribution of  income after n rounds of the game`):
fi:
 
if nops([args])=1 and op(1,[args])=GFE then
print(`GFE(Pmat,s): Given the transition matrix in a Markovian`):
print(`game of chance, Pmat (as a list of lists of polynomials`):
print(`in some variable), finds the rational function in s,`):
print(`the coeff. of s^n of `):
print(`whose Taylor expansion (w.r.t. s=0) gives, `):
print(`the expected win`):
print(`after n rounds of the game`):
fi:
 
if nops([args])=1 and op(1,[args])=GFV then
print(`GFV(Pmat,t,s): the rational function that in s,`):
print(`when Taylor-expanded w.r.t. s=0, has coeff. s^n`):
print(`that equals to the second moment of the gain after n rounds`):
print(`of the game of chance given by Pmat (using the variable t)`):
fi:
 
if nops([args])=1 and op(1,[args])=GWordGame then
print(`GWordGame(W,A,B,M,pr2,s,e): the generating function, in s,`):
print(`of the expected win of game W[1]W[2]..W[nops(W)]`):
print(`for generalized Parrondo A,B given by prob. pr2 (formerly 3/4)`):
print(`and general M (formerly M=3)`):
print(`and the deviation e (when e=0 it is fair)`):
print(` For example, try: GWordGame([A,B,B],A,B,3,3/4,s,0)`):
fi:
 
if nops([args])=1 and op(1,[args])=GWordGame1 then
print(`GWordGame1(W,A,B,M,pr2,e): the asymptotic `):
print(`of the expected win of game per move of  W[1]W[2]..W[nops(W)]`):
print(`for generalized Parrondo A,B given by prob. pr2 (formerly 3/4)`):
print(`and general M (formerly M=3)`):
print(`and the deviation e (when e=0 it is fair)`):
print(` For example, try: GWordGame1([A,B,B],A,B,3,3/4,1/1000)`):
fi:
 
if nops([args])=1 and op(1,[args])=Juan then
print(`Juan(s,p1,p2,M): the generating function, in s,`):
print(`for the expected win after n moves `):
print(`for a generalized game B with probabilty of win`):
print(`p1 if the capital is divisible by M and`):
print(`p2 if it is not divisible by M`):
fi:
 
if nops([args])=1 and op(1,[args])=Juan1 then
print(`Juan1(p1,p2,M): the asyptotic win per move`):
print(`for a generalized game B with probabilty of win`):
print(`p1 if the capital is divisible by M and`):
print(`p2 if it is not divisible by M`):
fi:
 
if nops([args])=1 and op(1,[args])=WordGame then
print(`WordGame(W,A,B,s,e): the generating function, in s,`):
print(`of the expected win of game W[1]W[2]..W[nops(W)]`):
fi:
 
if nops([args])=1 and op(1,[args])=WordGame0 then
print(`WordGame0(W,A,B,s,t,e): the generating function, in s,`):
print(`of the probability distibution of game W[1]W[2]..W[nops(W)]`):
fi:
 
if nops([args])=1 and op(1,[args])=WordGame1 then
print(`WordGame1(W,A,B,e): the per move asymptotics,`):
print(`of the expected win of game W[1]W[2]..W[nops(W)]`):
fi:
 
end:
 
 
#GFk(Pmat,s): Given the transition matrix in a Markovian
#game of chance, Pmat (as a list of lists of polynomials
#in some variable), finds the rational function in s,
#the coeff. of s^n of 
#whose Taylor expansion (w.r.t. s=0) gives, 
#the generating polynomial (in t) of the probability
#distribution of  income after n rounds of the game
GFk:=proc(Pmat,s)
local P1,L1,C1,I1,n,i1:
with(linalg):
n:=nops(Pmat):
P1:=convert(Pmat,matrix):
L1:=matrix(1,n,[1,seq(0,i1=1..n-1)]):
C1:=matrix(n,1,[seq(1,i1=1..n)]):
I1:=array(identity,1..n,1..n):
normal(multiply(multiply(L1,inverse(evalm(I1-s*P1))),C1)[1,1]):
end:
 
 
#GFki(Pmat): Given the transition matrix in a Markovian
#game of chance, Pmat (as a list of lists of polynomials
#in some variable), finds the generating polynomial
#of all outcomes
GFki:=proc(Pmat)
local P1,L1,C1,n,i1:
with(linalg):
n:=nops(Pmat):
P1:=convert(Pmat,matrix):
L1:=matrix(1,n,[1,seq(0,i1=1..n-1)]):
C1:=matrix(n,1,[seq(1,i1=1..n)]):
normal(multiply(L1,multiply(P1,C1))[1,1]):
end:
 
#GFE(Pmat,t,s): the rational function that in s,
#when Taylor-expanded w.r.t. s=0, has coeff. s^n
#that equals to the expeced gain after n rounds
#of the game of chance given by Pmat (using the variable t)
GFE:=proc(Pmat,t,s):
normal(subs(t=1,diff(GFk(Pmat,s),t))):
end:
 
 
#GFV(Pmat,t,s): the rational function that in s,
#when Taylor-expanded w.r.t. s=0, has coeff. s^n
#that equals to the second moment of the gain after n rounds
#of the game of chance given by Pmat (using the variable t)
GFV:=proc(Pmat,t,s):
normal(subs(t=1,diff(diff(t*GFk(Pmat,s),t),t))):
end:
 
#AGame(s,e): the generating function, in s,
#of the expected win of game A
AGame:=proc(s,e)
local A,t:
A:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
GFE(A,t,s):
end:
 
 
#AGame0(s,t,e): the generating function, in s,
#of the probabilty distribution in t, of win of game A
AGame0:=proc(s,t,e)
local A:
A:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
GFk(A,s):
end:
 
 
#AGame1(e): The expected win per move of game A
AGame1:=proc(e)
local A,t:
A:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
EWperM(A,t):
end:
 
#BGame(s,e): the generating function, in s,
#of the expected win of game B
BGame:=proc(s,e)
local B,t:
B:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
GFE(B,t,s):
end:
 
#BGame0(s,t,e): the generating function, in s,
#of the expected win of game B
BGame0:=proc(s,t,e)
local B:
B:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
GFk(B,s):
end:
 
#Juan(s,p1,p2,M): the generating function, in s,
#for a generalized game B with probabilty of win
#p1 if the capital is divisible by M and
#p2 if it is not divisible by M
Juan:=proc(s,p1,p2,M)
local lu,i1,gu,i,t:
lu:=[0,p1*t,seq(0,i1=1..M-3),(1-p1)/t]:
gu:=[lu]:
 
for i from 2 to M-1 do
lu:=[seq(0,i1=1..i-2),(1-p2)/t,0,p2*t,seq(0,i1=1..M-i-1)]:
gu:=[op(gu),lu]:
od:
lu:=[p2*t,seq(0,i1=1..M-3),(1-p2)/t,0]:
 
gu:=[op(gu),lu]:
GFE(gu,t,s):
end:
 
#Juan1(p1,p2,M): the expected win
#for a generalized game B with probabilty of win
#p1 if the capital is divisible by M and
#p2 if it is not divisible by M
Juan1:=proc(p1,p2,M)
local t,lu,i1,gu,i:
lu:=[0,p1*t,seq(0,i1=1..M-3),(1-p1)/t]:
gu:=[lu]:
 
for i from 2 to M-1 do
lu:=[seq(0,i1=1..i-2),(1-p2)/t,0,p2*t,seq(0,i1=1..M-i-1)]:
gu:=[op(gu),lu]:
od:
lu:=[p2*t,seq(0,i1=1..M-3),(1-p2)/t,0]:
gu:=[op(gu),lu]:
EWperM(gu,t):
end:
 
 
#BGame1(s,e): The expected win per move of game B
BGame1:=proc(e)
local B,t:
B:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
EWperM(B,t):
end:
 
#AABBGame(s,e): the generating function, in s,
#of the expected win of game AABB
AABBGame:=proc(s,e)
local A,B,C,t,i1,j1:
with(linalg):
A:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
B:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
A:=convert(A,matrix):
B:=convert(B,matrix):
C:=multiply(A,multiply(A,multiply(B,B))):
C:=[ seq([seq(C[i1,j1],j1=1..3)],i1=1..3)]:
GFE(C,t,s):
end:
 
 
#AABBGame1(e): the expected win per move of game AABB
#normalized by dividing by 4
AABBGame1:=proc(e)
local A,B,C,t,i1,j1:
with(linalg):
A:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
B:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
A:=convert(A,matrix):
B:=convert(B,matrix):
C:=multiply(A,multiply(A,multiply(B,B))):
C:=[ seq([seq(C[i1,j1],j1=1..3)],i1=1..3)]:
normal(EWperM(C,t)/4):
end:
 
 
#ABGameR(s,p,e): the generating function, in s,
#of the expected win of random combination of
#A's and B's with bias p towards A
ABGameR:=proc(s,p,e)
local A,B,R,t,i1,j1:
with(linalg):
A:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
B:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
A:=convert(A,matrix):
B:=convert(B,matrix):
R:=evalm(p*A+(1-p)*B):
R:=[ seq([seq(R[i1,j1],j1=1..3)],i1=1..3)]:
GFE(R,t,s):
end:
 
 
#ABGameR1(p,e): the asymptotic expected win, per move
#of a truly random combination of
#A's and B's with bias p towards A
ABGameR1:=proc(p,e)
local A,B,R,t,i1,j1:
with(linalg):
A:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
B:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
A:=convert(A,matrix):
B:=convert(B,matrix):
R:=evalm(p*A+(1-p)*B):
R:=[ seq([seq(R[i1,j1],j1=1..3)],i1=1..3)]:
normal(EWperM(R,t)):
end:
 
#EWperM(Pmat,t): the asymptotic expected win per move
#in a game of chance governed by Pmat (in terms of t)
EWperM:=proc(Pmat,t)
local s,gu:
gu:=GFE(Pmat,t,s):
subs(s=1,normal((s-1)^2*gu)):
end:
 
#EVperM(Pmat,t): the asymptotic variance win per move
#in a game of chance governed by Pmat (in terms of t)
EVperM:=proc(Pmat,t)
local s,gu:
gu:=GFV(Pmat,t,s):
subs(s=1,normal((s-1)^3*gu))/2-EWperM(Pmat,t)^2:
end:
 
 
#G2Game(p,q,s): the generetic function, in s of
#expected win per move
#in a game of chance governed by a two-state game of chance using 
#letters p and q
G2Game:=proc(p,q,s)
local A,t:
A:=
[[p[1,1]*t+q[1,1]/t,p[1,2]*t+(1-p[1,1]-q[1,1]-p[1,2])/t],
[p[2,1]*t+q[2,1]/t,p[2,2]*t+(1-p[2,1]-q[2,1]-p[2,2])/t]]: 
GFE(A,t,s)
end:
 
#G2Game1(p,q): the asymptotic expected win per move
#in a game of chance governed by a two-state game of chance using 
#letters p and q
G2Game1:=proc(p,q)
local A,t:
A:=
[[p[1,1]*t+q[1,1]/t,p[1,2]*t+(1-p[1,1]-q[1,1]-p[1,2])/t],
[p[2,1]*t+q[2,1]/t,p[2,2]*t+(1-p[2,1]-q[2,1]-p[2,2])/t]]: 
EWperM(A,t)
end:
 
#Bestp(e): for a given `bias', e,
#finds the  p such that in  a random
#A-B game with probability p of picking
#A and prob. 1-p of of picking B,
#then the axpected pay-off per move is
#maximized, followed by that value
 
Bestp:=proc(e)
local e1,gu,mu,p,hakhi,zokhe,i,lu,muam:
 
gu:=ABGameR1(p,e1):
gu:=normal(gu):
mu:=diff(gu,p):
mu:=numer(mu):
mu:=subs(e1=e,mu):
gu:=subs(e1=e,gu):
lu:=[fsolve(mu,p)]:
 
zokhe:=0:
hakhi:=subs(p=0,gu):
muam:=subs(p=1,gu):
 
if muam>hakhi then
hakhi=muam:
zokhe:=1:
fi:
 
for i from 1 to nops(lu) do
if lu[i]>0 and lu[i]<1 then
muam:=subs(p=lu[i],gu):
 
if muam>hakhi then
hakhi:=muam:
zokhe:=lu[i]:
fi:
fi:
od:
 
hakhi,zokhe:
 
end:
 
 
#WordGame(W,A,B,s,e): the generating function, in s,
#of the expected win of game W[1]W[2]..W[nops(W)]
WordGame:=proc(W,A,B,s,e)
local A1,B1,C,t,i1,j1,i:
with(linalg):
A1:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
B1:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
 
A1:=convert(A1,matrix):
B1:=convert(B1,matrix):
 
 
if W[1]=A then
 C:=A1:
else
 C:=B1:
fi:
 
for i from 2 to nops(W) do
 
if W[i]=A then
 C:=multiply(C,A1):
else
 C:=multiply(C,B1):
fi:
od:
 
C:=[ seq([seq(C[i1,j1],j1=1..3)],i1=1..3)]:
GFE(C,t,s):
end:
 
 
#WordGame0(W,A,B,s,t,e): the generating function, in s,
#of the expected win of game W[1]W[2]..W[nops(W)]
WordGame0:=proc(W,A,B,s,t,e)
local A1,B1,C,i1,j1,i:
with(linalg):
A1:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
B1:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
 
A1:=convert(A1,matrix):
B1:=convert(B1,matrix):
 
 
if W[1]=A then
 C:=A1:
else
 C:=B1:
fi:
 
for i from 2 to nops(W) do
 
if W[i]=A then
 C:=multiply(C,A1):
else
 C:=multiply(C,B1):
fi:
od:
 
C:=[ seq([seq(C[i1,j1],j1=1..3)],i1=1..3)]:
GFk(C,s):
end:
 
 
#WordGame1(W,A,B,e): the asymptotic expected
#win per move of game W[1]W[2]..W[nops(W)]
#where W is a word in A and B given in terms of a list
WordGame1:=proc(W,A,B,e)
local A1,B1,C,t,i1,j1,i:
with(linalg):
A1:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
B1:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
 
A1:=convert(A1,matrix):
B1:=convert(B1,matrix):
 
 
if W[1]=A then
 C:=A1:
else
 C:=B1:
fi:
 
for i from 2 to nops(W) do
 
if W[i]=A then
 C:=multiply(C,A1):
else
 C:=multiply(C,B1):
fi:
od:
 
C:=[ seq([seq(C[i1,j1],j1=1..3)],i1=1..3)]:
normal(EWperM(C,t)/nops(W)):
end:
 
 
#WordGamei(W,A,B,e): the expected
#win per move of playing W[1]W[2]..W[nops(W)]
#only ONCE
#where W is a word in A and B given in terms of a list
WordGamei:=proc(W,A,B,e)
local A1,B1,C,t,i1,j1,i,gu:
with(linalg):
A1:=
[[0,(1/2-e)*t,(1/2+e)/t],[(1/2+e)/t,0,(1/2-e)*t],[(1/2-e)*t,(1/2+e)/t,0]]:
B1:=[[0,(1/10-e)*t,(9/10+e)/t],[(1/4+e)/t,0,(3/4-e)*t],
[(3/4-e)*t,(1/4+e)/t,0]]:
 
A1:=convert(A1,matrix):
B1:=convert(B1,matrix):
 
 
if W[1]=A then
 C:=A1:
else
 C:=B1:
fi:
 
for i from 2 to nops(W) do
 
if W[i]=A then
 C:=multiply(C,A1):
else
 C:=multiply(C,B1):
fi:
od:
 
C:=[ seq([seq(C[i1,j1],j1=1..3)],i1=1..3)]:
 
gu:=GFki(C):
subs(t=1,diff(gu,t))/nops(W):
end:
 
#Words(A,B,n): The set of all words of length n
#in A and B
Words:=proc(A,B,n)
local gu,mu,i:
if n=1 then
RETURN({[A],[B]}):
fi:
mu:=Words(A,B,n-1):
 
gu:={}:
for i from 1 to nops(mu) do
gu:=gu union {[op(mu[i]),A],[op(mu[i]),B]}:
od:
gu:
end:
 
 
#BestWord(n,A,B,e): best word of length n with its
#accompanying value, for bias e
BestWord:=proc(n,A,B,e)
local gu,i,hakhi,muam,zokhe:
gu:=Words(A,B,n):
 
zokhe:={gu[1]}:
hakhi:=WordGame1(gu[1],A,B,e):
 
for i from 2 to nops(gu) do
muam:=WordGame1(gu[i],A,B,e):
 
if muam>hakhi then
 hakhi:=muam:
 zokhe:={gu[i]}:
fi:
if muam=hakhi then
 zokhe:=zokhe union {gu[i]}:
fi:
 
od:
 
zokhe,hakhi:
 
end:
 
 
#BestWordi(n,A,B,e): best word of length n (for playing ONCE)
#with its accompanying value, for bias e
BestWordi:=proc(n,A,B,e)
local gu,i,hakhi,muam,zokhe:
gu:=Words(A,B,n):
 
zokhe:={gu[1]}:
hakhi:=WordGamei(gu[1],A,B,e):
 
for i from 2 to nops(gu) do
muam:=WordGamei(gu[i],A,B,e):
 
if muam>hakhi then
 hakhi:=muam:
 zokhe:={gu[i]}:
fi:
if muam=hakhi then
 zokhe:=zokhe union {gu[i]}:
fi:
 
od:
 
zokhe,hakhi:
 
end:
 
 
 
#GWordGame(W,A,B,M,pr2,s,e): the generating function, in s,
#of the expected win of game W[1]W[2]..W[nops(W)]
#for generalized Parrondo A,B given by prob. pr2 (formerly 3/4)
#and general M (formerly M=3)
#and the deviation e (when e=0 it is fair)
GWordGame:=proc(W,A,B,M,pr2,s,e)
local A1,B1,C,t,i1,j1,i,lu,p1,p2,mu,GU,gu:
with(linalg):
 
p1:=solve(Juan1(p1,p2,M),p1):
p1:=subs(p2=pr2,p1):
p2:=pr2:
p1:=p1-e:
p2:=p2-e:
 
mu:=[0,(1/2-e)*t,seq(0,i1=1..M-3),(1/2+e)/t]:
lu:=[0,p1*t,seq(0,i1=1..M-3),(1-p1)/t]:
gu:=[lu]:
GU:=[mu]:
for i from 2 to M-1 do
lu:=[seq(0,i1=1..i-2),(1-p2)/t,0,p2*t,seq(0,i1=1..M-i-1)]:
mu:=[seq(0,i1=1..i-2),(1/2+e)/t,0,(1/2-e)*t,seq(0,i1=1..M-i-1)]:
gu:=[op(gu),lu]:
GU:=[op(GU),mu]:
od:
lu:=[p2*t,seq(0,i1=1..M-3),(1-p2)/t,0]:
gu:=[op(gu),lu]:
mu:=[(1/2+e)*t,seq(0,i1=1..M-3),(1/2-e)/t,0]:
GU:=[op(GU),mu]:
 
A1:=convert(GU,matrix):
B1:=convert(gu,matrix):
 
 
if W[1]=A then
 C:=A1:
else
 C:=B1:
fi:
 
for i from 2 to nops(W) do
 
if W[i]=A then
 C:=multiply(C,A1):
else
 C:=multiply(C,B1):
fi:
od:
 
C:=[ seq([seq(C[i1,j1],j1=1..3)],i1=1..3)]:
GFE(C,t,s):
end:
 
 
#GWordGame1(W,A,B,M,pr2,e): the expected win per move
#of the expected win of game W[1]W[2]..W[nops(W)]
#for generalized Parrondo A,B given by prob. pr2 (formerly 3/4)
#and general M (formerly M=3)
#and the deviation e (when e=0 it is fair)
GWordGame1:=proc(W,A,B,M,pr2,e)
local A1,B1,C,t,i1,j1,i,lu,p1,p2,mu,GU,gu:
with(linalg):
 
p1:=solve(Juan1(p1,p2,M),p1):
p1:=subs(p2=pr2,p1):
p2:=pr2:
p1:=p1-e:
p2:=p2-e:
 
mu:=[0,(1/2-e)*t,seq(0,i1=1..M-3),(1/2+e)/t]:
lu:=[0,p1*t,seq(0,i1=1..M-3),(1-p1)/t]:
gu:=[lu]:
GU:=[mu]:
for i from 2 to M-1 do
lu:=[seq(0,i1=1..i-2),(1-p2)/t,0,p2*t,seq(0,i1=1..M-i-1)]:
mu:=[seq(0,i1=1..i-2),(1/2+e)/t,0,(1/2-e)*t,seq(0,i1=1..M-i-1)]:
gu:=[op(gu),lu]:
GU:=[op(GU),mu]:
od:
lu:=[p2*t,seq(0,i1=1..M-3),(1-p2)/t,0]:
gu:=[op(gu),lu]:
mu:=[(1/2+e)*t,seq(0,i1=1..M-3),(1/2-e)/t,0]:
GU:=[op(GU),mu]:
 
A1:=convert(GU,matrix):
B1:=convert(gu,matrix):
 
 
if W[1]=A then
 C:=A1:
else
 C:=B1:
fi:
 
for i from 2 to nops(W) do
 
if W[i]=A then
 C:=multiply(C,A1):
else
 C:=multiply(C,B1):
fi:
od:
 
C:=[ seq([seq(C[i1,j1],j1=1..3)],i1=1..3)]:
normal(EWperM(C,t)/nops(W)):
end:
 
 
 
#GBestWord(n,A,B,M,pr2,e): best word of length n with its
#accompanying value, for Parrondo M and p_2=pr2 for bias e
GBestWord:=proc(n,A,B,M,pr2,e)
local gu,i,hakhi,muam,zokhe:
gu:=Words(A,B,n):
 
zokhe:={gu[1]}:
 
hakhi:=GWordGame1(gu[1],A,B,M,pr2,e):
 
for i from 2 to nops(gu) do
muam:=GWordGame1(gu[i],A,B,M,pr2,e):
 
if muam>hakhi then
 hakhi:=muam:
 zokhe:={gu[i]}:
fi:
if muam=hakhi then
 zokhe:=zokhe union {gu[i]}:
fi:
 
od:
 
zokhe,hakhi:
 
end: