######################################################################
##Nash.txt: Save this file as  Nash.txt                              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read Nash.txt<Enter>                               #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Feb. 3, 2016
 
print(`Created: Feb. 3,  2016`):
print(` This is Nash.txt `):
print(`to compute Nash Equilibrium`):
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(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the VERBOSE procedures type ezraV();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezraV:=proc()

if args=NULL then
 print(` The Verbose procedures are:  BR1v, BR2v, BRpV, MaxSetV, MNEv, PNEv  `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  IsDominatedR1, IsDominatedR, IsDominatedC1, IsDominatedC  `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: BR1, BR2, BRp, DS, MaxSet, MNE, PNE, PO, RaG `):
 print(` `):



elif nops([args])=1 and op(1,[args])=BR1 then
print(`BR1(M): inputs an m by n matrix of pairs of payoffs, where nops(M)=m, and player 1 has m options`):
print(`and n=nops(M[1]) and player n has n options, and M[i,j] is the pair [PayOffFor1, PayOffFor2]`):
print(`outputs the list of length n containing the sets of best reponses of Player 1 to each of the`):
print(` n strategies of Player 2. `):
print(` Try `):
print(` BR1([[[2,8],[3,5]],[[8,2],[1,7]]]); `):

elif nops([args])=1 and op(1,[args])=BR1v then
print(`BR1v(M): verbose form of BR1(M) `):
print(`inputs an m by n matrix of pairs of payoffs, where nops(M)=m, and player 1 has m options`):
print(`and n=nops(M[1]) and player n has n options, and M[i,j] is the pair [PayOffFor1, PayOffFor2]`):
print(`outputs the list of length n containing the sets of best reponses of Player 1 to each of the`):
print(` n strategies of Player 2. `):
print(` Try `):
print(` BR1v([[[2,8],[3,5]],[[8,2],[1,7]]]); `):

elif nops([args])=1 and op(1,[args])=BR2 then
print(`BR2(M): inputs an m by n matrix of pairs of payoffs, where nops(M)=m, and player 1 has m options`):
print(`and n=nops(M[1]) and player n has n options, and M[i,j] is the pair [PayOffFor1, PayOffFor2]`):
print(`outputs the list of length m containing the sets of best reponses of Player 2 to each of the`):
print(` m strategies of Player 1. `):
print(` Try `):
print(` BR2([[[2,8],[3,5]],[[8,2],[1,7]]]); `):

elif nops([args])=1 and op(1,[args])=BR2v then
print(`BR2v(M): verbose form of BR2(M) `):
print(`BR2v(M): inputs an m by n matrix of pairs of payoffs, where nops(M)=m, and player 1 has m options`):
print(`and n=nops(M[1]) and player n has n options, and M[i,j] is the pair [PayOffFor1, PayOffFor2]`):
print(`outputs the list of length m containing the sets of best reponses of Player 2 to each of the`):
print(` m strategies of Player 1. `):
print(` Try `):
print(` BR2v([[[2,8],[3,5]],[[8,2],[1,7]]]); `):

elif nops([args])=1 and op(1,[args])=BRp then
print(`BRp(f,p,q): inputs a bilinear expression in p,q and outputs the best p for a given q`):
print(`Try:`):
print(` BRp(PO([[[2,8],[3,5]],[[8,2],[1,7]]],p,q)[1],p,q); `):

elif nops([args])=1 and op(1,[args])=BRpV then
print(`BRpV(f,p,q):  verbose form of BRp(f,p,q) `):
print(`BRpV(f,p,q): inputs a bilinear expression in p,q and outputs the best p for a given q`):
print(`Try:`):
print(` BRpV(PO([[[2,8],[3,5]],[[8,2],[1,7]]],p,q)[1],p,q); `):



elif nops([args])=1 and op(1,[args])=DS then
print(`DS(M): inputs a finite 2-player game (given as pay-offs matrix M where M[i][j]=[RowPayOff,ColumnPayOff] if Row played i and Column played j`):
print(`and outputs the pair of sets [DR,DC] where DR is the set of dominated strategies for player Row and DS for player Column.`):
print(` Try: `):
print(` DS([[[2,2],[0,3]],[[3,0],[1,1]]]); `):

elif nops([args])=1 and op(1,[args])=IsDominatedC then
print(`IsDominatedC(M,i,j) :inputs a finite game M, and a strategy i  for Player Column, decides whether i is dominated by some other strategy`):
print(`Try:`):
print(`IsDominatedC([ [[2,2],[0,3]],[[3,0],[1,1]]],1);`):

elif nops([args])=1 and op(1,[args])=IsDominatedC1 then
print(`IsDominatedC1(M,i,j) :inputs a finite game M, and two strategies i and j for Player Column, decides whether i is dominated by j`):
print(`Try:`):
print(`IsDominatedC1([[[2,2],[0,3]],[[3,0],[1,1]]],1,2);`):

elif nops([args])=1 and op(1,[args])=IsDominatedR then
print(`IsDominatedR(M,i,j) :inputs a finite game M, and a strategy i  for Player Row, decides whether i is dominated by some other strategy`):
print(`Try:`):
print(`IsDominatedR([ [[2,2],[0,3]],[[3,0],[1,1]] ],1);`):

elif nops([args])=1 and op(1,[args])=IsDominatedR1 then
print(`IsDominatedR1(M,i,j) :inputs a finite game M, and two strategies i and j for Player Row, decides whether i is dominated by j`):
print(`Try:`):
print(`IsDominatedR1([[[2,2],[0,3]],[[3,0],[1,1]]],1,2);`):

elif nops([args])=1 and op(1,[args])=MaxSet then
print(`MaxSet(L): inputs a list of numbers, L, outputs the set of locations where it is the max.`):
print(`Try:`):
print(`MaxSet([3,4,2,4,1]);`):

elif nops([args])=1 and op(1,[args])=MaxSetV then
print(`MaxSetV(L):inputs a list of numbers, L, outputs the set of locations where it is the max, Verbose version`):
print(`Try:`):
print(`MaxSetV([3,4,2,4,1]);`):

elif nops([args])=1 and op(1,[args])=MNE then
print(`MNE(M): inputs a 2 by 2 matrix M of payoffs and outputs a mixed Nash Equilibrum`):
print(`Try: `):
print(` MNE([[[4,4],[0,3]],[[3,0],[2,2]]]); `):
print(` MNE([[[1,8],[7,9]],[[3,4],[5,5]]]); `):
print(` MNE([[[2,8],[3,5]],[[8,2],[1,7]]]); `):
print(`For Prisoner's Dilemma `):
print(` MNE([[[2,2],[0,3]],[[3,0],[1,1]]]); `):
print(`For BoS `):
print(` MNE([[[2,1],[0,0]],[[0,0],[1,2]]]); `):
print(`For Matching Pennies`):
print(` MNE([[[1,-1],[-1,1]],[[-1,1],[1,-1]]]); `):

elif nops([args])=1 and op(1,[args])=MNEv then
print(`MNEv(M):  verbose form of MNE(M) `):
print(`MNEv(M): inputs a 2 by 2 matrix M of payoffs and outputs a mixed Nash Equilibrum`):
print(`Try: `):
print(` MNEv([[[4,4],[0,3]],[[3,0],[2,2]]]); `):
print(` MNEv([[[1,8],[7,9]],[[3,4],[5,5]]]); `):
print(` MNEv([[[2,8],[3,5]],[[8,2],[1,7]]]); `):
print(`For Prisoner's Dilemma `):
print(` MNEv([[[2,2],[0,3]],[[3,0],[1,1]]]); `):
print(`For BoS `):
print(` MNEv([[[2,1],[0,0]],[[0,0],[1,2]]]); `):
print(`For Matching Pennies`):
print(` MNEv([[[1,-1],[-1,1]],[[-1,1],[1,-1]]]); `):

elif nops([args])=1 and op(1,[args])=PNE then
print(`PNE(M): inputs a matrix M of payoffs and outputs the set of pure Nash Equilibira`):
print(`Try: `):
print(` PNE([[[4,4],[0,3]],[[3,0],[2,2]]]); `):
print(` PNE([[[1,8],[7,9]],[[3,4],[5,5]]]); `):
print(` PNE([[[2,8],[3,5]],[[8,2],[1,7]]]); `):
print(`For Prisoner's Dilemma `):
print(` PNE([[[2,2],[0,3]],[[3,0],[1,1]]]); `):
print(`For BoS `):
print(` PNE([[[2,1],[0,0]],[[0,0],[1,2]]]); `):

elif nops([args])=1 and op(1,[args])=PNEv then
print(`PNEv(M): verbose form of PNE(M) `):
print(`PNEv(M): inputs a matrix M of payoffs and outputs the set of pure Nash Equilibira`):
print(`Try: `):
print(` PNEv([[[4,4],[0,3]],[[3,0],[2,2]]]); `):
print(` PNEv([[[1,8],[7,9]],[[3,4],[5,5]]]); `):
print(` PNEv([[[2,8],[3,5]],[[8,2],[1,7]]]); `):
print(`For Prisoner's Dilemma `):
print(` PNEv([[[2,2],[0,3]],[[3,0],[1,1]]]); `):
print(`For BoS `):
print(` PNEv([[[2,1],[0,0]],[[0,0],[1,2]]]); `):

elif nops([args])=1 and op(1,[args])=PO then
print(`PO(M,p,q): inputs a 2 by 2 game whose entries are [PayOffFor1, PayOffFor 2] outputs the`):
print(` two expected pay-offs if player 1 chooses the first strategy with prob. p `):
print(` and player 2 chooses the first strategy with prob. q. Try: `):
print(` PO([[[2,8],[3,5]],[[8,2],[1,7]]],p,q); `):
print(`For Prisoner's Dilemma `):
print(` PO([[[2,2],[0,3]],[[3,0],[1,1]]],p,q); `):
print(`For BoS `):
print(` PO([[[2,1],[0,0]],[[0,0],[1,2]]],p,q); `):
print(`For Matching Pennies `):
print(` PO([[[1,-1],[-1,1]],[[-1,1],[1,-1]]],p,q); `):

elif nops([args])=1 and op(1,[args])=RaG then
print(`RaG(m,n,K): a random game where Player 1 has m options and Player 2 has n options, and the payoffs are from 1 to K Try:`):
print(`RaG(5,3,10);`):


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


#MaxSet(L): inputs a list of numbers, L, outputs the set of locations where it is the max.
#Try:
#MaxSet([3,4,2,4,1]);
MaxSet:=proc(L) local i,alufim,si:

si:=L[1]:
alufim:={1}:

for i from 2 to nops(L) do
 if L[i]>si then
   alufim:={i}:
   si:=L[i]:
 elif L[i]=si then
      alufim:=alufim union {i}:
fi:
od:

alufim:
end:

  

#BR1(M): inputs an m by n matrix of pairs of payoffs, where nops(M)=m, and player 1 has m options
#and n=nops(M[1]) and player n has n options, and M[i,j] is the pair [PayOffFor1, PayOffFor2]
#outputs the list of length n containing the sets of best reponses if Player 1 to each of the
#n strategies of Player 2.
#Try
#BR1([[[2,8],[3,5]],[[8,2],[1,7]]]);
BR1:=proc(M) local m,n,i,j,gu:
m:=nops(M):
n:=nops(M[1]):
gu:=[]:
for j from 1 to n do
 gu:=[op(gu),MaxSet([seq(M[i][j][1],i=1..m)])]:
od:
gu:
end:

#BR2(M): inputs an m by n matrix of pairs of payoffs, where nops(M)=m, and player 1 has m options
#and n=nops(M[1]) and player n has n options, and M[i,j] is the pair [PayOffFor1, PayOffFor2]
#outputs the list of length m containing the sets of best reponses if Player 2 to each of the
#m strategies of Player 1.
#Try
#BR2([[[2,8],[3,5]],[[8,2],[1,7]]]);
BR2:=proc(M) local m,n,i,j,gu:
m:=nops(M):
n:=nops(M[1]):
gu:=[]:
for i from 1 to m do
 gu:=[op(gu),MaxSet([seq(M[i][j][2],j=1..n)])]:
od:
gu:
end:


#PNE(M): inputs a matrix M of payoffs and outputs the set of pure Nash Equilibira
PNE:=proc(M) local gu1,gu2,m,n, i,j,mu:
m:=nops(M):
n:=nops(M[1]):

gu1:=BR1(M):
gu2:=BR2(M):

mu:={}:


for i from 1 to m do
 for j from 1 to n do
   if member(i,gu1[j]) and  member(j,gu2[i])  then
    mu:=mu  union {[i,j]}:
   fi:
 od:
od:

mu:

end:


#RaG(m,n,K): a random game where Player 1 has m options and Player 2 has n options, and the payoffs are from 1 to K Try:
#RaG(5,3,10);
RaG:=proc(m,n,K) local ra,i,j:
ra:=rand(1..K):

[seq([seq([ra(),ra()],j=1..n)],i=1..m)]:
end:

#PO(M,p,q): inputs a 2 by 2 game whose entries are [PayOffFor1, PayOffFor 2] outputs the
#two expected pay-offs if player 1 chooses the first strategy with prob. p
#and player 2 chooses the first strategy with prob. q. Try:
#PO([[[2,8],[3,5]],[[8,2],[1,7]]],p,q);
PO:=proc(M,p,q):

if not ( nops(M)=2  and nops(M[1])=2 ) then
RETURN(FAIL):
fi:

expand([M[1][1][1]*p*q+M[1][2][1]*p*(1-q)+M[2][1][1]*(1-p)*q+M[2][2][1]*(1-p)*(1-q),
M[1][1][2]*p*q+M[1][2][2]*p*(1-q)+M[2][1][2]*(1-p)*q+M[2][2][2]*(1-p)*(1-q)]):

end:



#BRp(f,p,q): inputs a bilinear expression in p,q and outputs the best p for a given q
#Try:
#BRp(PO([[[2,8],[3,5]],[[8,2],[1,7]]],p,q),p,q);
BRp:=proc(f,p,q) local gu,q0:
gu:=coeff(f,p):

if min(subs(q=0,gu),subs(q=1,gu))>0 then
     RETURN([[1,1],[0,1]]):
fi:

if max(subs(q=0,gu),subs(q=1,gu))<0 then
     RETURN([[0,0],[0,1]]):
fi:

q0:=solve(gu,q):

if subs(q=0,gu)<0 then
[ [[0,0],[0,q0]],[[0,1],[q0,q0]],[[1,1],[q0,1]] ]:
else
[ [[1,1],[0,q0]],[[0,1],[q0,q0]],[[0,0],[q0,1]] ]:
fi:

end:




#MNE(M): inputs a matrix M of payoffs and outputs a Mixed Nash Equilibira

MNE:=proc(M) local p,q, gu,mu1,mu2:


if not ( nops(M)=2  and nops(M[1])=2 ) then
print(`must be a 2 by 2 matrix`):
RETURN(FAIL):
fi:

gu:=PO(M,p,q):

mu1:=BRp(gu[1],p,q):
mu2:=BRp(gu[2],q,p):

if nops(mu1)=3 and nops(mu2)=3 then
 RETURN([mu1[2][2][1],mu2[2][2][1]]):
else
 RETURN(FAIL):
fi:

end:

#IsDominatedR1(M,i,j) :inputs a finite game M, and two strategies i and j for Player Row, decides whether i is dominated by j
#Try:
#IsDominatedR1([[[2,2],[0,3]],[[3,0],[1,1]]],1,2);
IsDominatedR1:=proc(M,i,j)  local m,n,i1:
m:=nops(M):
n:=nops(M[1]):

if not (1<=i and i<=m and 1<=j and j<=m and i<>j) then
 print(`Bad input `):
 RETURN(FAIL):
fi:

if min(seq(M[j][i1][1]-M[i][i1][1],i1=1..n))>0 then
 true:
else
 false:
fi:

end:




#IsDominatedR(M,i) :inputs a finite game M, and a strategy i for Player Row, decides whether i is dominated by  some other strategy
#Try:
#IsDominatedR([[[2,2],[0,3]],[[3,0],[1,1]]],1);
IsDominatedR:=proc(M,i)  local j,m:
m:=nops(M):

for j from 1 to m do
 if j<>i and IsDominatedR1(M,i,j) then
   RETURN(true):
 fi:
od:

false:

end:


#IsDominatedC1(M,i,j) :inputs a finite game M, and two strategies i and j for Player Column, decides whether i is dominated by j
#Try:
#IsDominatedC1([[[2,2],[0,3]],[[3,0],[1,1]]],1,2);
IsDominatedC1:=proc(M,i,j)  local m,n,i1:
m:=nops(M):
n:=nops(M[1]):

if not (1<=i and i<=n and 1<=j and j<=n and i<>j) then
 print(`Bad input `):
 RETURN(FAIL):
fi:

if min(seq(M[i1][j][2]-M[i1][i][2],i1=1..m))>0 then
 true:
else
 false:
fi:

end:




#IsDominatedC(M,i) :inputs a finite game M, and a strategy i for Player Column, decides whether i is dominated by  some other strategy
#Try:
#IsDominatedC([[[2,2],[0,3]],[[3,0],[1,1]]],1);
IsDominatedC:=proc(M,i)  local j,n:
n:=nops(M[1]):

for j from 1 to n do
 if j<>i and IsDominatedC1(M,i,j) then
   RETURN(true):
 fi:
od:

false:

end:



#DS(M): inputs a finite 2-player game (given as pay-offs matrix M where M[i][j]=[RowPayOff,ColumnPayOff] if Row played i and Column played j
#and outputs the pair of sets [DR,DC] where DR is the set of dominated strategies for player Row and DS for player Column.
#Try
#DS([[[2,2],[0,3]],[[3,0],[1,1]]]);
DS:=proc(M) local m,n,DR,DC,i,j:
m:=nops(M):
n:=nops(M[1]):
DR:={}:
DC:={}:

for i from 1 to m do
 if IsDominatedR(M,i) then
     DR:=DR union {i}:
 fi:
od:


for j from 1 to n do
 if IsDominatedC(M,j) then
     DC:=DC union {j}:
 fi:
od:

[DR,DC]:

end:


###start Verbose procedures


#MaxSetV(L): Verbose version of MaxSet(L)
#Try:
#MaxSetV([3,4,2,4,1]);
MaxSetV:=proc(L) local i,alufim,si:

si:=L[1]:
alufim:={1}:

for i from 2 to nops(L) do
 if L[i]>si then
   alufim:={i}:
   si:=L[i]:
 elif L[i]=si then
      alufim:=alufim union {i}:
fi:
od:

print(`among the numbers in the list`, L, `the following places have the largest value`, alufim):
alufim:
end:

#BR1v(M): Verbose version of BR1(M). Inputs an m by n matrix of pairs of payoffs, where nops(M)=m, and player 1 has m options
#and n=nops(M[1]) and player n has n options, and M[i,j] is the pair [PayOffFor1, PayOffFor2]
#outputs the list of length n containing the sets of best reponses if Player 1 to each of the
#n strategies of Player 2.
#Try
#BR1v([[[2,8],[3,5]],[[8,2],[1,7]]]);
BR1v:=proc(M) local m,n,i,j,gu,lu,vu,M1:
m:=nops(M):
n:=nops(M[1]):
gu:=[]:
print(`Recall that the matrix of our game is`):
print(matrix(M)):

print(`Player Column has`, n, `possible strategies. Let's find Row's best response for each of them `):

print(`Since we are only interested right now in Row's pay-offs, let's list them here `):
M1:=[seq([seq(M[i][j][1],j=1..n)],i=1..m)]:
print(matrix(M1)):

for j from 1 to n do
print(`For Column's strategy Number`, j, `the respective pay-offs for each of Row's`, m, `strategies are `):
 lu:=[seq(M[i][j][1],i=1..m)]:
print(lu):
 vu:=MaxSetV(lu):
print(`Hence, Row's best response(s) is(are)`, vu):
 gu:=[op(gu),vu]:
od:
print(`To sum up, the table of Row's best responses for each of Column's`, n, `strategy are, respectively `):
print(gu):
gu:
end:


#BR2v(M): Verbose version of BR2(M). Inputs an m by n matrix of pairs of payoffs, where nops(M)=m, and player 1 has m options
#and n=nops(M[1]) and player n has n options, and M[i,j] is the pair [PayOffFor1, PayOffFor2]
#outputs the list of length n containing the sets of best reponses if Player 1 to each of the
#n strategies of Player 2.
#Try
#BR2v([[[2,8],[3,5]],[[8,2],[1,7]]]);
BR2v:=proc(M) local m,n,i,j,gu,lu,vu,M1:
m:=nops(M):
n:=nops(M[1]):
gu:=[]:
print(`Recall that the matrix of our game is`):
print(matrix(M)):

print(`Player Row has`, m, `possible strategies. Let's find Column's best response for each of them `):

print(`Since we are only interested right now in Column's pay-offs, let's list them here `):
M1:=[seq([seq(M[i][j][2],j=1..n)],i=1..m)]:
print(matrix(M1)):

for i from 1 to m do
print(`For Row's strategy Number`, i, `the respective pay-offs for each of Column's`, n, `strategies are `):
 lu:=[seq(M[i][j][2],j=1..n)]:
print(lu):
 vu:=MaxSetV(lu):
print(`Hence, Column's best response(s) is(are)`, vu):
 gu:=[op(gu),vu]:
od:
print(`To sum up, the table of Column's best responses for each of Row's`, m, `strategy are, respectively `):
print(gu):
gu:
end:

#PNEv(M): Verbose form of PNE(M)
#inputs a matrix M of payoffs and outputs the set of pure Nash Equilibira
#Try:
PNEv:=proc(M) local gu1,gu2,m,n, i,j,mu:
m:=nops(M):
n:=nops(M[1]):

gu1:=BR1v(M):
gu2:=BR2v(M):

mu:={}:


for i from 1 to m do
     print(`We are considering Strategy Number `, i, `of player Row`):
     print(`we saw above that the set of Column's best reponses to it `, gu2[i]):

     for j from 1 to n do
           print(`We are considering Strategy Number `, j, `of player Column`):

         if   not member(j,gu2[i])  then
          print(`Since strategy`, j, `is NOT in the above set, we alredy know that`, [i,j], `is NOT a Pure Nash Equilibrium` ):
         else
         
          print(`Since strategy`, j, `is in the above set there is hope, but now we have to check the second condition`):
          print(`Recall that the set of best responses of Row to Column's strategy`, j, ` is `, gu1[j]):

          if not member(i,gu1[j])  then

        print(`Since strategy`, i, `is NOT in the above set, we alredy know that`, [i,j], `is NOT a Pure Nash Equilibrium` ):
         else
         print(`Yea!, Row's Strategy`, i, `belongs to the set of Best Reponses of Column's Strategy`, j):
         print(`and vice-versa `):
         print(`Yea!, Column's Strategy`, j, `belongs to the set of Best Reponses of Row's Strategy`, i):
         print(` Hence `, [i,j], `is a Pure Nash Equilibrium `):
          mu:=mu union {[i,j]}:
         fi:
     fi:

 od:
od:

print(`To sum up we found that the set of pure Nash Equilibria for the above game is`, mu):
mu:

end:


#BRpV(f,p,q): verbose form of BRp(f,p,q): 
#inputs a bilinear expression in p,q and outputs the best p for a given q
#Try:
#BRpV(PO([[[2,8],[3,5]],[[8,2],[1,7]]],p,q),p,q);
BRpV:=proc(f,p,q) local gu,q0:
print(`Our function of `, p,q , `is `, f):
print(`We have to find, for every`, q , `between  0 and 1 the value(s) of`, p, `that will make it largest `):
gu:=coeff(f,p,1):

print(`The coefficient of p is the following expression linear function of q`, gu):

if min(subs(q=0,gu),subs(q=1,gu))>0 then

 print(`Since this is always positive (for all q from 0 to 1), the best respose is p=1, no matter what q is `):
     RETURN([[1,1],[0,1]]):
fi:

if max(subs(q=0,gu),subs(q=1,gu))<0 then
 print(`Since this is always negative (for all q), the best respose is p=0, no matter what q is `):
     RETURN([[0,0],[0,1]]):
fi:

q0:=solve(gu,q):

print(`The expression`, gu, `has a zero when `, q=q0):

if subs(q=0,gu)<0 then
 print(`when q is striclty less than`, q0, `this expression is negative, hence the best response is p=0 `):
 print(`when q equals`, q0, `this expression is zero, hence every possible p is a best response `):
 print(`when q is strictly larger than `, q0, `this is positive, hence the best response is p=1 `):

RETURN([ [[0,0],[0,q0]],[[0,1],[q0,q0]],[[1,1],[q0,1]] ]):
else
 print(`when q is striclty less than`, q0, `this expression is positive, hence the best response is p=1 `):
 print(`when q equals`, q0, `this expression is zero, hence every possible p is a best response `):
 print(`when q is strictly larger than `, q0, `this expression is negative, hence the best response is p=0 `):



RETURN([ [[1,1],[0,q0]],[[0,1],[q0,q0]],[[0,0],[q0,1]] ]):
fi:

end:

#MNEv(M): verbose form of MNE(M)
#inputs a matrix M of payoffs and outputs a Mixed Nash Equilibira

MNEv:=proc(M) local p,q, gu,mu1,mu2,p0,q0,lu1,lu2:


if not ( nops(M)=2  and nops(M[1])=2 ) then
print(`must be a 2 by 2 matrix`):
RETURN(FAIL):
fi:

print(`Our matrix game is`):
print(matrix(M)):

print(`Assume that player Row plays the Top strategy with probability p (and hence the Bottom strategy with prob. 1-p`):
print(`also assume that player Column plays the  Left strategy with probability q (and hence the Right strategy with prob. 1-q`):

gu:=PO(M,p,q):

lu1:=p*q*M[1][1][1]+p*(1-q)*M[1][2][1]+(1-p)*q*M[2][1][1]+(1-p)*(1-q)*M[2][2][1]:

print(`The expected pay-off of Row is`, lu1, `which simplifies to`, gu[1]):

lu2:=p*q*M[1][1][2]+p*(1-q)*M[1][2][2]+(1-p)*q*M[2][1][2]+(1-p)*(1-q)*M[2][2][2]:
print(`The expected pay-off of Column is`, lu2, `that simplifies to`, gu[2]):



print(`Let's consider Row's best response to Column's paying Left with prob. `, q):

mu1:=BRpV(gu[1],p,q):
print(``):
print(`---------------------------------------------`):
print(``):

print(`Now Let's consider Column's best response to Row's paying Top with prob. `, p):
mu2:=BRpV(gu[2],q,p):

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


if nops(mu1)=3 and nops(mu2)=3 then

p0:=mu1[2][2][1]:
q0:=mu2[2][2][1]:
print(`Since when p equals `, p0, `every q is OK and `):
print(`when q equals `, q0, `every p is OK and `):
print(`The mixed Nash equilibrium is`):
print([p0,q0]):
print(``):
print(`in other words Row should play Strategy Top with probability`, p0):
print(`and hence Strategy Bottom with probability`, 1-p0):
print(`and  Column should play Strategy left with probability`, q0):
print(`and hence Strategy Right with probability`, 1-q0):
print(``):
print(`With this mixed strategy the expected pay-off of player Row is`, subs({p=p0,q=q0},lu1) ):
print(`and the expected pay-off of player Column is`, subs({p=p0,q=q0},lu2) ):

 RETURN([mu1[2][2][1],mu2[2][2][1]]):
else
 print(`There are no mixed Nash Equilibrium, the pure ones are`):
 print(PNE(M)):
 RETURN(FAIL):
fi:

end: