######################################################################
## GMIP.txt Save this file as  GMIP.txt  to use it,                  #
# stay in the                                                        #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read  `GMIP.txt` <Enter>                          #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
##  DoronZeil at gmail dot com                                       # 
######################################################################


print(`First Written: Jan. 16, 2025: tested for Maple 2020 `):
print(`Version  : Jan. 16, 2025 `):
print():
print(`This is GMIP.txt, A Maple package`):
print(`accompanying Shalosh B. Ekhad and Doron Zeilberger's  article: `):
print(` Experimenting with the Amazing Garsia-Milne Involution Principle`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/mamarim/mamarimhtml/gmip.html `):
print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/GMIP.txt .`):
print(`Please report all bugs to: DoronZeil at gmail dot com .`):
print():
print(`For general help, and a list of the MAIN functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name);" `):
print(`For a list of the supporting functions type: ezra1();`):
print():
print(`-------------------`):


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




ezraS:=proc()
if args=NULL then

print(`The Story procedures are:  GMstory, LimMomsPaper, Paper `):
print(``):

else
ezra(args):

fi:


end:



ezraSi:=proc()
if args=NULL then

print(`The Simulation procedures are:  Alpha, RandL, Simu, SimuS `):
print(``):

else
ezra(args):

fi:


end:

ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are: MaxC, MomsN, MomsS, MomsSam, Pars, Pars1, PLgfS, PLgfXs, plotDist, k, SMomsNam, TotC`):
print(``):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` GMIP.txt: A Maple package for  studying and experimentin with the Garsia-Milne Involution Principle`):

print(`The MAIN procedures are:  FC, FCnn, FCnn1, FCs, GM, GM1, GMdet, LimMoms, PLgf, PLgfX, Pnk, SMomsSam,  SMomsSsam`):


print(``):

elif nargs=1 and args[1]=Alpha then
print(`Alpha(f,x,K): Given a probab. generating function finds the average, variance, and scaled moments about the mean up to the K-th. Try:`):
print(`Alpha((1+x)^1000,x,8);`):

elif nargs=1 and args[1]=FC then
print(`FC(f,c,x): The probability generating function, according to complexity of all scenarios with f faithful man and women, and c cheaing men and women`):
print(`Try:`):
print(`FC(5,5,x);`):


elif nargs=1 and args[1]=FCnn then
print(`FCnn(n,x): same as FC(n,n,x) but using the recurrence produced by the Almkvist-Zeilberger algorithm. Try:`):
print(`FCnn(50,x);`):


elif nargs=1 and args[1]=FCnn1 then
print(`FCnn1(n,x): same as FnnC(n,x) but using little memory. Try:`):
print(`FCnn1(50,x);`):

elif nargs=1 and args[1]=FCs then
print(`FCs(f,c,x): The probability generating function, according to complexity of a SINGLE faithful man all scenarios with f faithful man and women, and c cheating men and women`):
print(`Try:`):
print(`FCs(5,5,x);`):

elif nargs=1 and args[1]=GM then
print(`GM(n,L): Given n and a "mistress list" L for the first k:=nops(L) cheating men outputs the list of length n-nops(L) for the matching of the faithful men. Try:`):
print(`GM(4,[3,4]);`):

elif nargs=1 and args[1]=GM1 then
print(`GM1(n,L,x): imputs a positive integer n, and a list of length k:=nops(L), L, where it is assumed that Mr. i has an affair with Mrs. L[i], and Mr. 1, ..., k`):
print(`all have affairs. Finds what faithful woman to match to Mr. x if x is faithful, i.e. x>k. Try, also outputs the  path, the last member is the faitfhul woman x is matched to. Try:`):
print(`GM1(4,[2,4],3);`):


elif nargs=1 and args[1]=GMdet then
print(`GMdet(n,L): Given n and a "mistress list" L for the first k:=nops(L) cheating men outputs the list of lists, of length n-nops(L) for the matching of the path leading from a faithful man. `):
print(`to a faithful woman. Try:`):
print(`GMdet(4,[3,4]);`):

elif nargs=1 and args[1]=GMstory then
print(`GMstory(n,L): Verbose form of GMdet(n,L). Try:`):
print(`GMstory(4,[3,4]);`):

elif nargs=1 and args[1]=LenVec then
print(`LenVec(n,L): the vector of length n-nops(L) with the path-length for the faithful men. Try:`):
print(`LenVec(10,[2,1,5,2,9]);`):

elif nargs=1 and args[1]=LimMoms then
print(`LimMoms(k,K): The list of length K whose first entry is the limit of the expectation, the second of the variance, and the next  moments about the mean up to the K-th`):
print(`for the r.v. length of matching path for a random faithful man where there are n faithful men (and women) and k*n cheating men (and women). Try:`):
print(`LimMoms(k,10);`):

elif nargs=1 and args[1]=LimMomsPaper then
print(`LimMomsPaper(k,K): a verbose form of LimNoms(k,K). Try:`):
print(`LimMomsPaper(k,10);`):

elif nargs=1 and args[1]=MaxC then
print(`MaxC(n,L): the maximum complexity of the complement of L. Try:`):
print(`MaxC(10,[2,1,5,2,9]);`):

elif nargs=1 and args[1]=MomsN then
print(`MomsN(f,c,K): The first K moments for the comlexity of the Garsia-Milne involution principle with f faithful men and women and c cheating men and women`):
print(`for NUMERIC f and c. Try:`):
print(`MomsN(5,5,4);`):

elif nargs=1 and args[1]=MomsS then
print(`MomsS(f,c,K): The first K moments for the comlexity of the Garsia-Milne involution principle with f faithful men and women and c cheating men and women`):
print(`for Symbolic f and c. Try:`):
print(`MomsS(f,c,4);`):

elif nargs=1 and args[1]=MomsSam then
print(`MomsSam(f,c,K): The first K moments ABOUT THE MEAN for the TOTAL comlexity of the Garsia-Milne involution principle with f faithful men and women and c cheating men and women`):
print(`for Symbolic f and c. Try:`):
print(`MomsSam(f,c,4);`):


elif nargs=1 and args[1]=MomsSsam then
print(`MomsSsam(f,c,K): The first K moments ABOUT THE MEAN for the INDVIDUAL comlexity of one (random)  faithful man of the Garsia-Milne involution principle with f faithful men and women and c cheating men and women`):
print(`for Symbolic f and c. Try:`):
print(`MomsSsam(f,c,4);`):

elif nargs=1 and args[1]=Paper then
print(`Paper(K): A paper about the average, variance, and K first scaled moments of the r.v. "complexity" of the Garsia-Milne Invoultion Principle`):
print(`for random scenarios with f faithful men (and women) and c faithful men (and women). Try:`):
print(`Paper(4);`):

elif nargs=1 and args[1]=Pars then
print(`Pars(k): the partitions with smallest number more than 1 that show up in PLgf(n,k) after removing the 1s together with the smallest n where they show up, Try:`):
print(`Pars(5);`):

elif nargs=1 and args[1]=Pars1 then
print(`Pars1(n,k): all the partitions that show up in PLgf(n,k)`):

elif nargs=1 and args[1]=PLgf then
print(`PLgf(n,k,t): the polynomial in t[1],..., t[n-k] whose coefficient of t[1]^a[1]*...*t[n-k]^a[n-k] is the number of mapping from the cheating men (assumed to be 1..k) into`):
prnt(`all women called 1,...,n, such that the Mr. k+1 has a path of length a[1], ... Mr. n has a path of length a[n-k] towards the faithful woman he is mapped to by the`):
print(`Garsia-Milne involution principle. Try:`):
print(`PLgf(6,4,t);`):


elif nargs=1 and args[1]=PLgfX then
print(`PLgfX(n,k,x): the polynomial in x whose coefficient of x^r is the number of mapping from the cheating men (assumed to be 1..k) into`):
prnt(`all women called 1,...,n, such that the total lengths of the path from faithful men to their match is r, via`):
print(`Garsia-Milne involution principle. Try:`):
print(`PLgfX(6,4,x);`):

elif nargs=1 and args[1]=PLgfXs then
print(`PLgfXs(n,k,x): Same as PLgfX(n,k,x) but much slower, done directly, just for checking. Try`):
print(`PLgfXs(6,4,x);`):


elif nargs=1 and args[1]=PLgfS then
print(`PLgfS(n,k,t): Like PLgf(n,k,t) but done directly. For checking purposes only. Try:`):
print(`PLgfS(6,4,t);`):

elif nargs=1 and args[1]=plotDist then
print(`plotDist(f,x,K1,K2): Given a prob. gen. function f(x) that has a `):
print(`Taylor series`):
print(`for a discrete r.v. plots it`):

elif nargs=1 and args[1]=Pnk then
print(`Pnk(n,k): all lists of distict elements from {1,...,n} with k elements. Try:`):
print(`Pnk(4,2);`):

elif nargs=1 and args[1]=RandL then
print(`RandL(n,k): a random list of length k with distinct elements from {1,...,n}. In other words a random injection from the k cheating men {1,...,k} into the set of n women {1,...,n}. Try:`):
print(`RandL(10,5);`):


elif nargs=1 and args[1]=Simu then
print(`Simu(f,c,K,N): doing the Garsia-Milne involution principle with f faithful men (and f faithful women) and c cheating men (and women) K times`):
print(`It returns the empirical average, variance and up to the N-th moment of the r.v. TOTAL comlexity. Try:`):
print(`Simu(100,100,300,4);`):


elif nargs=1 and args[1]=SimuS then
print(`SimuS(f,c,K,N): doing the Garsia-Milne involution principle with f faithful men (and f faithful women) and c cheating men (and women) K times`):
print(`It returns the empirical average, variance and up to the N-th moment of the r.v.  comlexity for ONE faithful man. Try:`):
print(`SimuS(100,100,300,4);`):

elif nargs=1 and args[1]=SMomsNam then
print(`SMomsNam(f,c,K): For NUMERIC f and c, The expectation, variance and SCALED moments about the mean for the r.v. "complexity of a random arragement with f faithful men and women and c cheating men and women, up to the K-th moment.Try:`):
print(`SMomsNam(f,c,6)`):

elif nargs=1 and args[1]=SMomsSam then
print(`SMomsSam(f,c,K): FOR SYMBOLIC f and c, The expectation, variance and SCALED moments about the mean for the r.v. "complexity of a random arragement with f faithful men and women and c cheating men and women, up to the K-th moment.Try:`):
print(`SMomsSam(f,c,6)`):


elif nargs=1 and args[1]=SMomsSsam then
print(`SMomsSsam(f,c,K): FOR SYMBOLIC f and c, The expectation, variance and SCALED moments about the mean for the r.v. The INDIVIDUAL complexity (length of path) for a random faithful man`):
print(` with a random arragement with f faithful men and women and c cheating men and women, up to the K-th moment.Try:`):
print(`SMomsSsam(f,c,6)`):

elif nargs=1 and args[1]=TotC then
print(`TotC(n,L): the total complexity of the complement of L. Try:`):
print(`TotC(10,[2,1,5,2,9]);`):


else
 print(`There is no such thing as`, args):

fi:


end:




with(combinat):


#plotDist(f,x,K): Given a prob. gen. function f(x) that has a 
#Taylor series
#for a discrete r.v.
#plots its normalized version (X-mu)/sig between mu-K1*sig and
#mu+K2*sig
#For example, try:
#plotDist((1+x)^40,x,5,5);
plotDist:=proc(f,x,K1,K2) local mu,f1,lu,gu,sig,i,j1:
f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
gu:=f1/x^mu:
sig:=sqrt(subs(x=1,x*diff(x*diff(gu,x),x))):

lu:=taylor(f1,x=0,trunc(mu)+K2*trunc(sig)+10):

lu:=[seq([i,coeff(lu,x,i)],i=max(0,trunc(mu-K1*sig))..trunc(mu+K2*sig))]:

lu:=evalf([seq([(lu[j1][1]-mu)/sig,lu[j1][2]*sig],j1=1..nops(lu))]):

plot(lu):
#plot(lu,scaling=constrained):

end:

#Pnk(n,k): all lists of distict elements from {1,...,n} with k elements. Try:
#Pnk(4,2);
Pnk:=proc(n,k) local gu,i,gu1:
option remember:
gu:=choose([seq(i,i=1..n)],k):
{seq(op(permute(gu1)),gu1 in gu)}:
end:

#GM1(n,L,x): imputs a positive integer n, and a list of length k:=nops(L), L, where it is assumed that Mr. i has an affair with Mrs. L[i], and Mr. 1, ..., k
#all have affairs. Finds what faithful woman to match to Mr. x if x is faithful, i.e. x>k. Try, also outputs the  path, the last member is the faitfhul woman x is matched to. Try:
#GM1(4,[2,4],3);
GM1:=proc(n,L,x) local k,FW,P,L1,i,y,CW:
k:=nops(L):
if not (x>k and x<=n) then
 print(x, `is not faithful`):
 RETURN(FAIL):
fi:

CW:={op(L)}:
FW:={seq(i,i=1..n)} minus CW:


if member(x, FW) then
 RETURN([x]):
fi:

for i from 1 to k do
 L1[L[i]]:=i:
od:


P:=[x]:
y:=x:
while member(y,CW) do
y:=L1[y]:
P:=[op(P),y]:
od:
P:
end:

#GM(n,L): Given n and a "mistress list" L for the first k:=nops(L) cheating men outputs the list of length n-nops(L) for the matching of the faithful men. Try:
#GM(4,[3,4]);
GM:=proc(n,L) local k,x:
k:=nops(L):

[seq(GM1(n,L,x)[-1],x=k+1..n)]:
end:


#GMdet(n,L): Given n and a "mistress list" L for the first k:=nops(L) cheating men outputs the list of lists of length n-nops(L) for the matching of the faithful men, each with the alternate path. Try:
#GMdet(4,[3,4]);
GMdet:=proc(n,L) local k,x:
k:=nops(L):

[seq(GM1(n,L,x),x=k+1..n)]:
end:



#RandL(n,k): a random list of length k with distinct elements from {1,...,n}. In other words a random injection from the k cheating men {1,...,k} into the set of n women {1,...,n}. Try:
#RandL(10,5);
RandL:=proc(n,k) local lu,pi,i:
lu:=convert(randcomb(n,k),list):
pi:=randperm(k):
[seq(lu[pi[i]],i=1..nops(pi))]:
end:



#LenVec(n,L): the vector of length n-nops(L) with the path-length for the faithful men. Try:
#LenVec(10,[2,1,5,2,9]);
LenVec:=proc(n,L) local i:
[seq(nops(GM1(n,L,i)),i=nops(L)+1..n)]:
end:

#TotC(n,L): the total complexity of the complement of L. Try:
#TotC(10,[2,1,5,2,9]);
TotC:=proc(n,L) local i:
add(nops(GM1(n,L,i)),i=nops(L)+1..n):
end:


#MaxC(n,L): the maximum complexity of the complement of L. Try:
#MaxC(10,[2,1,5,2,9]);
MaxC:=proc(n,L) local i:
max(seq(nops(GM1(n,L,i)),i=nops(L)+1..n)):
end:

#PLgfS(n,k,t): the polynomial in t[1],..., t[n-k] whose coefficient of t[1]^a[1]*...*t[n-k]^a[n-k] is the number of mapping from the cheating men (assumed to be 1..k) into
#all women called 1,...,n, such that the Mr. k+1 has a path of length a[1], ... Mr. n has a path of length a[n-k] towards the faithful woman he is mapped to by the
#Garsia-Milne involution principle. Try:
#PLgf(6,4,t);
PLgfS:=proc(n,k,t) local gu,gu1,i1:
gu:=Pnk(n,k):

add(mul(t[i1]^nops(GM1(n,gu1,k+i1)),i1=1..n-k),gu1 in gu):
end:

#Pars1(n,k): all the partitions that show up in PLgf(n,k)
Pars1:=proc(n,k) local t,gu,i,i1,lu,lu1:
gu:=PLgf(n,k,t):
lu:={seq([seq(degree(op(i,gu),t[i1]),i1=1..n-k)],i=1..nops(gu))}:
{seq(sort(lu1),lu1 in lu)}:
end:


Chop1:=proc(L) local i:
for i from 1 to nops(L) while L[i]=1 do od:
[op(i..nops(L),L)]:
end:

#Pars(k): the partitions with smallest number more than 1 that show up in PLgf(n,k) after removing the 1s together with the smallest n where there all whos up
Pars:=proc(k) local n,lu,lu1:
for n from k+1 while nops(Pars1(n+1,k))>nops(Pars1(n,k)) do od:
lu:=Pars1(n,k):
[n,{seq(Chop1(lu1),lu1 in lu)}]:
end:



#PLgfXs(n,k,x): the polynomial in x whose coefficient of x^i is the number of mapping from the cheating men (assumed to be 1..k) into
#all women called 1,...,n, such that the total path-length is i
#Garsia-Milne involution principle. Try:
#PLgfXs(6,4,X);
PLgfXs:=proc(n,k,x) local gu,t,i1:

gu:=PLgfS(n,k,t):
subs({seq(t[i1]=x,i1=1..n-k)},gu):
end:


#PLgf(n,k,t): Same as PLgfS(n,k,t) but the clever way! Try
#PLgf(6,4,t);
PLgf:=proc(n,k,t) local gu,z,i:
gu:=mul(t[i]/(1-t[i]*z),i=1..n-k)/(1-z):
expand(k!*coeff(taylor(gu,z=0,k+1),z,k)):
end:



#PLgfX(n,k,x): Same as PLgfX(n,k,x) but the clever way! Try
#PLgfX(6,4,t);
PLgfX:=proc(n,k,x) local gu,z:
gu:=(x/(1-x*z))^(n-k)/(1-z):
expand(k!*coeff(taylor(gu,z=0,k+1),z,k)):
end:


#FC(f,c,x): The probability generating function, according to complexity of all scenarios with f faithful man and women, and c cheating men and women
#Try:
#FC(5,5,x);
FC:=proc(f,c,x) local z,gu:
gu:=c!*f!/(c+f)!*(x/(1-x*z))^f/(1-z):
coeff(taylor(gu,z=0,c+1),z,c):
end:

#MomsN(f,c,K): The first K moments for the comlexity of the Garsia-Milne involution principle with f faithful men and women and c cheating men and women
#for NUMERIC f and c. Try:
#MomsN(5,5,4);
MomsN:=proc(f,c,K) local x,F,i,gu:
F:=FC(f,c,x):
gu:=[]:

for i from 1 to K do
 F:=x*diff(F,x):
 gu:=[op(gu),subs(x=1,F)]:
od:

gu:

end:


#MomsS(f,c,K): Explicit expressions for the first K moments for the comlexity of the Garsia-Milne involution principle with f faithful men and women and c cheating men and women
#for SYMBOLIC f and c. Try:
#MomsS(5,5,4);
MomsS:=proc(f,c,K) local x,F,i,gu,z,w,lu,lu1,j:
F:=c!*f!/(c+f)!*(x/(1-x*z))^f/(1-z):
gu:=[]:

for i from 1 to K do
 F:=x*diff(F,x):
 gu:=[op(gu),subs(x=1,F)]:
od:

gu:=normal(simplify(gu,symbolic)):

gu:=[seq(gu[i]/(-1+z)^(-f-1-i),i=1..K)]:
gu:=simplify(gu,symbolic):

gu:=expand(subs(z=1+w,gu)):
lu:=[]:

for i from 1 to K do
 #lu1:=add(coeff(gu[i],w,j)*(-1)^c*binomial(-f-i-1+j,c),j=0..degree(gu[i],w)):
   lu1:=add(coeff(gu[i],w,j)*(-1)^(f+i+1-j)*binomial(f+i+c-j,c),j=0..degree(gu[i],w)):
 lu1:=convert(lu1,factorial)/(-1)^(2*f):
 lu1:=simplify(lu1,symbolic):
 lu:=[op(lu),lu1]:
od:

lu:=simplify(lu,factorial):
 lu:=simplify(lu,symbolic):
end:



#MomsSam(f,c,K): The expectation, variance and moments about the mean for the r.v. "complexity of a random arragement with f faithful men and women and c cheating men and women, up to the K-th moment.Try:
#MomsSam(f,c,K)
MomsSam:=proc(f,c,K) local gu,mu,s,r:
gu:=MomsS(f,c,K):
mu:=gu[1]:
normal([mu,seq(add((-mu)^s*binomial(r,s)*gu[r-s],s=0..r-1)+(-mu)^r,r=2..K)]):
end:



#SMomsSam(f,c,K): The expectation, variance and SCALED moments about the mean for the r.v. "complexity of a random arragement with f faithful men and women and c cheating men and women, up to the K-th moment.Try:
#SMomsSam(f,c,K)
SMomsSam:=proc(f,c,K) local gu,i:
gu:=MomsSam(f,c,K):
[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..K)]:
end:


#SMomsNam(f,c,K): The expectation, variance and SCALED moments about the mean for the r.v. "complexity of a random arragement with f faithful men and women and c cheating men and women, up to the K-th moment.Try:
#SMomsNam(f,c,K)
SMomsNam:=proc(f,c,K) local gu,i:
gu:=MomsNam(f,c,K):
[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..K)]:
end:



#MomsNam(f,c,K): The expectation, variance and moments about the mean for the r.v. "complexity of a random arragement with f faithful men and women and c cheating men and women, up to the K-th moment.Try:
#For numeric f and c
#MomsNam(f,c,K)
MomsNam:=proc(f,c,K) local gu,mu,s,r:
gu:=MomsN(f,c,K):
mu:=gu[1]:
normal([mu,seq(add((-mu)^s*binomial(r,s)*gu[r-s],s=0..r-1)+(-mu)^r,r=2..K)]):
end:




#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:


#Alpha(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#Alpha(((1+x)/2)^100,x,4);
Alpha:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:


#Simu(f,c,K): doing the Garsia-Milne involution principle with f faithful men (and f faithful women) and c cheating men (and women) K times
#It returns the empirical average, variance and up to the N-th moment of the r.v. TOTAL complexity. Try:
#Simu(100,100,300,4);
Simu:=proc(f,c,K,N) local t,gu,i:
gu:=add(t^TotC(f+c,RandL(f+c,c)),i=1..K)/K:
evalf(Alpha(gu,t,N)):
end:


#SimuS(f,c,K): doing the Garsia-Milne involution principle with f faithful men (and f faithful women) and c cheating men (and women) K times
#It returns the empirical average, variance and up to the N-th moment of the r.v. Complexity for ONE faithful man. Try:
#SimuS(100,100,300,4);
SimuS:=proc(f,c,K,N) local t,gu,i:
gu:=add(t^nops(GM1(f+c,RandL(f+c,c),c+1)),i=1..K)/K:
evalf(Alpha(gu,t,N)):
end:

#Paper(K): A paper about the average, variance, and K first scaled moments of the r.v. "complexity" of the Garsia-Milne Invoultion Principle
#for random scenarios with f faithful men (and women) and c faithful men (and women). Try:
#Paper(4);
Paper:=proc(K) local i,gu,f,c,n,lu,x:
gu:=MomsSam(f,c,K):
print(``):
print(`Explicit formulas for the first`, K, `moments for the Complexity of a Garsia-Milne Involution mapping`):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`Suppose that there are f+c married couples, f of the men are faitfhul (and so are f of the women)  while c of the men (and c of the women) are cheating`):
print(``):
print(` The c cheating men pick  uniformly c cheating women, and then we use the Involution Principle to have a map the f faithful men to the f faithful women`):
print(``):
print(`Let X(f,c) be the sum of the lengths of the alternating paths`):

print(``):
print(`The expectation of X is`):
print(``):
print(gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):
print(``):

print(``):
print(`The variance of X is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[2]):
print(``):

for i from 3 to K do

print(``):
print(`The `, i, `-th  moment about the mean is`):
print(``):
print(gu[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[i]):
print(``):
od:

print(`Now let's see what if there are as many cheating men as faitfhul men, i.e. f=c=n`):
print(``):
gu:=subs({f=n,c=n},gu):
print(``):

print(``):
print(`The expectation of X is`):
print(``):
print(gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):
print(``):

print(``):
print(`The variance of X is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[2]):
print(``):

print(`As n goes to infinity it is`):
print(``):
lprint(limit(gu[2],n=infinity)):
print(``):

for i from 3 to K do

print(``):
print(`The `, i, `-th  moment about the mean is`):
print(``):
print(gu[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[i]):
print(``):

print(`As n goes to infinity it is`):
print(``):
lprint(limit(gu[i],n=infinity)):
print(``):

od:

print(`To sum up the expectation, variance and the  moments about the mean up to the`, K, `are `):
print(``):
lprint(gu):
print(``):
print(`the limits of the variance etc. are`):
print(``):
lu:=[seq(limit(gu[i],n=infinity),i=2..K)]:
lprint(lu):
print(``):
print(`This sequence is similar (with the odd entries with a minus sign) of OEIS sequence A53841 https://oeis.org/A052841, whose exponential generating function is`):
print(``):
print(exp(2*x)/(2*exp(x)-1)):
print(``):
print(``):

if  [seq(i!*coeff(taylor(exp(2*x)/(2*exp(x)-1),x=0,K+1),x,i),i=2..K)]=lu then
print(`that according to Svante Janson (see comment there) are the central limits of Ge(-1/2) `):
fi:

print(`-----------------------------------`):
print(``):
print(`This ends this paper that took`, time(), `second to produce`):
end:


#FCnn(n,x): a fast way to compute FC(n,n,x) using the recurrence outputted by the Almkvist-Zeilberger algorithm. Try:
#[seq(K(i,x),i=1..20)];
FCnn:=proc(n,x) 
option remember:
if n=1 then
 1/2*x^2+1/2*x:
elif n=2 then
1/2*x^4+1/3*x^3+1/6*x^2:
else
expand(normal((1/2*x*(8*n^2*x^3-12*n^2*x^2-16*n*x^3+2*n^2*x+26*n*x^2+6*x^3+n^2-7*n*x-12*x^2-2*
n+6*x)/(-1+x)/(2*n-1)/(2*n*x-n-3*x+2))*FCnn(n-1,x)+
(1/2*x^3*(n-1)*(2*n*x-n-x+1)/(-1+x)/(2*n
-1)/(2*n*x-n-3*x+2))*FCnn(n-2,x))):
fi:
end:

#FCnn1(n,x): a fast way to compute FC(n,n,x) using the recurrence and little memory
FCnn1:=proc(m,x) local A,B,C,n:

 A:=1/2*x^2+1/2*x :
B:=1/2*x^4+1/3*x^3+1/6*x^2:

if m=1 then
  RETURN(A):
elif m=2 then
  RETURN(B):
else
for n from 3 to m do
C:=expand(normal((1/2*x*(8*n^2*x^3-12*n^2*x^2-16*n*x^3+2*n^2*x+26*n*x^2+6*x^3+n^2-7*n*x-12*x^2-2*
n+6*x)/(-1+x)/(2*n-1)/(2*n*x-n-3*x+2))*B+
(1/2*x^3*(n-1)*(2*n*x-n-x+1)/(-1+x)/(2*n
-1)/(2*n*x-n-3*x+2))*A)):
A:=B:
B:=C:
od:
RETURN(C):
fi:
end:




####start for the complexity of a single faithful man

#FCs(f,c,x): The probability generating function, according to complexity of a SINGLE faithful man all scenarios with f faithful man and women, and c cheating men and women
#Try:
#FCs(5,5,x);
FCs:=proc(f,c,x) local z,gu:
gu:=c!*f!/(c+f)!*x/(1-x*z)/(1-z)^f:
coeff(taylor(gu,z=0,c+1),z,c):
end:




#MomsSs(f,c,K): Explicit expressions for the first K moments for the comlexity of a SINGLE faithful man in the Garsia-Milne involution principle with f faithful men and women and c cheating men and women
#for SYMBOLIC f and c. Try:
#MomsSs(5,5,4);
MomsSs:=proc(f,c,K) local x,F,i,gu,z,w,lu,lu1,j:
F:=c!*f!/(c+f)!*(x/(1-x*z))/(1-z)^f:
gu:=[]:

for i from 1 to K do
 F:=x*diff(F,x):
 gu:=[op(gu),subs(x=1,F)]:
od:

gu:=normal(simplify(gu,symbolic)):


gu:=[seq(gu[i]*(1-z)^(1+f+i),i=1..K)]:
gu:=simplify(gu,symbolic):


gu:=expand(subs(z=1-w,gu)):
lu:=[]:

for i from 1 to K do
   lu1:=add(coeff(gu[i],w,j)*binomial(f+i-j+c,c),j=0..degree(gu[i],w)):
 lu1:=convert(lu1,factorial)   :
 lu1:=simplify(lu1,symbolic):
 lu:=[op(lu),lu1]:
od:

lu:=simplify(lu,factorial):
 lu:=simplify(lu,symbolic):
end:


#MomsSsam(f,c,K): The expectation, variance and moments about the mean for the r.v. "complexity of a random arragement with f faithful men and women and c cheating men and women, up to the K-th moment.Try:
#MomsSsam(f,c,K)
MomsSsam:=proc(f,c,K) local gu,mu,s,r:
gu:=MomsSs(f,c,K):
mu:=gu[1]:
normal([mu,seq(add((-mu)^s*binomial(r,s)*gu[r-s],s=0..r-1)+(-mu)^r,r=2..K)]):
end:



#SMomsSsam(f,c,K): The expectation, variance and SCALED moments about the mean for the r.v. "complexity of a random arragement with f faithful men and women and c cheating men and women, up to the K-th moment.Try:
#SMomsSsam(f,c,K)
SMomsSsam:=proc(f,c,K) local gu,i:
gu:=MomsSsam(f,c,K):
[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..K)]:
end:




#GMstory(n,L): Verbose form of GMdet(n,L). Try
#GMstory(4,[3,4]);
GMstory:=proc(n,L) local k,gu,i,j,gu1:
k:=nops(L):
gu:=GMdet(n,L):

print(`Once upon a time there was a village with`, n,` married coupled. Let's call them Mr. i and Mrs. i, where Mr. i and Mrs. i are married to each other, where i goes from 1 to`, n):
print(``):
print(`It so happened that Mr. i for i=1 to i=`, k, `are cheaters and they each have one mistress (that may be their own wife)`):
print(``):
for i from 1 to k do
print(``):
print(`The Mistress of Mr.`, i, ` is Mrs.`, L[i]):
print(`and hence the Lover of Mrs.`, L[i], ` is Mr. `, i):
od:

print(`The `,k, `cheating men and their mistresses refuse to go to Church on Sunday, they rather go to the pub`):
print(``):
print(`There must be a way to match the`, n-k, `faitfhul men, Mr.`, {seq(i,i=k+1..n)}, `with the`, n-k, `faitfhul women Mrs. `):
print(``):
print(convert({seq(i,i=1..n)} minus {op(L)},list)):
print(``):

for i from 1 to nops(gu) do
print(``):
print(`----------------------------------`):
print(``):

gu1:=gu[i]:
print(`Mr. `, k+i, `first asks his wife`):
print(`My dear wife, will you go with me to Church?`):
if nops(gu1)=1 then
print(`She answers, of course dear! I am a faitfhul woman`):
print(`hence Mr.`, k+i, `and his wife Mrs. `, k+i, `happily go to Church together.`):
else
print(`Sorry, hubby, but I am going to the pub with my lover, Mr.`, gu1[2], `perhaps his wife Mrs.`, gu1[2], `is willling?`):
for j from 2 to nops(gu1)-1 do
print(`So Mr. `, k+i, `asks Mrs.`, gu1[j], `and she replies `):
print(`Sorry Mr.`, k+i, `I can't make it, I am going to the pub with my lover Mr. `, gu1[j+1], `why won't you ask his wife, Mrs.`, gu1[j+1]):
#print(`So he asks Mrs.`, gu1[j+1]):
od:
print(`So he asks Mrs.`, gu1[-1], `and she happily accepts, since she has no lover and rather go to Church`):
print(`So Mr,`, k+i, `goes to Church with Mrs.`, gu1[-1] ):
fi:
od:

print(`To sum up `):

print(`Regarding the sinners`):

for i from 1 to k do

print(`Mr.`, i, `goes to the pub with Mrs.`, L[i]):
od:

print(`Regarding faithful folks:`):

for i from 1 to n-k do
print(`Mr.`, k+i, `and Mrs.`, gu[i][-1], `happily go to Church together`):
od:
end:



#LimMoms(k,K): The list of length K whose first entry is the limit of the expectation, the second of the variance, and the next the scaled moments up to the mean
#for the r.v. length of matching path for a random faithful man where there are n faithful men (and wome) and kn cheating men (and women). Try:
#LimMoms(k,10);
LimMoms:=proc(k,K) local n,gu,i:
gu:=MomsSsam(n,k*n,K):
[seq(limit(gu[i],n=infinity),i=1..nops(gu))]:
end:

#LimMomsPaper(k,K): a verbose form of LimNoms(k,K). Try:
#LimMomsPaper(k,10);
LimMomsPaper:=proc(k,K) local gu,F,F1,t0,i,mu,mu1,x:
t0:=time():
gu:=LimMoms(k,K):

print(`Let k be a fixed constant and consdider a random configuration with n faitfful men (and n faithful women) and kn cheating men (and kn cheating women)`):
print(`consider the r.v. length of the path of any (specific) faithful man to the faithful woman matched to him via the Garsia-Milne`):
print(`Involution principal`):
print(``):
print(`Let's look at the limit as n goes to infinity`):
print(``):
lprint(`The expecation tends to`, gu[1]):
print(``):
lprint(`The variance tends to`, gu[2]):
print(``):
for i from 3 to K do
lprint(` The `, i, `-th moment about the mean is`, gu[i]):
print(``):
od:


print(`To sum up here they are`):

lprint(gu):

print(``):
mu:=subs(k=1,gu):

print(``):
print(`when k=1, i.e. there are as many faithful men as cheating men`):
print(``):
print(`The limits are`):
print(``):
lprint(mu):
print(``):
F:=1/exp(x)/(2-exp(x)):
F1:=taylor(F,x=0,K+3):

mu1:=[seq(i!*coeff(F1,x,i),i=2..K)]:

if [op(2..K,mu)]=mu1 then
 print(`Starting at the second moment about the mean, at least up to the`, K, `moment, the i-th moment is i! times the Taylor coefficient of x^i`):
print(``):
lprint(F):
print(``):
fi:

print(``):
print(`-----------------------`):
print(``):
print(`This ends this paper that took`, time()-t0, `seconds to make `):
print(``):
end:


#PaperS(K): A paper about the average, variance, and K first scaled moments of the r.v. "complexity" of an INDIVIDUAL faithful man in the Garsia-Milne Invoultion Principle
#for random scenarios with f faithful men (and women) and c faithful men (and women). Try:
#PaperS(4);
PaperS:=proc(K) local i,gu,f,c,n,lu,x:
gu:=MomsSsam(f,c,K):
print(``):
print(`Explicit formulas for the first`, K, `moments for the Complexity of a SINGLE faithful man in the Garsia-Milne Involution mapping`):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`Suppose that there are f+c married couples, f of the men are faitfhul (and so are f of the women)  while c of the men (and c of the women) are cheating`):
print(``):
print(` The c cheating men pick  uniformly c cheating women, and then we use the Involution Principle to have a map the f faithful men to the f faithful women`):
print(``):
print(`Let X(f,c) be the length of the alternating path leading from any faitfhul man (say the first) to the faitfhful woman that he gots matched to according to the `):
print(`Garsia-Milne Involution Principle`):

print(``):
print(`The expectation of X is`):
print(``):
print(gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):
print(``):

print(``):
print(`The variance of X is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[2]):
print(``):

for i from 3 to K do

print(``):
print(`The `, i, `-th  moment about the mean is`):
print(``):
print(gu[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[i]):
print(``):
od:

print(`Now let's see what if there are as many cheating men as faitfhul men, i.e. f=c=n`):
print(``):
gu:=subs({f=n,c=n},gu):
print(``):

print(``):
print(`The expectation of X is`):
print(``):
print(gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):
print(``):

print(``):
print(`The variance of X is`):
print(``):
print(gu[2]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[2]):
print(``):

print(`As n goes to infinity it is`):
print(``):
lprint(limit(gu[2],n=infinity)):
print(``):

for i from 3 to K do

print(``):
print(`The `, i, `-th  moment about the mean is`):
print(``):
print(gu[i]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[i]):
print(``):

print(`As n goes to infinity it is`):
print(``):
lprint(limit(gu[i],n=infinity)):
print(``):

od:

print(`To sum up the expectation, variance and the  moments about the mean up to the`, K, `are `):
print(``):
lprint(gu):
print(``):
print(`the limits of the variance etc. are`):
print(``):
lu:=[seq(limit(gu[i],n=infinity),i=2..K)]:
lprint(lu):
print(``):
print(`This sequence is the same as OEIS sequence A53841 https://oeis.org/A052841, whose exponential generating function is`):
print(``):
print(1/(exp(x)*(2-exp(x)))):

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

if  [seq(i!*coeff(taylor(1/(exp(x)*(2-exp(x))),x=0,K+1),x,i),i=2..K)]=lu then
print(`that according to Svante Janson (see comment there) are the central limits of Ge(1/2) `):
fi:

print(`-----------------------------------`):
print(``):
print(`This ends this paper that took`, time(), `second to produce`):
end: