Help11:=proc(): print(`Ex():RCT(n), FindDiv(M,W,pi), FindStable(M,W,K) , inv, GaleSh(M,W) `):end:

with(combinat):

Ex:=proc(): [ [2,3,4,1],[1,4,2,3],[4,3,1,2],[2,4,3,1]], [[1,4,2,3],[2,3,4,1],[4,2,1,3],[4,2,3,2]]:end:


#RCT(n): makes random-chice tables: pairs [M,W] where M and W are lists of lists of size n and M[i][j]:=The j-th bet woman of Mr. i and W[j][i]:=The i-th best man for Ms. j
RCT:=proc(n) local i:
[seq(randperm(n),i=1..n)], [seq(randperm(n),i=1..n)]:

end:


#FindDiv(M,W,pi): Given Ranking tables for men and women M,W, and  a current matching, tries to find a man and a woman who prefer to divorce their current spouses
#it outputs the new matching; If none exits, i.e. the matching is stable, it  return FAIL
FindDiv:=proc(M,W,pi) local n,i1,j1,i2,j2:
n:=nops(pi):
for i1 from 1 to n do

j1:=pi[i1]:

for i2 from i1+1 to n do

j2:=pi[i2]:


if M[i1][j2]<M[i1][j1] and W[j2][i1]<W[j2][i2] then
 RETURN([op(1..i1-1,pi),j2,op(i1+1..i2-1,pi),j1,op(i2+1..n,pi)]):
fi:

od:
od:

FAIL:

end:

#FindStable(M,W,K): Finds a stable matching with at most K iterations, followed by the number of iterations, starting with a random matching
FindStable:=proc(M,W,K) local n,pi,pi1,pi2,co:
n:=nops(M):
pi:=randperm(n):
pi1:=FindDiv(M,W,pi):


for co from 0 to K while pi1<>FAIL do
 pi2:=FindDiv(M,W,pi1):
co:=co+1:
pi:=pi1:
pi1:=pi2:
od:

if FindDiv(M,W,pi)<>FAIL then
print(`Could not find anything in`,K, ` propsoals` ):
RETURN( FAIL):
else
print(`The following is a stable matching`):
print(``):
print(pi):
print(``):
print(`It took`, co, `proposals `):

RETURN(pi,co):
fi:

end:


#inv(pi): The inv
inv:=proc(pi) local n,T,i:
n:=nops(pi):
for i from 1 to n do
T[pi[i]]:=i:
od:
[seq(T[i],i=1..n)]:
end:

#GaleSh(M,W): The Gale-Shapely algorithm where men are proposing
GaleSh:=proc(M,W) local n,i,j, Mrank,Wrank,CurH,CurW,SinM,SinW,i1,co:

n:=nops(M):

#co is a counter for the number of proposals
co:=0:

#Mrank[i]: Ranking by Mr. i of currently available women
for i from 1 to n do
 Mrank[i]:=inv(M[i]):
od:

#Wrank[j]: Ranking by Ms. j of currently available men
for j from 1 to n do
 Wrank[j]:=inv(W[j]):
od:


#CurH[i]: Mr. i's current wife (0 if he is single)
#CurW[j]: Ms. j's current husband (0 if she is single)
for i from 1 to n do
CurW[i]:=0:
CurH[i]:=0:
od:

#SimM: currently set of single men
SinM:={seq(i,i=1..n)}:

#SimW: currently set of single women
SinW:={seq(i,i=1..n)}:

while SinM<>{} do
#i is any currently single man
i:=SinM[1]:

#i proposes to his first currently available choice Mrank[i][1]:
j:=Mrank[i][1]:
i1:=CurH[j]:
co:=co+1:
print(` I am Mr.`, i, `Ms. `, j, `would you marry me?`):

if i1=0 then
 print(`Sure, since I am currently single`):
 SinM:=SinM minus {i}:
 SinW:=SinW minus {j}:
 CurW[i]:=j:
 CurH[j]:=i:

elif W[j][i1]<W[j][i] then
 print(`How dare you!, I am currently married to Mr.`, i1, `and I like him better than you!`):
  Mrank[i]:=[op(2..nops(Mrank[i]),Mrank[i])]:
else
 print(`Sure, while I am currently married to Mr. `, i1, `I prefer you, so I will dump him and marry you`):
 SinM:=SinM union {i1} minus {i}:
  Mrank[i1]:=[op(2..nops(Mrank[i1]),Mrank[i1])]:
 CurW[i]:=j:
 CurH[j]:=i:
fi:
od:

print(`The Stable matching is`):

for i from 1 to n do
print(` Mr.`, i, ` is married to  Ms.`, CurW[i]):
od:

print(`It took`, co, `proposals `):

[seq(CurW[i],i=1..n)],co:

end: