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

with(combinat):
print(`First Written: Jan.  2024: tested for Maple 2020 `):
print(`Version  : Jan. 2024:  `):
print():
print(`This is Mourad.txt, A Maple package`):
print(`accompanying Shalosh B. Ekhad and Doron Zeilberger's  article: `):
print(` Efficient Weighted Counting of Multi-Set Derangements `):

print():
print(`The most current version is available on WWW at:`):
print(` http://sites.math.rutgers.edu/~zeilberg/tokhniot/Mourad.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():
 



ezra1:=proc()
if args=NULL then

print(`The SUPPORTING procedures are`):
print(` DerSeq, enkx, Lna, MSD, MSP, PerSeq, rf, Seqk, Stat, StatPer, Wder1 `):

else
ezra(args):

fi:


end:

ezra:=proc()
if args=NULL then

print(` Mourad.txt : A Maple package for  computing weighted (and unweighted) multi-set derangents, following in part the footsteps of Mourad Ismail  `):

print(`The MAIN procedures are`):
print(` CHk, Operk, RecStory, SeqkF, Wder`):



elif nargs=1 and args[1]=CHk then
print(`CHk(k,a,K): inputs a positive integer k, and a large integer K, and ouptuts a chapter about the sequence`):
print(`f_k(n,a):=weight-enumerator of multiset derangemets of 1(repeated k times), ..., n (repeated k times) according to a^NumberOfCycles Try:`):
print(`CHk(3,a,100);`):

elif nargs=1 and args[1]=DerSeq then
print(`DerSeq(a,N): The first N terms of the weight-enumerators of (usual) derangements according to the weight a^NumberOfCycles. Try:`):
print(`DerSeq(a,20);`):

elif nargs=1 and args[1]=enkx then
print(`enkx(n,k,x): The elementary symmetric function e_k(x[1], ..., x[n]). Try:`):
print(`enk(5,3,x);`):

elif nargs=1 and args[1]=Lna then
print(`Lna(n,a,x): the generalized Laguerre polynomials L_n^(a)(x). Try:`):
print(`Lna(4,a,x);`):


elif nargs=1 and args[1]=MSD   then
print(`MSD(L): Given a list  finds the set of all permutations of a a multi-set with nops(L) sets of sizes L[1],L[2], ..., L[nops(L)] elements respectively`):
print(`such that there are no incest.`):
print(`Try:`):
print(`MSD([2,2,2]);`):

elif nargs=1 and args[1]=MSP   then
print(` Given a sorted list w outputs all multi-set permutations with that w in TWO-LINE notation. Try:`):
print(`MSP([1,1,2,2,3,3]);`):


elif nargs=1 and args[1]=Operk then
print(`Operk(k,a,n,N): The operator annihilating the weight-enumerators sequence  of multi-set derangements of k^n, with weight a^NumberOfCycles. Try:`):
print(`Operk(3,a,n,N);`):


elif nargs=1 and args[1]=PerSeq then
print(`PerSeq(a,N): The first N terms of the weight-enumerators of permutations according to the weight a^NumberOfCycles. Try:`):
print(`PerSeq(a,20);`):

elif nargs=1 and args[1]=RecStory then
print(`RecStory(K): the recurrences for the weight-enumerators of multi-set derangemnts of 1^k...n^k for k from 1 to K. Try:`):
print(`RecStory(2);`):

elif nargs=1 and args[1]=rf then
print(`rf(a,k): the raising factorial a*(a+1)*...*(a+k-1). Try: `):
print(`rf(a,4);`):

elif nargs=1 and args[1]=Seqk then
print(`Seqk(k,K,a): the first K terms, starting with n=1 of the Wder([k$n],a). Very slow. Try:`):
print(`Seqk(3,8,a);`):


elif nargs=1 and args[1]=SeqkF then
print(`SeqkF(k,K,a): the first K terms, starting with n=1 of the Wder([k$n],a), the fast way, using the operator. Try:`):
print(`SeqkF(3,20,a);`):

elif nargs=1 and args[1]=Stat then
print(`Stat(F,t,K): Given a generating function F of t, for a r.v. finds the statistical information. Try:`):
print(`Stat((1+t)^n,t,4);`):

elif nargs=1 and args[1]=StatPer then
print(`StatPer(n,K): Same as Stat(PerSeq(a,n)[n],a,K) but via the formula. Try:`):
print(`Statper(20,6);`):

elif nargs=1 and args[1]=Wder then
print(`Wder(N,a): The sum of the weights of multi-set derangements for the list N, Try:`):
print(`Wder([3,4,5],a);`):

elif nargs=1 and args[1]=Wder1 then
print(`Wder1(N,a): Same as Wder(N,a) but via the generating function. Try:`):
print(`Wder1([3,4,5],a);`):

elif nargs=1 and args[1]=WderD then
print(`WderD(N,a): The sum of the weights of multi-set derangements for the list N done directly (just for checking). Try:`):
print(`WderD([2,1,2],a);`):

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

fi:


end:


Yafe:=proc(ope,N) local i,ope1,coe1,L: 
if ope=0 then
 RETURN(1,0):
fi:
ope1:=expand(ope):
L:=degree(ope1,N):
coe1:=coeff(ope1,N,L):
ope1:=normal(ope1/coe1):
ope1:=normal(ope1):
ope1:=
convert(
[seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`):
factor(coe1),ope1:
end:

 
#SeqFromRec(ope,n,N,Ini,K): Given the first L-1
#terms of the sequence Ini=[f(1), ..., f(L-1)]
#satisfied by the recurrence ope(n,N)f(n)=0
#extends it to the first K values
SeqFromRec:=proc(ope,n,N,Ini,K)
local ope1,gu,L,n1,j1:
ope1:=Yafe(ope,N)[2]:
L:=degree(ope1,N):
if nops(Ini)<>L then
 ERROR(`Ini should be of length`, L):
fi:

ope1:=expand(subs(n=n-L,ope1)/N^L):

gu:=Ini:

for n1 from nops(Ini)+1 to K do
gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)),j1=1..L)]:
gu:=normal(gu):
od:

gu:

end:


 #####Start from EKHAD
 
 
  
pashetd:=proc(p,N)
local gu,p1,mekh,i:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu,mekh):
end:

goremd:=proc(N,ORDER)
local i, gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*N^i:
od:
RETURN(gu):
end:

duisd:= proc(A,ORDER,y,n,N)
local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1:
KAMA:=20:
gorem:=goremd(N,ORDER):
 
 
conj:=gorem:
yakhas:=0:
 
for i1 from 0 to degree(conj,N) do
 
yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
 
 
shad:=pashetd(gorem,N):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
 
RETURN(shad[1],S):
end:

AZd:= proc(A,y,n,N)
local ORDER,gu,KAMA:
KAMA:=20:
 
for ORDER from 0 to KAMA do
gu:=duisd(A,ORDER,y,n,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
0:
end:
 
#End from EKAHD

#Stat(F,t,K): Given a generating function F of t, for a r.v. finds the statistical information. Try:
#Stat((1_t)^10000,4);
Stat:=proc(F,t,K) local gu,mu,sd,k,f:
f:=normal(F/subs(t=1,F)):
mu:=subs(t=1,diff(f,t)):
gu:=[mu]:
if K=1 then
 RETURN(gu):
fi:
f:=f/t^mu:
f:=t*diff(t*diff(f,t),t):
sd:=sqrt(subs(t=1,f)):
gu:=[op(gu),sd]:
for k from 3 to K do
 f:=t*diff(f,t):
gu:=[op(gu),subs(t=1,f)/sd^k]:
od:

gu:
end:

#rf(a,k): the rising factoria a*(a+1)*...*(a+k-1).
rf:=proc(a,k) local i:
mul(a+i,i=0..k-1):
end:

#Lna(n,a,x): the generalized Laguerre polynomials L_n^(a)(x). Try:
#Lna(4,a,x);
Lna:=proc(n,a,x) local k:
add((-1)^k*normal(rf(1+a,n)/rf(1+a,k))*x^k/k!/(n-k)!,k=0..n):
end:


#L1na(n,a,x): same as L(n,a,x) by via the generating function. Try
#L1na(4,a,x);
L1na:=proc(n,a,x) local t:
expand(coeff(taylor(exp(-x*t/(1-t))/(1-t)^(a+1),t=0,n+1),t,n)):
end:


#Wder(N,a): The sum of the weights of multi-set derangements for the list N, Try:
#Wder([3,4,5],a);
Wder:=proc(N,a) local i,x:
if nops(N)=1 then
 0:
else
expand(simplify(int(mul((-1)^N[i]*N[i]!*Lna(N[i],a-1,x),i=1..nops(N))*exp(-x)*x^(a-1),x=0..infinity)/(a-1)!)):
fi:
end:


#enkx(n,k,x): The elementary symmetric function e_k(x[1], ..., x[n]). Try:
#enk(5,3,x);
enkx:=proc(n,k,x) local i,t:
expand(coeff(mul(1+t*x[i],i=1..n),t,k)):
end:

#Wder1(N,a): Same as Wder(N,a) but via the generating function. Try:
#Wder1([3,4,5],a);
Wder1:=proc(N,a) local x,f,i,m:
m:=nops(N):
f:=1/(1-add((i-1)*enkx(m,i,x),i=2..m))^a:
for i from 1 to m do
f:=normal(coeff(taylor(f,x[i],N[i]+1),x[i],N[i])):
od:
mul(N[i]!,i=1..nops(N))*f:
end:


#Lna1(n,a,x): same as Lna(n,a,x) but via the generating function, just for checking. Try:
#Lna1(5,a,x);
Lna1:=proc(n,a,x) local u:
expand(coeff(taylor((1-u)^(-a-1)*exp(-x*u/(1-u)),u=0,n+1),u,n)):
end:


#DerSeq(a,N): The first N terms of the weight-enumerators of (usual) derangements according to the weight a^NumberOfCycles. Try:
#DerSeq(a,20);
DerSeq:=proc(a,N) local gu,i,x:
gu:=(1-x)^(-a)*exp(-a*x):
gu:=taylor(gu,x=0,N+1):
expand([seq(i!*coeff(gu,x,i),i=1..N)]):
end:



#Operk(k,a,n,N): The operator annihilating the weight-enumerators sequence  of multi-set derangements of k^n, with weight a^NumberOfCycles. Try:
#Operk(3,a,n,N);
Operk:=proc(k,a,n,N) local F,x:
option remember:
F:=(-1)^(k*n)*k!^n*Lna(k,a-1,x)^n*exp(-x)*x^(a-1):
AZd(F,x,n,N)[1]:
end:

#Seqk(k,K,a): the first K terms, starting with n=1 of the Wder([k$n],a). Very slow. Try:
#Seqk(3,8,a);
Seqk:=proc(k,K,a) local n1:
[seq(Wder([k$n1],a),n1=1..K)]:
end:


#SeqkF(k,K,a): the first K terms, starting with n=1 of the Wder([k$n],a) the fast way, using the operator. Try:
#SeqkF(3,30,a);
SeqkF:=proc(k,K,a) local n,N,ope,INI:
ope:=Operk(k,a,n,N):
INI:=Seqk(k,degree(ope,N),a):
expand(SeqFromRec(ope,n,N,INI,K)):
end:






#PerSeq(a,N): The first N terms of the weight-enumerators of permutations according to the weight a^NumberOfCycles. Try:
#PerSeq(a,20);
PerSeq:=proc(a,N) local gu,i,x:
gu:=(1-x)^(-a):
gu:=taylor(gu,x=0,N+1):
expand([seq(i!*coeff(gu,x,i),i=1..N)]):
end:



#MSP(w): Given a sorted list w outputs all multi-set permutations with that w in TWO-LINE notation
#Try:
#MSP([1,1,2,2,3,3]);
MSP:=proc(w) local gu,gu1,i1:
gu:=convert(permute(w),set):
{seq([seq([w[i1],gu1[i1]],i1=1..nops(w))],gu1 in gu)}:
end:

#IsDer(w): given a multi-set permutation in two line notation decides whether it is a derangement. Try
#IsDer([[1,2],[2,1]]);
IsDer:=proc(w) local i:
for i from 1 to nops(w) do
 if w[i][1][1]=w[i][2][1] then
   RETURN(false):
fi:
od:
true:
end:




#MSD(L): Given a list  finds the set of all permutations of a a multi-set with nops(L) sets of sizes L[1],L[2], ..., L[nops(L)] elements respectively
#such that there are no incest.
#Try:
#MSD([2,2,2]);
MSD:=proc(L) local lu,i1,gu,i,j,gu1,lu1:
lu:=[seq(seq([i,j],j=1..L[i]),i=1..nops(L))]:
gu:=permute(lu):

lu:={seq([seq([lu[i1],gu1[i1]],i1=1..nops(lu))],gu1 in gu)}:

gu:={}:

for lu1 in lu do
 if IsDer(lu1) then
  gu:=gu union {lu1}:
 fi:
od:
gu:

end:


#KhapesUp(w,b): Given a multi-set permutation w, and an integer b, in two-line notation finds the smallest i such w[i][1]=b
KhapesUp:=proc(w,b) local i:
for i from 1 to nops(w) while w[i][1]<>b do od:
if i=nops(w)+1 then
  RETURN(FAIL):
else
 RETURN(i):
fi:
end:


#ExtCyc(w): Given  a multi-permutation w in two-line notation, extracts the first  cycle, followed by the rest. Try:
#ExtCyc([[1,2],[1,2],[2,1],[2,1]]);
ExtCyc:=proc(w) local a,C1,b,i,C,i1,w1:
option remember:

if w=[] then
 RETURN([]):
fi:

C1:=[1]:
a:=w[1][1]:
b:=w[1][2]:

while a<>b do
i:=KhapesUp(w,b):
C1:=[op(C1),i]:
b:=w[i][2]:
od:


C:=[seq(w[C1[i1]],i1=1..nops(C1))]:
C1:=sort(convert({seq(i1,i1=1..nops(w))} minus convert(C1,set),list)):
w1:=[seq(w[C1[i1]],i1=1..nops(C1))]:
C:=[seq(C[i1][1],i1=1..nops(C))]:
[C,w1]:
end:


#CycD(w): The cycle decomposition of the multi-set permutation w try:
#CycD([[1,2],[1,2],[2,3],[2,3],[3,1],[3,1]]);
CycD:=proc(w) local gu,w1,lu,i:

w1:=w:
gu:=[]:

while w1<>[]  do
lu:=ExtCyc(w1):
gu:=[op(gu),lu[1]]:
w1:=lu[2]:
od:
gu:
end:


#WderD(N,a): The sum of the weights of multi-set derangements for the list N done directly (just for checking). Try:
#WderD([2,1,2],a);
WderD:=proc(N,a) local gu,i,gu1:
gu:=MSD(N):

add(a^nops(CycD(gu1)),gu1 in gu):

end:

#CHk(k,a,K): inputs a positive integer k, and a large integer K, and ouptuts a chapter about the sequence
#f_k(n):=weight-enumerator of multiset derangemets of 1(repeated k times), ..., n (repeated k times) according to a^NumberOfCycles Try:
#CHk(3,a,100);

CHk:=proc(k,a,K) local gu,ope,F,n,Sn,i:

if K<30 then
 print(K, `must be at least 30`):
RETURN(FAIL):
fi:
ope:=Operk(k,a,n,Sn):
gu:=SeqkF(k,K,a):

print(`On the Weight-Enumerator, according to the Number of Cycles, of permutations of a union of n disjoint sets, each with`, k, `members such that every member must go to a member of a different set`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let `, F[k](n), `be the polynomial in a whose coefficient of a^r is the number of permutations satisfying the condition stated in the title.`):
print(``):
print(`F[k](n)  satifies the following homogeneous linear recurrence equation of order`, degree(ope,Sn)):
print(``):
print(add(coeff(ope,Sn,i)*F[k](n+i),i=0..degree(ope,Sn))=0):
print(``):
print(`and in Maple notation`):
print(``):
lprint(add(coeff(ope,Sn,i)*F[k](n+i),i=0..degree(ope,Sn))=0):
print(``):
print(`Here are the first 30 terms, starting at n=1`):
print(``):
lprint([op(1..30,gu)]):
print(``):
print(`Just for fun here is the`, K, `-th term `):
print(``):
lprint(F[k](K)=gu[K]):
print(``):
print(`The average, standard deviation and the scaled moments up to the 6th are`):
print(``):
lprint(evalf(Stat(gu[K],a,6))):
print(``):

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

end:



#RecStory(K): the recurrences for the weight-enumerators of multi-set derangemnts of 1^k...n^k for k from 1 to K. Try:
#RecStory(2);
RecStory:=proc(K) local k,F,a,n,t0,ope,Sn,i:
t0:=time():
print(`Linear reccurrences for the weight-enumerators of multi-set derangements of 1^k...n^k for k from 1 to`, K):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):

for k from 1 to K do
ope:=Operk(k,a,n,Sn):
print(`Theorem number`, k, `: Let `, F[k](n), `be the weight-enumerator, according to the weight a^NumberOfCycles`):
print(`of permutations of a set that is the union of n disjoint sets, each with`, k, `elements such that no member can go to one of its set-mates`):
print(``):
print(`then we have the following linear recurrence.`):
print(``):
lprint(add(coeff(ope,Sn,i)*F[k](n+i),i=0..degree(ope,Sn))=0):
print(``):
od:
print(``):
print(`----------------------`):
print(``):
print(`This ends this article that took`, time()-t0, `seconds to generate `):
print(``):
end:


##start StatPer
ez:=proc(): print(` fk(p,k), mk(p,k), mamk(p,k), smk(n,k), StatPer(n,K) `):end:
with(combinat):
fk:=proc(p,k) local r:
option remember:
if k=0 then
1:
else
expand((-1)^(k-1)/k*add((-1)^r*fk(p,r)*p[k-r],r=0..k-1)):
fi:
end:

mk:=proc(p,k) local i:
add(stirling2(k,i)*i!*fk(p,i),i=0..k):
end:

#mamk(p,k)
mamk:=proc(p,k) local mu,r:
mu:=p[1]:
expand(add(mk(p,k-r)*binomial(k,r)*(-mu)^r,r=0..k)):
end:

#smk(n,k): the scaled k-th moment about the mean of the r.v. number of cycles in a permutation. Try:
#smk(30,4);
smk:=proc(n,k) local  p,gu,j,i1:
gu:=mamk(p,k)/mamk(p,2)^(k/2):

for j from 1 to k do
gu:=subs(p[j]=add(1/i1^j,i1=1..n),gu):
od:
gu:
end:



#StatPer(n,K): Same as Stat(PerSeq(a,n)[n],a,K) but via the formula. Try:
#Statper(20,6);
StatPer:=proc(n,K) local k,i1:
[add(1/i1,i1=1..n),sqrt(add(1/i1,i1=1..n)-add(1/i1^2,i1=1..n)),seq(smk(n,k),k=3..K)]:
end:
##End StatPer