Help:=proc(): print(`GM(n,pi), GMInfo(n,pi),CS(n,m),GMcom(n,m)`):end:

#GM(n,pi): Inputs a pos. integer n and
#a list pi listing
#the mistress of Mr. 1,Mr. 2, ..., Mr. m
#(where m=nops(pi))
#outputs the list sigma that finds mates
#for the faithful men Mr. m+1, ..., Mr. n
GM:=proc(n::posint,pi::list)
local i,m,HusbandOf, WifeOf,MistressOf,LoverOf,CM,CW,
FM,FW,w,sig:

m:=nops(pi):

CM:={seq(i,i=1..m)}:
CW:=convert(pi,set):
FM:={seq(i,i=1..n)} minus CM:
FW:={seq(i,i=1..n)} minus CW:



for i from 1 to n do
 HusbandOf[i]:=i:
 WifeOf[i]:=i:
od:


for i from 1 to m do
MistressOf[i]:=pi[i]:
LoverOf[pi[i]]:=i:
od:

for i from m+1 to n do
MistressOf[i]:=FAIL:
od:

for i from 1 to nops(FW) do
LoverOf[FW[i]]:=FAIL:
od:

sig:=[]:
 
for i from m+1 to n do
 w:=WifeOf[i]:
  
     while w in CW do
        w:=WifeOf[LoverOf[w]]:
     od:
  sig:=[op(sig),w]:
od:
sig:
  
end:

#GMInfo(n,pi): Inputs a pos. integer n and
#a list pi listing
#the mistress of Mr. 1,Mr. 2, ..., Mr. m
#(where m=nops(pi))
#outputs the list sigma that finds mates
#for the faithful men Mr. m+1, ..., Mr. n
GMInfo:=proc(n::posint,pi::list)
local i,m,HusbandOf, WifeOf,MistressOf,LoverOf,CM,CW,
FM,FW,w,sig,PathL,j:

m:=nops(pi):

CM:={seq(i,i=1..m)}:
CW:=convert(pi,set):
FM:={seq(i,i=1..n)} minus CM:
FW:={seq(i,i=1..n)} minus CW:



for i from 1 to n do
 HusbandOf[i]:=i:
 WifeOf[i]:=i:
od:


for i from 1 to m do
MistressOf[i]:=pi[i]:
LoverOf[pi[i]]:=i:
od:

for i from m+1 to n do
MistressOf[i]:=FAIL:
od:

for i from 1 to nops(FW) do
LoverOf[FW[i]]:=FAIL:
od:

sig:=[]:
 
PathL:=[]:

for i from m+1 to n do
 w:=WifeOf[i]:
j:=0:
  
     while w in CW do
        w:=WifeOf[LoverOf[w]]:
        j:=j+1:
     od:
  sig:=[op(sig),w]:
  PathL:=[op(PathL),j]:
od:
sig, PathL:
  
end:


#CS(n,m): the list of all cheating situations (a set of lists of
#length m
CS:=proc(n,m)
local S,i,T:
S:=choose(n,m):
T:=[]:

for i from 1 to nops(S) do
   T:=[op(T),op(permute(S[i]))]:
od:
T:


end:

#The list of all path-length of all possible cheating scenarios
GMcom:=proc(n,m)
local S,i:
S:=CS(n,m):

[seq(convert(GMInfo(n,S[i])[2],`+`)/(n-m),i=1..nops(S))]:

end:

GMav:=proc(n,m) local S:
S:=GMcom(n,m):
convert(S,`+`)/nops(S):

end: