######################################################################
## SMCper: Save this file as SMCper. To use it, stay in the          #
## same directory, get into Maple (by typing: maple <Enter> )        #
## and then type:  read SMCper : <Enter>                             #
## Then follow the instructions given there                          #
##                                                                   #
## Written by Doron Zeilberger, Rutgers University ,                 #
#  zeilberg@math.rutgers.edu.                                        # 
######################################################################

#Created: Sept. 2003-Jan. 2004
#This version:  Jan. 10, 2004
#SMCper: A Maple package to implement Doron Zeilberger's
#Symbolic Moment Calculus applied to patterns in permutations.
#Please report bugs to zeilberg@math.rutgers.edu



print(`Created: Sept 2003- Jan. 2004.`):
print(`This version: Jan. 10, 2003`):
print(`This is SMCper, a Maple package to implement`):
print(`Doron Zeilberger's Symbolic Moment Calculus`):
print(`applied to patterns in permutations `):
print(`as it is described in his article`):
print(` Symbolic Moment Calculus: I. Foundations and `):
print(`  Permutation-Patterns`):
print(``):
print(`Written by Doron Zeilberger, zeilberg@math.rutgers.edu`):
lprint(``):
print(`Please report bugs to zeilberg@math.rutgers.edu`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);`):
 print(``):
 print(`For a list of the supporting procedures type ezra1();, for help with`):
 print(`a specific supporting procedure, type ezra1(procedure_name);`):

 print(``):
with(combinat):
ezra1:=proc()
if args=NULL then

print(`The supporting procedures are`):
print(` ACorrSlow, ACovari,ACovariSlow, Ain, Bin`):
print(`CheckMishkalT1, CheckMishkalT2,Corr`):
 print(` CTU, kthMomentV,  kthMoment12n0`):
 print(`MishkalP, MishkalT, MishkalT0, MishkalT1, `):
 print(` MishkalT2, Profile3, Tupe, `):
print(` SetMishkalT, `):
fi:

if nops([args])=1 and op(1,[args])=ACorrSlow then
print(`ACorrSlow(Pat1,Pat2,n): Inputs two patterns (permutations not`):
print(`necessarily of the same length), Pat1, Pat2, and a symbol n`):
print(`outputs, the asymptotic (statistical) correlation of the`):
print(`random variables #Pat1 and #Pat2 defined on S_n (the set of`):
print(`permutations on {1, ..., n})`):
print(`Warning(This is the slow version just for testing ACorr`):
print(`You should always use ACorr).`):
print(`For example, try: ACorrSlow([1,2,3],[1,3,2],n);`):
fi:

if nops([args])=1 and op(1,[args])=ACovariSlow then
print(`ACovariSlow(Pat1,Pat2,n): The asymptotic covariance of`):
print(`Pat1 and Pat2 obtained by first computing the Covariance`):
print(`and then taking the leading term `):
print(`It is very slow, and should only be used for testing ACovari`):
fi:

if nops([args])=1 and op(1,[args])=Ain then
print(`Ain(i,n): E[binomial(X,i)] for the r.v. #12 over S_n (symbolic n)`):
fi:

if nops([args])=1 and op(1,[args])=Bin then
print(`Bin(i,n): E[X^i] for the r.v. #12 over S_n (symbolic n)`):
fi:


if nops([args])=1 and op(1,[args])=CheckMishkalT1 then
print(`CheckMishkalT1(m,n): checks procedure MishkalT1 for`):
print(`all possible inputs of patterns of length m and n`):
print(`(i.e. all m!*n! possibilities)`):
fi:

if nops([args])=1 and op(1,[args])=CheckMishkalT2 then
print(`CheckMishkalT2(m,n): checks procedure MishkalT2 for`):
print(`all possible inputs of patterns of length m and n`):
print(`(i.e. all m!*n! possibilities)`):
fi:


if nops([args])=1 and op(1,[args])=Corr then
print(`Corr(Pat1,Pat2,n): Inputs two patterns (permutations not`):
print(`necessarily of the same length), Pat1, Pat2, and a symbol n`):
print(`outputs the (statistical) correlation of the `):
print(`random variables #Pat1 and #Pat2 defined on S_n (the set of`):
print(`permutations on {1, ..., n})`):
print(`[Warning: it is very slow, to get the asymptotic correlation`):
print(`use ACorr`):
print(`For example, try: Corr([1,2,3],[1,3,2],n);`):
fi:

if nops([args])=1 and op(1,[args])=CTU then
print(`CTU(TU,Resh): Given a list of sets TU, covering, say`):
print(`{1, ..., t} of length k, say, and a list of patterns`):
print(`Resh, also of length k, such that nops(TU[i])=nops(Resh[i])`):
print(`finds the set of all permutations on {1,..t} such that`):
print(`the reduction of the sub-perm occupying the places`):
print(`indicated by TU[i] equals Resh[i] for all i in {1...t}`):
print(`For example, try: CTU([[1,2],[2,3]],[[1,2],[2,1]]);`):
fi:

if nops([args])=1 and op(1,[args])=kthMomentV then
print(`kthMomentV(Pat1,k,n): The k^th moment of the r.v. #Pat1 defined`):
print(` on S_n  (i.e. the average pf (#Pat1)^k.`):
print(` Warning: if k is larger than 2, it is very slow. `):
print(`For example, try kthMomentV([1,2],3,n);`):
fi:

if nops([args])=1 and op(1,[args])=kthMoment12Vn0 then
print(`kthMoment12Vn0(k0,n0): The k0^th moment of the r.v. #12, on S_n0`):
print(`(n0 an integer) a faster way`):
print(`(using generating functions) for testing `):
print(`subs(n=n0,kthMomentV([1,2],k0,n));`):
print(`For example, try: kthMoment12Vn0(3,5) `):
fi:

if nops([args])=1 and op(1,[args])=kthMoment12n0 then
print(`kthMoment12n0(k0,n0): The k0^th moment of the r.v. #12, on S_n0`):
print(`ABOUT the MEAN`):
print(`(n0 an integer) a faster way`):
print(`(using generating functions) for testing `):
print(`subs(n=n0,kthMoment([1,2],k,n));`):
print(`For example, try: kthMoment12n0(2,5) `):
fi:


if nops([args])=1 and op(1,[args])=MishkalP then
print(`MishkalP(Resh1,n,Depth1): The contribution to`):
print(`MishkalP(Resh1,n) from the Depth1 leading terms`):
print(`MishkalP(Resh1,n,gadol-katan+1) (where katan is the max of `):
print(`nops(Resh1[i]), (over all i)) and gadol is their sum`):
print(`should be the same as MishkalP(Resh1,n)`):
fi:

if nops([args])=1 and op(1,[args])=MishkalT then
print(`MishkalT(Resh1,t): Given a list of patterns Resh1 (of size k, say)`):
print(`and an integer t, finds how many pairs [[S1, ..., S_k],Pi] are there`):
print(`where the union of S_1,..., S_k is {1, ..., t} and the restriction`):
print(`of Pi to S_i reduces to Resh1[i] for i=1..k`):
print(`For example, try: MishkalT([[1,2],[2,1]],4);`);
fi:

if nops([args])=1 and op(1,[args])=MishkalT0 then
print(`MishkalT0(Resh1): A quick way to compute `):
print(`the (coeff. of the) leading contribution to Mishkal(Resh1,n)`):
print(`i.e. Mishkal(Resh1,n,sums(nops(Resh1[i]),i=1..nops(Resh1)))`):
print(`For example, try MishkalT0([[1,2],[1,2],[1,2]]);`):
fi:

if nops([args])=1 and op(1,[args])=MishkalT1 then
print(`MishkalT1(Resh1): A quick way to compute `):
print(`the (coeff. of the) 2nd-to-leading contribution to Mishkal(Resh1,n)`):
print(`i.e. Mishkal(Resh1,n,sum(nops(Resh1[i]),i=1..nops(Resh1))-1)`):
print(`For example, try MishkalT1([[1,2],[1,2],[1,2]]);`):
fi:

if nops([args])=1 and op(1,[args])=MishkalT2 then
print(`MishkalT2(Resh1): A quick way to compute `):
print(`the (coeff. of the) third-to-leading contribution to `):
print(`Mishkal(Resh1,n)`):
print(`i.e. Mishkal(Resh1,n,sums(nops(Resh1[i]),i=1..nops(Resh1))-2)`):
print(`So far Resh1 can only have length 2`):
print(`For example, try MishkalT2([[1,2],[1,2]]);`):
fi:

if nops([args])=1 and op(1,[args])=Profile3 then
print(`Profile3(Triplet): Given a triplet of sets`):
print(`outputs `):
print(`[nops(S1 intersect S2),nops(S1 intersect S3),nops(S2 intersect S3)]`):
print(`For example, try: Profile3([[1,2],[1,3],[2,3]]);`):
fi:

if nops([args])=1 and op(1,[args])=SetMishkalT then
print(`SetMishkalT(Resh1,t): Given a list of patterns Resh1 (of size k, say)`):
print(`and an integer t, finds the set of  pairs [[S1, ..., S_k],Pi] `):
print(`where the union of S_1,..., S_k is {1, ..., t} and the restriction`):
print(`of Pi to S_i reduces to Resh1[i] for i=1..k`):
print(`For example, try: SetMishkalT([[1,2],[2,1]],4);`);
fi:

if nops([args])=1 and op(1,[args])=Tupe then
print(`Tupe(Resh,S): Given a list of integers, Resh, and a set S`):
print(`Returns the set of tuples [S1,S2, ..., S_k] (where k=nops(Resh))`):
print(`such that the cardinality of S_i is Resh[i], and the union is S`):
print(`For example, try: Tupe([2,2],{1,2,3});`):
fi:
end:

ezra:=proc()
if args=NULL then

print(`Contains the following procedures: ACorr, `) :
print(`ACovari, AllACorr, AllACorrDiffSize, AllACorrDiffSizeV, AllACorrV, `):
print(`AllAThirdMom, AllkthMoments, AllkthMomentsE, AThirdMom, `):
print(`AllAVariance,AllAVarianceV, Covari, `):
print(` CovariE,`):
print(` EListE, EListES, kthMoment, kthMoment12, kthMomentE,Kurtosis, `):
print(`Mishkal `):
fi:

if nops([args])=1 and op(1,[args])=ACorr then
print(`ACorr(Pat1,Pat2): Inputs two patterns (permutations not`):
print(`necessarily of the same length), Pat1, Pat2.`):
print(`Outputs the asymptotic (statistical) correlation of the`):
print(`random variables #Pat1 and #Pat2 defined on S_n (the set of`):
print(` permutations of size n) `):
print(`(it is always a number independent of n)`):
print(`If it is O(1/n) then it outputs 0`):
print(`For example, try: ACorr([1,2,3],[1,3,2]);`):
fi:

if nops([args])=1 and op(1,[args])=ACovari then
print(`ACovari(Pat1,Pat2,n): The asymptotic covariance of`):
print(`the r.v. #Pat1 and #Pat2  defined on S_n, where`):
print(`Pat1, Pat2 are patterns, and n is a symbol`):
print(`for example, try: ACovari([1,2],[2,1],n);`):
fi:

if nops([args])=1 and op(1,[args])=AllkthMoments then
print(`AllkthMoments(a,k,n): The list of the formulas for the kth Moments`):
print(`over S_n (n a symbol) for patterns of size a`):
print(`For example, try AllkthMoments(3,2,n);`):
fi:

if nops([args])=1 and op(1,[args])=AllkthMomentsE then
print(`AllkthMomentsE(a,k,n): The list of the formulas for the kth Moments`):
print(`With the empirical (yet rigorous!) method `):
print(`over S_n (n a symbol) for patterns of size a`):
print(`For example, try AllkthMomentsE(3,2,n);`):
fi:

if nops([args])=1 and op(1,[args])=AThirdMom then
print(`AThirdMom(Pat1,n): The asymptotics of E[(X-E[X])^3]`):
print(`where X is the r.v. #Pat1 on S_n, where`):
print(`Pat1 is a pattern, and n is a symbol.`):
print(`For example, try: AThirdMom([1,2],n);`):
fi:

if nops([args])=1 and op(1,[args])=AllACorr then
print(`AllACorr(k): The list of all Asymptotic Correlations exactly and`):
print(`in floating-point, for all pairs of patterns of`):
print(`length k, in increasing order. For example, try: AllACorr(3);`):
fi:

if nops([args])=1 and op(1,[args])=AllACorrDiffSize then
print(`AllACorrDiffSize(k1,k2): The list of all Asymptotic Correlations`):
print(` exactly and in floating-point, for all pairs of patterns, one of`):
print(`length k1, and the other of length k2, For example, `):
print(`try: AllACorrDiffSize(2,3);`):
fi:

if nops([args])=1 and op(1,[args])=AllACorrDiffSizeV then
print(`AllACorrDiffSizeV(k1,k2): Verbose version of AllACorrDiffSize(k1,k2)`):
print(` (q,v,)`):
print(`For example, try: AllACorrDiffSizeV(2,3);`):
fi:

if nops([args])=1 and op(1,[args])=AllACorrV then
print(`A verbose version of AllACorr(k) (q.v)`):
print(`For example try: AllACorrV(3);`):
fi:

if nops([args])=1 and op(1,[args])=Covari then
print(`Covari(Pat1,Pat2,n): Inputs permutations (patterns) Pat1 and Pat2`):
print(`(of any length) and a symbol n, outputs the (polynomial) `):
print(` experssion for the Covariance of #Pat1 and #Pat2 on S_n`):
print(`For example, type:  Covari([1,2],[1,3,2],n); `):
fi:

if nops([args])=1 and op(1,[args])=CovariE then
print(`CovariE(Pat1,Pat2,n): The Covariance of #Pat1 and #Pat2 on S_n`):
print(`By the empirical method`):
print(`For example, type:  CovariE([1,2],[1,3,2],n); `):
fi:

if nops([args])=1 and op(1,[args])=EListE then
print(`EListE(Pats,n0): Given a list of patterns, and an integer n0`):
print(`outputs the expectation `):
print(`of the random variable (on the set of permutations of length n0)`):
print(`Product of `): 
print(`(#Occurences Pats[i]-Ave(%)), for i=1..nops(Pats).`):
print(``):
print(``):
print(`For example, to find the variance of the r.v. number of occurences`):
print(`of the pattern 123 in S_5 type: EListE([[1,2,3],[1,2,3]],5)`):
print(``):
print(`Another example, to find the covariance of the r.v.s `):
print(` "number of occurences of 123" and "number of occurences of 132" `):
print(`in S_4`):
print(`type: EListE([[1,2,3],[1,3,2]],4)`):
fi:

if nops([args])=1 and op(1,[args])=EListES then
print(`EListES(Pats,n): Given a list of patterns, and a symbol n`):
print(`outputs the expression for the expectation `):
print(`of the random variable (on the set of permutations of length n)`):
print(`Product of `): 
print(`(#Occurences Pats[i]-Ave(%)), for i=1..nops(Pats).`):
print(``):
print(`This procedure uses a completely experimental, empirical`):
print(` method, yet the answer is`):
print(`rigorous. The only drawback is that the running time is super-exp`):
print(`(in the sum of the lengths of the patters)`):
print(``):
print(`For example, to find the variance of the r.v. number of occurences`):
print(`of the pattern 123, type: EListES([[1,2,3],[1,2,3]],n)`):
print(``):
print(`Another example, to find the covariance of the r.v.s `):
print(` "number of occurences of 123" and "number of occurences of 132" `):
print(`type: EListES([[1,2,3],[1,3,2]],n)`):
fi:

if nops([args])=1 and op(1,[args])=AllAThirdMom then
print(`AllAThirdMom(k): The list of patterns of length k`):
print(`according to their asymptotic third-moment (about the mean)`):
print(` of the r.v. #Pat,`):
print(`(defined on S^n) from smallest`):
print(`to biggest (divided by n^(3k-2))`):
print(`For example, try AllAThirdMom(3);`):
fi:

if nops([args])=1 and op(1,[args])=kthMoment then
print(`kthMoment(Pat1,k,n): The k^th moment ABOUT THE MEAN of`):
print(` the r.v. #Pat1 defined`):
print(` on S_n  (i.e. the average op  (#Pat1-E[#Pat1])^k.`):
print(` Warning: if k is larger than 2, it is very slow. `):
print(`For example, try kthMoment([1,2],3,n);`):
fi:

if nops([args])=1 and op(1,[args])=kthMoment12 then

print(`kthMoment12(k,n): E[(X-E[X])^k] for the r.v. #12 over S_n `):
print(`(symbolic n) `):
print(` [Like kthMoment([1,2],k,n), but much faster (using the g.f`):
print(`of the number of inversions)`):
print(`Try kthMoment12(4,n)`):

fi:

if nops([args])=1 and op(1,[args])=kthMomentE then
print(`kthMomentE(Pat1,k,n): Like kthMoment(Pat1,k,n) (q.v.) `):
print(`but using the Empirical (yet rigorous!) method`):
print(` Warning: if k*nops(Pat1) is larger than 9, it is very slow. `):
print(`For example, try kthMomentE([1,2],3,n);`):
fi:

if nops([args])=1 and op(1,[args])=Kurtosis then
print(`Kurtosis(Pat1,n): The  Kurtosis of the r.v. #Pat1 defined on S_n`):
print(`Warning: For n>2 it takes too much memory and time.`):
print(`For example, try: Kurtosis([1,2],n):`):
fi:

if nops([args])=1 and op(1,[args])=Mishkal then
print(`Mishkal(Resh1,n): Given a list of patterns Resh1`):
print(`and a symbol n, computes the Expection of the r.v.`):
print(` #Resh1[1](sigma) x #Resh1[2](sigma) x ... x#Resh1[k](sigma)`):
print(` over S_n `):
print(`For example, try: Mishkal([[1,2],[2,1]],n);`);
fi:

if nops([args])=1 and op(1,[args])=AllAVariance then
print(`AllAVariance(k): The list of patterns of length k`):
print(`according to their asymptotic variance of the r.v. #Pat,`):
print(`(defined on S^n) from smallest`):
print(`to biggest (divided by n^(2k-1))`):
print(`For example, try AllAVariance(3);`):
fi:

if nops([args])=1 and op(1,[args])=AllAVarianceV then
print(`AllAVarianceV(k): Verbose version of AllAVariance(k) (q.v.)`):
print(`For example, try AllAVarianceV(3);`):
fi:

end:

#empirical part
ez:=proc():print(`Kama1(pi,sig), EListE(Pats,n), PerLi(n) `):
print(`EListES(Pats,n), TotalWtProf110(Resh)`):
end:

#redu(perm): Given a permutation on a set of positive integers
#finds its reduced form

redu:=proc(perm)
local gu,i,ra,mu:
gu:=sort(perm):

for i from 1 to nops(gu) do
ra[gu[i]]:=i:
od:

mu:=[]:

for i from 1 to nops(perm) do
mu:=[op(mu),ra[perm[i]]]:
od:
mu:
end:


#IV(n,k) all increasing vectors of length
#k of the form 1<=i_1< ...<i_k <=n 

IV:=proc(n,k)
local i,gu,mu:

if k=0 then
RETURN({[]}):
fi:

if n<k then
RETURN({}):
fi:

gu:=IV(n-1,k):
mu:=IV(n-1,k-1):

for i from 1 to nops(mu) do
 gu:=gu union {[op(op(i,mu)),n]}:
od:

gu:
end:

#Kama1(pi,sig): Counts how many times the pattern sig appears in
#the perm pi
Kama1:=proc(pi,sig)
local n,k,mu,mu1,co,i1,i:
n:=nops(pi):
k:=nops(sig):

mu:=IV(n,k):

co:=0:

for i from 1 to nops(mu) do
mu1:=mu[i]:

if redu([seq(pi[mu1[i1]],i1=1..nops(mu1))])=sig then
 co:=co+1:
fi:

od:
co:
end:





PerLi:=proc(n)
local p,p1,d,A,i,Don,m,j:
for i from 1 to n+1 do
p[i]:=i:
p1[i]:=i:
d[i]:=-1:
od:
p[0]:=n+1:
A:={seq(i,i=2..n)}:

Don:=false:

while not Don do
 print([seq(p[i],i=1..n)]):

 if A<>{} then
     m:=max(op(A)):
     j:=p1[m]:
     p[j]:=p[j+d[m]]:
     p[j+d[m]]:=m:
     p1[m]:=p1[m]+d[m]:
     p1[p[j]]:=j:
       if m<p[j+2*d[m]] then
           d[m]:=-d[m]:
            A:=A minus {m}:
            A:=A union {seq(i,i=m+1..n)}:
        fi:
   else
        Don:=true:
  fi:
od:
end:


EListE:=proc(Pats,n)
local p,p1,d,A,i,Don,m,j,pi,gu:
gu:=0:
for i from 1 to n+1 do
p[i]:=i:
p1[i]:=i:
d[i]:=-1:
od:
p[0]:=n+1:
A:={seq(i,i=2..n)}:

Don:=false:

while not Don do
pi:=[seq(p[i],i=1..n)]:
gu:=gu+convert([seq(
Kama1(pi,Pats[j])-binomial(n,nops(Pats[j]))/nops(Pats[j])!,
j=1..nops(Pats))],`*`):
 if A<>{} then
     m:=max(op(A)):
     j:=p1[m]:
     p[j]:=p[j+d[m]]:
     p[j+d[m]]:=m:
     p1[m]:=p1[m]+d[m]:
     p1[p[j]]:=j:
       if m<p[j+2*d[m]] then
           d[m]:=-d[m]:
            A:=A minus {m}:
        fi:
            A:=A union {seq(i,i=m+1..n)}:
   else
        Don:=true:
  fi:
od:
gu/n!:
end:


GuessPol:=proc(lu,x)
local i1,deg,a,eq,var,P:
deg:=nops(lu)-1:

P:=0:
var:={}:
for  i1 from 0 to deg do
P:=P+a[i1]*x^i1:
var:=var union {a[i1]}:
od:
eq:={}:
for i1 from 1 to nops(lu) do
eq:=eq union {subs(x=i1-1,P)-lu[i1]}:
od:
var:=solve(eq,var):

subs(var,P):
end:



EListES:=proc(Pats,n)
local k,i:
option remember:
k:=convert([seq(nops(Pats[i]),i=1..nops(Pats))],`+`):
factor(GuessPol([seq(EListE(Pats,i),i=0..k)],n)):
end:


###end of empirical (yet rigorous) part





#Subsets1(S,k): all the k-subsets of S
Subsets1:=proc(S,k) local i,gu1,gu2,S1,N: option remember:

if k=0 then
RETURN({{}}):
fi:

if nops(S)=0 or nops(S)<k or k<0 then
RETURN({}):
fi:
N:=max(op(S)):
S1:=S minus {N}:
gu1:=Subsets1(S1,k):
gu2:=Subsets1(S1,k-1):
gu2:={seq(gu2[i] union {N},i=1..nops(gu2))}:
gu1 union gu2:
end:

#CTU(TU,Resh): Given a list of sets TU, covering, say
#{1, ..., t} of length k, say, and a list of patterns
#Resh, also of length k, such that nops(TU[i])=nops(Resh[i])
#finds the set of all permutations on {1,..t} such that
#the reduction of the sub-perm occupying the places
#indicated by TU[i] equals Resh[i] for all i in {1...t}
CTU:=proc(TU,Resh)
local k,i,t,S,TU1,L1,Resh1,gu,pla,akha1,akha2,cand,GU,j,mu,pi,sig,i1:

k:=nops(TU):

if nops(Resh)<>k then
 ERROR(`Bad input`):
fi:

if {seq(evalb(nops(TU[i])=nops(Resh[i])),i=1..k)}<>{true} then
 ERROR(`Bad input`):
fi:

S:={seq(op(TU[i]),i=1..nops(TU))}:
t:=nops(S):

if S<>{seq(i,i=1..t)} then
 ERROR(`Bad input`):
fi:

if k=1 then
RETURN({Resh[1]}):
fi:

TU1:=[op(1..k-1,TU)]:
akha1:=TU[k]:
Resh1:=[op(1..k-1,Resh)]:
akha2:=Resh[k]:
L1:={seq(op(TU1[i]),i=1..k-1)}:
pla:=sort(convert(S minus L1,list));

TU1:=Redu1(TU1):
gu:=CTU(TU1,Resh1):
mu:=permute([seq(i,i=1..t)],nops(pla)):

GU:={}:
for i from 1 to nops(gu) do
 sig:=gu[i]:
 for j from 1 to nops(mu) do
  pi:=mu[j]:

    cand:=Insert1(sig,pla,pi):
   if redu([seq(cand[akha1[i1]],i1=1..nops(akha1))])=akha2 then
        
         GU:=GU union {cand}:
   fi:
  od:
od:

#The three lines below maybe uncommented in order to test this procedure
#if {seq(CheckComp(TU,Resh,GU[i]),i=1..nops(GU))}<>{{true}} then
#ERROR(`Something is wrong!`):
#fi:
GU:
end:

Yofi:=proc(Resh)
local i:
[seq(sort(convert(Resh[i],list)),i=1..nops(Resh))]:
end:

#Tup(Resh,S): Given a list of integers Resh, let's call
#[k_1, ..., k_r], and a set S
#finds all tuples  [S_1, ..., S_r],
#where S_1, ..., S_r are subsets of S and their union
#equals S and nops(S_i)=k_i
Tup:=proc(Resh,S)
local i,gu,r,mu,lu,j,Resh1,gu1,k:
option remember:
r:=nops(Resh):
if r=1 then
if nops(S)=Resh[1] then
  RETURN({[S]}):
else
 RETURN({})
fi:
fi:

Resh1:=[op(1..r-1,Resh)]:
lu:=Subsets1(S,Resh[r]):

gu:={}:

for i from 1 to nops(lu) do
mu:=powerset(lu[i]):

for j from 1 to nops(mu) do
gu1:= Tup(Resh1,S minus mu[j]):
  gu:=gu union {seq([op(gu1[k]),lu[i]],k=1..nops(gu1))}:
od:
od:


gu:
end:


Tupe:=proc(Resh,S)
local gu,i,GU,lu:
GU:={}:
gu:=Tup(Resh,S):
for i from 1 to nops(gu) do
lu:=gu[i]:
lu:=Yofi(lu):
GU:=GU union {lu}:
od:
GU:
end:


Mishkal:=proc(Resh1,n)
local mu,gu,i,t,j,mish,mish1,katan,gadol:
option remember:
mu:=[seq(nops(Resh1[i]),i=1..nops(Resh1))]:
katan:=max(op(mu)):
gadol:=convert(mu,`+`):
mish:=0:
for t from katan to gadol do
gu:=Tupe(mu,{seq(i,i=1..t)}):
  mish1:=0:
 for j from 1 to nops(gu) do
   mish1:=mish1+nops(CTU(gu[j], Resh1)):
 od:
 mish:=expand(mish+expand(binomial(n,t))/t!*mish1):
od:
factor(expand(mish)):
end:



#redu(perm): Given a permutation on a set of positive integers
#finds its reduced form

redu:=proc(perm)
local gu,i,ra,mu:
gu:=sort(perm):

for i from 1 to nops(gu) do
ra[gu[i]]:=i:
od:

mu:=[]:

for i from 1 to nops(perm) do
mu:=[op(mu),ra[perm[i]]]:
od:
mu:

end:


#Redu1(TU): Given a list TU of sets covering a set S
#finds the isomorphic
Redu1:=proc(TU) local L1,S,k,i,j:
L1:=sort(convert({seq(op(TU[i]),i=1..nops(TU))},list)):
k:=nops(L1):

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

[seq([seq(S[TU[i][j]],j=1..nops(TU[i]))],i=1..nops(TU))]:
end:


#Insert1(sig,pla,pi): Given a permutation sig, a list of places
#pla, and a gen. permutation pi, constructs a new permutation
#such that the places indicated by pla are exactly occopied by
#pi, and the remiander reduces to sig.
#For example, try: Insert1([1,2,3],[4,5],[1,3]);
Insert1:=proc(sig,pla,pi)
local t,i,L1,L2,T,ku:
t:=nops(sig)+nops(pla):
if {op(pi)} minus {seq(i,i=1..t)}<>{} then
ERROR(`Bad input`):
fi:

L1:=sort(convert({seq(i,i=1..t)} minus {op(pla)},list)):
L2:=sort(convert({seq(i,i=1..t)} minus {op(pi)},list)):


T:=Expand1(L1,L2,sig):


for i from 1 to nops(pla) do
 T[pla[i]]:=pi[i]:
od:

ku:=[seq(T[i],i=1..t)]:

ku:
end:


#Expand1(L1,L2,pi): given an increasing list of integers L1, 
#and another increasing list of integers L2, of the same
#length (say k), and a permutation pi of length , outputs a table, T,
#such that redu(T[L1[1]], ..., T[L1[k]])=pi
#For example, try Expand1([3,5],[1,8],[2,1]);

Expand1:=proc(L1,L2,pi)
local T,i,k:
k:=nops(pi):

if not (k=nops(L1) and k=nops(L2)) then
 ERROR(`Bad input`):
fi:

for i from 1 to k do
 T[L1[i]]:=L2[pi[i]]:
od:

op(T):
end:

#Covari(Pat1,Pat2,n): The Covariance of #Pat1 and #Pat2 on S_n
Covari:=proc(Pat1,Pat2,n) local k1,k2:
k1:=nops(Pat1):k2:=nops(Pat2):
factor(expand(Mishkal([Pat1,Pat2],n)-expand(binomial(n,k1)/k1!)*
expand(binomial(n,k2)/k2!))):
end:

#CovariE(Pat1,Pat2,n): The Covariance of #Pat1 and #Pat2 on S_n
#by the emprical method
CovariE:=proc(Pat1,Pat2,n):
EListES([Pat1,Pat2],n):
end:

MishkalT:=proc(Resh1,t)
local mu,gu,i,j,mish1:
option remember:
mu:=[seq(nops(Resh1[i]),i=1..nops(Resh1))]:

gu:=Tupe(mu,{seq(i,i=1..t)}):
  mish1:=0:
 for j from 1 to nops(gu) do
   mish1:=mish1+nops(CTU(gu[j], Resh1)):
 od:
mish1:
end:


#MishkalP(Resh1,n,Depth1): The contribution to
#MishkalP(Resh1,n) from the Depth1 leading terms
#MishkalP(Resh1,n,gadol-katan+1) (where katan is the max of 
#nops(Resh1[i]), (over all i)) and gadol is their sum
#should be the same as MishkalP(Resh1,n)

MishkalP:=proc(Resh1,n,Depth1)
local mu,i,t,mish,katan,gadol:
option remember:
mu:=[seq(nops(Resh1[i]),i=1..nops(Resh1))]:
katan:=min(op(mu)):
gadol:=convert(mu,`+`):
mish:=0:
for t from gadol by -1 to max(gadol-Depth1+1,katan) do
 mish:=expand(mish+expand(binomial(n,t))/t!*MishkalT(Resh1,t)):
od:
factor(expand(mish)):
end:

#MishkalT0(Resh1): A quick way to compute 
#the (coeff. of the) leading contribution to Mishkal(Resh1,n)
#i.e. Mishkal(Resh1,n,sums(nops(Resh1[i]),i=1..nops(Resh1)))
#For example, try MishkalT0([[1,2],[1,2],[1,2]]);
MishkalT0:=proc(Resh1)
local t,mu,i:
mu:=[seq(nops(Resh1[i]),i=1..nops(Resh1))]:
t:=convert(mu,`+`):
(t!/convert([seq(mu[i]!,i=1..nops(mu))],`*`))^2:
end:
 
MTC:=proc()
local gu,i:
gu:=[args]:
for i from 1 to nops(gu) do
if gu[i]<0 then
 ERROR(`All arguments must be non-negative `):
fi:
od:
convert(gu,`+`)!/
convert([seq(gu[i]!,i=1..nops(gu))],`*`):
end:



#CheckMishkalT1(m,n): checks procedure MishkalT1 for
#all possible inputs of patterns of length m and n
#(i.e. all m!*n! possibilities)
CheckMishkalT1:=proc(m,n)
local gu1,gu2,i,j:
gu1:=permute(m):
gu2:=permute(n):
{seq(seq(evalb(MishkalT1([gu1[i],gu2[j]])=MishkalT([gu1[i],gu2[j]],m+n-1)),
i=1..nops(gu1)),j=1..nops(gu2))}:
end:

#ACovariSlow(Pm 
ACovariSlow:=proc(Pat1,Pat2,n) local gu,d:
gu:=expand(Covari(Pat1,Pat2,n)):d:=degree(gu,n):coeff(gu,n,d)*n^d:
end:

#ACovari(Pat1,Pat2,n): The asymptotic covariance of
#the r.v. #Pat1 and #Pat2  defined on S_n, where
#Pat1, Pat2 are patterns, and n is a symbol
#for example, try: ACovari([1,2],[2,1],n);
ACovari:=proc(Pat1,Pat2,n) local gu,a,b,d:
a:=nops(Pat1):b:=nops(Pat2):
gu:=expand(
MishkalT0([Pat1,Pat2])*expand(binomial(n,a+b))/(a+b)!+
MishkalT1([Pat1,Pat2])*expand(binomial(n,a+b-1))/(a+b-1)!
-expand(binomial(n,a))/a!*expand(binomial(n,b))/b!):
d:=degree(gu,n):coeff(gu,n,d)*n^d:
end:


#ACovariB(Pat1,Pat2,n): The Better asymptotic covariance of
#the r.v. #Pat1 and #Pat2  (with two leading terms) defined on S_n, where
#Pat1, Pat2 are patterns, and n is a symbol
#for example, try: ACovari([1,2],[2,1],n);
ACovariB:=proc(Pat1,Pat2,n) local gu,a,b,d:
a:=nops(Pat1):b:=nops(Pat2):
gu:=expand(
MishkalT0([Pat1,Pat2])*expand(binomial(n,a+b))/(a+b)!+
MishkalT1([Pat1,Pat2])*expand(binomial(n,a+b-1))/(a+b-1)!+
MishkalT2([Pat1,Pat2])*expand(binomial(n,a+b-2))/(a+b-2)!
-expand(binomial(n,a))/a!*expand(binomial(n,b))/b!):
d:=degree(gu,n):coeff(gu,n,d)*n^d+coeff(gu,n,d-1)*n^(d-1)
:
end:


#Corr(Pat1,Pat2,n): Inputs two patterns (permutations not
#necessarily of the same length), Pat1, Pat2, and a symbol n
#outputs,
#The (statistical) correlation of the
#random variables #Pat1 and #Pat2 defined on S_n (the set of
#permutations on {1, ..., n})
#For example, try: Corr([1,2,3],[1,3,2],n);
Corr:=proc(Pat1,Pat2,n):
Covari(Pat1,Pat2,n)/
sqrt(Covari(Pat1,Pat1,n)*Covari(Pat2,Pat2,n)):
end:

#ACorrSlow(Pat1,Pat2,n): Inputs two patterns (permutations not
#necessarily of the same length), Pat1, Pat2, and a symbol n
#outputs,
#The asymptotic (statistical) correlation of the
#random variables #Pat1 and #Pat2 defined on S_n (the set of
#permutations on {1, ..., n})
#Warning(This is the slow version just for testing ACorr
#You should always use ACorr).
#For example, try: ACorrSlow([1,2,3],[1,3,2],n);
ACorrSlow:=proc(Pat1,Pat2,n) local gu,lu1,lu2,d,Cgu,Clu1,Clu2,d1,d2:
gu:=expand(Covari(Pat1,Pat2,n)):

lu1:=expand(Covari(Pat1,Pat1,n)):
lu2:=expand(Covari(Pat2,Pat2,n)):
d:=degree(gu,n):
d1:=degree(lu1,n):d2:=degree(lu2,n):
Cgu:=coeff(gu,n,d):
Clu1:=coeff(lu1,n,d1):Clu2:=coeff(lu2,n,d2):
Cgu/sqrt(Clu1*Clu2)*n^(d-(d1+d2)/2):

end:

#ACorr(Pat1,Pat2): Inputs two patterns (permutations not
#necessarily of the same length), Pat1, Pat2, 
#and outputs,
#The asymptotic (statistical) correlation of the
#random variables #Pat1 and #Pat2 defined on S_n (the set of
#permutations on {1, ..., n}) (it is always a number independent of n)
#For example, try: ACorr([1,2,3],[1,3,2]);
ACorr:=proc(Pat1,Pat2) local n,gu,lu1,lu2,d,Cgu,Clu1,Clu2,d1,d2:
gu:=expand(ACovari(Pat1,Pat2,n)):

lu1:=expand(ACovari(Pat1,Pat1,n)):
lu2:=expand(ACovari(Pat2,Pat2,n)):
d:=degree(gu,n):
d1:=degree(lu1,n):d2:=degree(lu2,n):
Cgu:=coeff(gu,n,d):
Clu1:=coeff(lu1,n,d1):Clu2:=coeff(lu2,n,d2):
if d-(d1+d2)/2<0 then
RETURN(0):
elif d-(d1+d2)/2=0 then
RETURN(Cgu/sqrt(Clu1*Clu2)):
else
ERROR(`Something is wrong`):
fi:

end:


#AllAVariance(k): The list of patterns of length k
#according to their asymptotic variance of the r.v. #Pat,
#(defined on S^n) from smallest
#to biggest (divided by n^(2k-1))
#For example, try AllAVariance(3); 
AllAVariance:=proc(k)
local gu,T,i,pi,lu,S,n,mu:
option remember:
gu:=Reps1(k):
mu:={}:

for i from 1 to nops(gu) do
pi:=gu[i]:
lu:=ACovari(pi,pi,n)/n^(2*k-1):
T[pi]:=lu:
mu:=mu union {lu}:
od:

for i from 1 to nops(mu) do
S[mu[i]]:={}:
od:

for i from 1 to nops(gu) do
pi:=gu[i]:
S[T[pi]]:=S[T[pi]] union {pi}:
od:
mu:=convert(mu,list):
mu:=sort(mu):
[seq([S[mu[i]],mu[i],evalf(mu[i])],i=1..nops(mu))]:
end:


#AllAVarianceV(k): A Verbose version of
#AllAVariance(k) (q.v.)
AllAVarianceV:=proc(k)
local gu,i:
gu:=AllAVariance(k):

print(`There are`, nops(gu), `different asymptotic variances for`):
print(`the random variable #Occurrences of the pattern pi, of length`, k):
print(`defined over all permutations of length n`):
print(`The ratio between the biggest and the smallest is`):
print(evalf(gu[nops(gu)][2]/gu[1][2])):

for i from 1 to nops(gu) do
print(`The `, i, `-th smallest Asymptotic variance belongs to the `):
print(`permutations in the set (and their obvious images)`):
print(gu[i][1]):
print(`whose asymptotic variance is  n to the power`, (2*k-1), ` times ` ):
print(gu[i][2], `which is roughly`, gu[i][3]):
od:
end:

#AllACorr(k): The list of all Asymptotic Correlation exactly and
#in floating-point, for all pairs of patterns of
#length k. For example, try: AllACorr(3);
AllACorr:=proc(k)
local gu,i,mu,ku,Tva,S,V:

mu:=Reps2(k):

gu:=[]:

for i from 1 to nops(mu) do
  ku:=ACorr(op(mu[i])):
  gu:=[op(gu),[mu[i],ku,evalf(ku)]]:
od: 

Tva:={seq(gu[i][3],i=1..nops(gu))}:

for i from 1 to nops(Tva) do
S[Tva[i]]:={}:
od;

for i from 1 to nops(gu) do
 S[gu[i][3]]:=S[gu[i][3]] union {gu[i][1]}:
 V[gu[i][3]]:=gu[i][2]:
od:

Tva:=sort(convert(Tva,list)):

[seq([S[Tva[i]],V[Tva[i]],Tva[i]],i=1..nops(Tva))]:

end:

#AllACorrV(k): A verbose version of AllACorr(k) (q.v.)
AllACorrV:=proc(k) local gu,i:
gu:=AllACorr(k):
print(`There are `, nops(gu) , `different values that the asymptotic `):
print(`Correlation of pairs of patterns of length`, k, ` can take. `):
print(`Here there are in increasing order`):

for i from 1 to nops(gu) do
print(`The `, i, `-th  smallest asymptotic correlation is`, gu[i][2]):
print(`which is approximately`, gu[i][3]):

print(`It is realized by the following pairs (and their obvious images)`):
print(gu[i][1]):
od:
end:

#AllACorrDiffSize(k1,k2): The list of all Asymptotic Covariances exactly and
#in floating-point, for all pairs of patterns, one of
#length k1, and the other of length k2, For example, 
#try: AllACorrDiffSize(2,3);
AllACorrDiffSize:=proc(k1,k2)
local gu,i,j,mu,ku,Tva,S,V:

mu:=Reps2DiffSize(k1,k2):

gu:=[]:
for i from 1 to nops(mu) do
ku:=ACorr(op(mu[i])):
gu:=[op(gu),[mu[i],ku,evalf(ku)]]:
od:

Tva:={seq(gu[i][3],i=1..nops(gu))}:

for i from 1 to nops(Tva) do
S[Tva[i]]:={}:
od;

for i from 1 to nops(gu) do
 V[gu[i][3]]:=gu[i][2]:
 S[gu[i][3]]:=S[gu[i][3]] union {gu[i][1]}:
od:

Tva:=sort(convert(Tva,list)):

[seq([S[Tva[i]],V[Tva[i]],Tva[i]],i=1..nops(Tva))]:

end:


#AllACorrDiffSizeV(k1,k2): A verbose version of AllACorrDiffSize(k1,k2) (q.v.)
AllACorrDiffSizeV:=proc(k1,k2) local gu,i:
gu:=AllACorrDiffSize(k1,k2):
print(`There are `, nops(gu) , `different values that the asymptotic `):
print(`Correlation of pairs of patterns one of length`, k1):
print(`the other of length`, k2, ` can take `):
print(`Here there are in increasing order`):

for i from 1 to nops(gu) do
print(`The `, i, `-th  smallest Asymptotic Correlation is`, gu[i][2]):
print(`which is appx.`, gu[i][3]):
print(`It is realized by the following pairs (and their images)`):
print(gu[i][1]):
od:
end:


#MishkalT1(Resh1): A quick way to compute 
#the (coeff. of the) second-to-leading contribution to Mishkal(Resh1,n)
#i.e. Mishkal(Resh1,n,sums(nops(Resh1[i]),i=1..nops(Resh1))-1)
#For example, try MishkalT1([[1,2],[1,2],[1,2]]);
MishkalT1:=proc(Resh1) local pi,sig,a,b,r,s,schum,mu,i,j,i1,T,t,schum1,i2:

schum:=0:
for i from 1 to nops(Resh1) do
for j from i+1 to nops(Resh1) do
mu:=[seq(nops(Resh1[i1]),i1=1..i-1),
seq(nops(Resh1[i1]),i1=i+1..j-1),
seq(nops(Resh1[i1]),i1=j+1..nops(Resh1)),nops(Resh1[i])+nops(Resh1[j])-1
]:
t:=convert(mu,`+`):
T:=(t!/convert([seq(mu[i2]!,i2=1..nops(mu))],`*`))^2:

pi:=Resh1[i]:
sig:=Resh1[j]:
a:=nops(pi):
b:=nops(sig):
schum1:=0:
for r from 1 to a do
 for s from 1 to b do
      schum1:=schum1+binomial(r+s-2,r-1)*binomial(a+b-r-s,a-r)*
          binomial(pi[r]+sig[s]-2,pi[r]-1)*
          binomial(a+b-pi[r]-sig[s],a-pi[r]):

od:
od:         

schum:=schum+schum1*T:
od:
od:


schum: 

end:


#AThirdMom(Pat1,n): The asymptotics of E[(X-E[X])^3]
#the r.s. #Pat1 and #Pat2  defined on S_n, where
#Pat1, Pat2 are patterns, and n is a symbol
#for example, try: AThirdMom([1,2],n);
AThirdMom:=proc(Pat1,n) local gu,a,d,gu3,gu2,gu1:
a:=nops(Pat1):
gu3:=expand(
MishkalT0([Pat1,Pat1,Pat1])*expand(binomial(n,3*a))/(3*a)!+
MishkalT1([Pat1,Pat1,Pat1])*expand(binomial(n,3*a-1))/(3*a-1)!+
MishkalT2([Pat1,Pat1,Pat1])*expand(binomial(n,3*a-2))/(3*a-2)!):
gu1:=expand(binomial(n,a))/a!:
gu2:=ACovari2(Pat1,Pat1,n):
gu:=expand(expand(gu3-3*gu2*gu1-gu1^3)):
d:=3*a-2:coeff(gu,n,d)*n^d:
end:





SetMishkalT:=proc(Resh1,t)
local mu,gu,i,j,mish1,k,ku:
option remember:
mu:=[seq(nops(Resh1[i]),i=1..nops(Resh1))]:
mish1:={}:
gu:=Tupe(mu,{seq(i,i=1..t)}):
  mish1:={}:
 for j from 1 to nops(gu) do
    ku:=CTU(gu[j], Resh1):
  for k from 1 to nops(ku) do
   mish1:=mish1 union {[gu[j],ku[k]]}:
 od:
 od:
mish1:
end:

#CheckMishkalT2(m,n): checks procedure MishkalT2 for
#all possible inputs of patterns of length m and n
#(i.e. all m!*n! possibilities)
ChecMkishkalT2:=proc(m,n)
local gu1,gu2,i,j:
gu1:=permute(m):
gu2:=permute(n):

{seq(seq(
evalb(MishkalT2([gu1[i],gu2[j]])=
MishkalT([gu1[i],gu2[j]],m+n-2)),
i=1..nops(gu1)),j=1..nops(gu2))}:
end:


#ACovari2(Pat1,Pat2,n): The asymptotic covariance of
#the r.v. #Pat1 and #Pat2 (with one extra term) defined on S_n, where
#Pat1, Pat2 are patterns, and n is a symbol
#for example, try: ACovari2([1,2],[2,1],n);
ACovari2:=proc(Pat1,Pat2,n) local gu,a,b,d:
a:=nops(Pat1):b:=nops(Pat2):
gu:=expand(
MishkalT0([Pat1,Pat2])*expand(binomial(n,a+b))/(a+b)!+
MishkalT1([Pat1,Pat2])*expand(binomial(n,a+b-1))/(a+b-1)!+
MishkalT2b([Pat1,Pat2])*expand(binomial(n,a+b-2))/(a+b-2)!
-expand(binomial(n,a))/a!*expand(binomial(n,b))/b!):
d:=degree(gu,n):coeff(gu,n,d)*n^d+coeff(gu,n,d-1)*n^(d-1):
end:


#MishkalT2a(Resh1): A quick way to compute 
#the contribution to
#the (coeff. of the) third-to-leading contribution to Mishkal(Resh1,n)
#from triplets sharing one location
#So far Resh1 can only have length 3
#For example, try MishkalT2a([[1,2],[1,2],[1,2]]);
MishkalT2a:=proc(Resh1) local pi,sig,mu,a,b,c,r,s,t,schum:

if nops(Resh1)<>3 then
ERROR(`nops(Resh1) must be 3, the general case is not yet implemented`):
fi:

pi:=Resh1[1]:
sig:=Resh1[2]:
mu:=Resh1[3]:
a:=nops(pi):
b:=nops(sig):
c:=nops(mu):
schum:=0:

for r from 1 to a do
 for s from 1 to b do
  for t from 1 to c do
      schum:=schum+
          MTC(r-1,s-1,t-1)*
          MTC(a-r,b-s,c-t)*
          MTC(pi[r]-1,sig[s]-1,mu[t]-1)*
          MTC(a-pi[r],b-sig[s],c-mu[t]):
od:
od:
od:         

schum: 

end:



#MishkalT2b(Resh1): A quick way to compute 
#the (coeff. of the) third-to-leading contribution to Mishkal(Resh1,n)
#(coming from interaction between two patterns only)
#i.e. Mishkal(Resh1,n,sums(nops(Resh1[i]),i=1..nops(Resh1))-2)
#For example, try MishkalT2b([[1,2],[1,2],[1,2]]);
MishkalT2b:=proc(Resh1) local pi,sig,a,b,r1,r2,s1,s2,
schum,mu,i,j,i1,T,t,schum1,i2:

schum:=0:
for i from 1 to nops(Resh1) do
for j from i+1 to nops(Resh1) do
mu:=[seq(nops(Resh1[i1]),i1=1..i-1),
seq(nops(Resh1[i1]),i1=i+1..j-1),
seq(nops(Resh1[i1]),i1=j+1..nops(Resh1)),nops(Resh1[i])+nops(Resh1[j])-2
]:
t:=convert(mu,`+`):
T:=(t!/convert([seq(mu[i2]!,i2=1..nops(mu))],`*`))^2:

pi:=Resh1[i]:
sig:=Resh1[j]:
a:=nops(pi):
b:=nops(sig):
schum1:=0:



for r1 from 1 to a do
for r2 from r1+1 to a do
 for s1 from 1 to b do
  for s2 from s1+1 to b do
     if (pi[r1]<pi[r2] and sig[s1]<sig[s2]) then
      schum1:=schum1+binomial(r1+s1-2,r1-1)*
                   binomial(r2+s2-r1-s1-2,r2-r1-1)*
                   binomial(a+b-r2-s2,a-r2)*
          binomial(pi[r1]+sig[s1]-2,pi[r1]-1)*
          binomial(pi[r2]-pi[r1]+sig[s2]-sig[s1]-2,pi[r2]-pi[r1]-1)*
          binomial(a+b-pi[r2]-sig[s2],a-pi[r2]):
      fi:

     if (pi[r2]<pi[r1] and sig[s2]<sig[s1]) then
      schum1:=schum1+binomial(r1+s1-2,r1-1)*
                   binomial(r2+s2-r1-s1-2,r2-r1-1)*
                   binomial(a+b-r2-s2,a-r2)*
          binomial(pi[r2]+sig[s2]-2,pi[r2]-1)*
          binomial(pi[r1]-pi[r2]+sig[s1]-sig[s2]-2,pi[r1]-pi[r2]-1)*
          binomial(a+b-pi[r1]-sig[s1],a-pi[r1]):
      fi:

od:
od:         
od:
od:


schum:=schum+schum1*T:
od:
od:


schum: 

end:



#kthMomentV(Pat1,k,n): The k^th moment of the r.v. #Pat1 defined on S_n
#Warning: if k is larger than 2, it is very slow.
#For example, try kthMomentV([1,2],3,n);
kthMomentV:=proc(Pat1,k,n)
if k=0 then
1:
else
Mishkal([Pat1$k],n):
fi:
end:

#kthMoment(Pat1,k,n): The k^th moment ABOUT THE MEAN
#of the r.v. #Pat1 defined on S_n
#Warning: if k is larger than 2, it is very slow.
#For example, try kthMoment([1,2],3,n);
kthMoment:=proc(Pat1,k,n) local r,m,r1:  
r:=nops(Pat1):
m:=expand(binomial(n,r))/r!:
factor(expand(convert([seq((-1)^(k-r1)*binomial(k,r1)*
m^(k-r1)*kthMomentV(Pat1,r1,n),r1=0..k)],`+`))):
end:


 
#Kurtosis(Pat1,n): The  Kurtosis of the r.v. #Pat1 defined on S_n
#For example, try: Kurtosis([1,2],n):
Kurtosis:=proc(Pat1,n):kthMoment(Pat1,4,n)/kthMoment(Pat1,2,n)^2:end:


#kthMoment12Vn0(k0,n0): The k0^th moment of the r.v. #12, on S_n0
#(n0 an integer) a faster way
#(using generating functions) for testing subs(n=n0,kthMomentV([1,2],k,n));
#For example, try: kthMoment12Vn0(k,n0) 
kthMoment12Vn0:=proc(k0,n0)
local gu,q,i,j:
gu:=1:
for i from  1 to n0 do
gu:=expand(gu*normal((1-q^i)/(1-q))):
od:

for j from 1 to k0 do
gu:=expand(q*diff(gu,q)):
od:
subs(q=1,gu)/n0!:
end:

#kthMoment12n0(k0,n0): The k^th moment ABOUT THE MEAN
#of the r.v. #12 defined on S_n0 (n0 an integer)
#For example, try kthMoment12n0([1,2],3,5);
kthMoment12n0:=proc(k,n0) local  m,r1: 
m:=expand(binomial(n0,2))/2!:
expand(convert([seq((-1)^(k-r1)*binomial(k,r1)*
m^(k-r1)*kthMoment12Vn0(r1,n0),r1=0..k)],`+`)):
end:

#Ain(i,n): E(binomial(X,i)) for the r.v. #12 over S_n (symbolic n)
Ain:=proc(i,n) local gu,j,N:
option remember :

if i=0 then
 RETURN(1):
fi:
gu:=0:

for j from 1 to i do
gu:=expand(gu+expand(binomial(n,j+1))/n*expand(subs(n=n-1,Ain(i-j,n)))):
od:
factor(expand(subs(N=n,sum(gu,n=1..N)))):
end:


Bin:=proc(i,n) local k:
if i=0 then
 RETURN(1):
fi:
factor(expand(convert(
[seq(stirling2(i, k)* k!* Ain(k,n), k = 0 .. i)], `+`))):
end:

Cin:=proc(i,n)  local r1,m: m:=expand(binomial(n,2))/2!:
factor(expand(convert([seq((-1)^(i-r1)*binomial(i,r1)*
m^(i-r1)*Bin(r1,n),r1=0..i)],`+`))):
end:

##begin waste-basket
###This should be deleted eventually it is the special case
#of MishkalT2b with Resh1 of length 2
#MishkalT2b2(Resh1): A quick way to compute 
#the (coeff. of the) third-to-leading contribution to Mishkal(Resh1,n)
#i.e. Mishkal(Resh1,n,sums(nops(Resh1[i]),i=1..nops(Resh1))-2)
#So far Resh1 can only have length 2
#For example, try MishkalT2b2([[1,2],[1,2]]);
MishkalT2b2:=proc(Resh1) local pi,sig,a,b,r1,s1,r2,s2,schum:

if nops(Resh1)<>2 then
ERROR(`nops(Resh1) must be 2, the general case is not yet implemented`):
fi:

pi:=Resh1[1]:
sig:=Resh1[2]:
a:=nops(pi):
b:=nops(sig):
schum:=0:
for r1 from 1 to a do
for r2 from r1+1 to a do
 for s1 from 1 to b do
  for s2 from s1+1 to b do
     if (pi[r1]<pi[r2] and sig[s1]<sig[s2]) then
      schum:=schum+binomial(r1+s1-2,r1-1)*
                   binomial(r2+s2-r1-s1-2,r2-r1-1)*
                   binomial(a+b-r2-s2,a-r2)*
          binomial(pi[r1]+sig[s1]-2,pi[r1]-1)*
          binomial(pi[r2]-pi[r1]+sig[s2]-sig[s1]-2,pi[r2]-pi[r1]-1)*
          binomial(a+b-pi[r2]-sig[s2],a-pi[r2]):
      fi:

     if (pi[r2]<pi[r1] and sig[s2]<sig[s1]) then
      schum:=schum+binomial(r1+s1-2,r1-1)*
                   binomial(r2+s2-r1-s1-2,r2-r1-1)*
                   binomial(a+b-r2-s2,a-r2)*
          binomial(pi[r2]+sig[s2]-2,pi[r2]-1)*
          binomial(pi[r1]-pi[r2]+sig[s1]-sig[s2]-2,pi[r1]-pi[r2]-1)*
          binomial(a+b-pi[r1]-sig[s1],a-pi[r1]):
      fi:

od:
od:         
od:
od:

schum: 

end:

MishkalT2:=proc(Resh1):
MishkalT2a(Resh1)+MishkalT2b(Resh1)+MishkalT2c(Resh1):
end:

#Profile3(Triplet): Given a triplet of sets
#outputs [nops(S1 intersect S2),nops(S1 intersect S3),nops(S2 intersect S3)]
#For example, do Profile3([[1,2],[1,3],[2,3]]);
Profile3:=proc(Tri)
local S1,S2,S3:
S1:=convert(Tri[1],set):
S2:=convert(Tri[2],set):
S3:=convert(Tri[3],set):
[nops(S1 intersect S2),nops(S1 intersect S3),nops(S2 intersect S3)]:
end:

ez:=proc():print(`S110(Resh1), SubProf110(Tri)`):
print(`AllProf110(Resh1), Comp(k,n), PotProf110(Resh), Bdok110(Resh)`):
print(`IsCompat(Resh,Pair1), EligProf110(Resh), Bdok110A`):
print(`TW110(Resh), Bdok110B(Resh), MishkalT2c,BdokMish(Resh1):`):
end:


##New Stuff
#S110(Resh1): the set of triplets with profile 110
#(i.e. |S1 intersect S2|=1, |S1 intersect S3|=1)
#and 12 before 13 for example try:S110([[1,2],[1,2],[1,2]]);
S110:=proc(Resh1)
local gu,i,mu,t,mu1,mu2,S1,S2,S3,ku:
option remember:
t:=convert([seq(nops(Resh1[i]),i=1..nops(Resh1))],`+`)-2:
mu:=SetMishkalT(Resh1,t):

gu:={}:

for i from 1 to nops(mu) do
 if Profile3(mu[i][1])=[1,1,0] then
   ku:=mu[i]:
    S1:=convert(ku[1][1],set):
    S2:=convert(ku[1][2],set):
    S3:=convert(ku[1][3],set):
  mu1:=S1 intersect S2:
  mu2:=S1 intersect S3:
 if  mu1[1] <mu2[1] then
     gu:=gu union {mu[i]}:
  
 fi:
 fi:
od:
gu:
end:

#SubProf110(Tri): Given a pair
#[triplet,perm] [[S1,S2,S3],Per] with [12,13] profile
#outputs [[x,y],[x1,y1]] where x,y,x1,y1, are triplets of integers
#For example, try SubProf([[1,2],[1,3],[2,4]],[1,2,3,4]]);
SubProf110:=proc(Tri)
local S,T,S1S2,S1S3,a,b,A,B,i,x,y,x1,y1,pi,i1,tem:
for i from 1 to 3 do
S[i]:=convert(Tri[1][i],set):
od:

pi:=Tri[2]:

S1S2:=S[1] intersect S[2]:
S1S3:=S[1] intersect S[3]:

if nops(S1S2)<>1 or nops(S1S3)<>1 then
 ERROR(`Wrong input`):
fi:

a:=S1S2[1]:b:=S1S3[1]:

A:=pi[a]:
B:=pi[b]:

if A>B then
tem:=A:
A:=B:
B:=tem:
fi:

if not (a<b) then
 ERROR(`Wrong input`):
fi:

for i from 1 to 3 do
T[i]:={seq(pi[S[i][i1]],i1=1..nops(S[i]))}
od:


for i from 1 to 3 do
x[i]:=0: y[i]:=0: x1[i]:=0: y1[i]:=0:
od:

for i from 1 to a-1 do
 if member(i,S[1]) then
    x[1]:=x[1]+1:
 elif member(i,S[2]) then
    x[2]:=x[2]+1:
 else
    x[3]:=x[3]+1:
 fi:
od:

for i from a+1 to b-1 do
 if member(i,S[1]) then
    y[1]:=y[1]+1:
 elif member(i,S[2]) then
    y[2]:=y[2]+1:
 else
    y[3]:=y[3]+1:
 fi:
od:



for i from 1 to A-1 do
 if member(i,T[1]) then
    x1[1]:=x1[1]+1:
 elif member(i,T[2]) then
    x1[2]:=x1[2]+1:
 else
    x1[3]:=x1[3]+1:
 fi:
od:

for i from A+1 to B-1 do
 if member(i,T[1]) then
    y1[1]:=y1[1]+1:
 elif member(i,T[2]) then
    y1[2]:=y1[2]+1:
 else
    y1[3]:=y1[3]+1:
 fi:
od:

[
 [
   [x[1],x[2],x[3]],
   [y[1],y[2],y[3]]
 ]
,
 [
   [x1[1],x1[2],x1[3]],
   [y1[1],y1[2],y1[3]]
 ]
]:

end:

AllProf110:=proc(Resh1) local lu, i:
lu:=S110(Resh1):
{seq(SubProf110(lu[i])[1],i=1..nops(lu))},
{seq(SubProf110(lu[i])[2],i=1..nops(lu))},
{seq(SubProf110(lu[i]),i=1..nops(lu))}
:
end:

#Comp(k,n): All the set of k-vectors of non-negative integers
#that sum to n. Try Comp(2,4);
Comp:=proc(k,n)
local lu,gu,m,i:
option remember:
if k=0 then
 if n=0 then
  RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:
gu:={}:
for m from 0 to n do
 lu:=Comp(k-1,n-m):
 gu:= gu union {seq([op(lu[i]),m],i=1..nops(lu))}:
od:
gu:
end:

#PotProf110(Resh): All the potential first-components of the subprofile
#of a triplet Resh
PotProf110:=proc(Resh)
local a1,a2,a3,lu1,lu2,lu3,i1,i2,i3,gu:
a1:=nops(Resh[1]):
a2:=nops(Resh[2]):
a3:=nops(Resh[3]):
lu1:=Comp(3,a1-2):
lu2:=Comp(3,a2-1):
lu3:=Comp(3,a3-1):
gu:={}:
for i1 from 1 to nops(lu1) do
for i2 from 1 to nops(lu2) do
for i3 from 1 to nops(lu3) do
 gu:=gu union {[[lu1[i1][1],lu2[i2][1],lu3[i3][1]],
[lu1[i1][2],lu2[i2][2],lu3[i3][2]]]}:
od:
od:
od:
gu:
end:

Bdok110:=proc(Resh1):evalb(PotProf110(Resh1)=AllProf110(Resh1)[1]):end:

rev:=proc(pi) local i:[seq(pi[nops(pi)-i+1],i=1..nops(pi))]:end:
Rev:=proc(Resh) local i:[seq(rev(Resh[i]),i=1..nops(Resh))]:end:

#IsCompat(Resh,Pair1): Given a triplet of patterns, and a pair Pair1
#compositions for the first component, finds all the x1y1
#finds whether is is compatible  with them
#For example, try:Compat([[1,2],[1,2],[1,2]],[[0,1,0],[0,0,1]]);
IsCompat:=proc(Resh,Pair1)
local xy,x1y1,x,y,x1,y1:
xy:=Pair1[1]:
x1y1:=Pair1[2]:
x:=xy[1]:
y:=xy[2]:
x1:=x1y1[1]:
y1:=x1y1[2]:


if Resh[1][x[1]+1]<Resh[1][x[1]+y[1]+2] then
 if (x1[1]=Resh[1][x[1]+1]-1 and 
  y1[1]=Resh[1][x[1]+y[1]+2]-Resh[1][x[1]+1]-1 and
  x1[2]=Resh[2][x[2]+1]-1 and
  x1[3]+y1[3]=Resh[3][x[3]+y[3]+1]-1) then
   RETURN(true):
  else
   RETURN(false):
 fi: 
else
 if (x1[1]=Resh[1][x[1]+y[1]+2]-1 and
   y1[1]=Resh[1][x[1]+1]-Resh[1][x[1]+y[1]+2]-1 and
   x1[2]+y1[2]=Resh[2][x[2]+1]-1 and
  x1[3]=Resh[3][x[3]+y[3]+1]-1) then
   RETURN(true):
  else
   RETURN(false):
  fi:
fi:


end:


#EligProf110(Resh): All the Eligible subprofiles
#of a triplet Resh
EligProf110:=proc(Resh)
local lu,i,j,gu,lu1,lu2,Pair1:

lu:=PotProf110(Resh):
gu:={}:

for i from 1 to nops(lu) do
  lu1:=lu[i]:
 
   for j from 1 to nops(lu) do
    lu2:=lu[j]:
    Pair1:=[lu1,lu2]:
     if IsCompat(Resh,Pair1) then
       gu:=gu union  {Pair1}:
     fi:
   od:
od:
gu:

end:

Bdok110A:=proc(Resh1):evalb(EligProf110(Resh1)=AllProf110(Resh1)[3]):end:


#TW110(Resh): the total number of triplets of sets 
#[S1,S2,S3] with two distinct common points (e.g. S1 and S2 intersect
# at a point and S1 and S3 intersect at another point)
#rresponding to the triplet of patterns Resh 
#For example, try TW110([[1,2],[1,2],[1,2]]);
TW110:=proc(Resh)
local lu,i,x,y,z,x1,y1,z1,a1,a2,a3,gu:
option remember:
lu:=EligProf110(Resh):
a1:=nops(Resh[1]):
a2:=nops(Resh[2]):
a3:=nops(Resh[3]):
gu:=0:

for i from 1 to nops(lu) do
 x:=lu[i][1][1]:y:=lu[i][1][2]:
 z:=[a1-2-x[1]-y[1],a2-1-x[2]-y[2],a3-1-x[3]-y[3]]:
 x1:=lu[i][2][1]:y1:=lu[i][2][2]:
 z1:=[a1-2-x1[1]-y1[1],a2-1-x1[2]-y1[2],a3-1-x1[3]-y1[3]]:
 gu:=gu+MTC(op(x))*MTC(op(y))*MTC(op(z))*MTC(op(x1))*MTC(op(y1))*MTC(op(z1)):
od:
gu:
end:


Bdok110B:=proc(Resh1):evalb(nops(S110(Resh1))=TW110(Resh1)):end:

#MishkalT2c(Resh1): A quick way to compute 
#the (coeff. of the) third-to-leading contribution to Mishkal(Resh1,n)
#(coming from interaction between three patterns)
#where one of the members of Resh1 intersects each of the other two
#For example, try MishkalT2c([[1,2],[1,2],[1,2]]);
MishkalT2c:=proc(Resh1) local a,b,c:
a:=Resh1[1]:b:=Resh1[2]:c:=Resh1[3]:
TW110([a,b,c])+TW110([a,c,b])+
TW110([b,a,c])+TW110([b,c,a])+
TW110([c,a,b])+TW110([c,b,a]):
end:

BdokMish:=proc(Resh1) local i,ku:
ku:=convert([seq(nops(Resh1[i]),i=1..nops(Resh1))],`+`):
evalb(MishkalT2(Resh1)=nops(SetMishkalT(Resh1,ku-2))):
end:

#AllAThirdMom(k): The list of patterns of length k
#according to their asymptotic third-moment-about-the mean
# of the r.v. #Pat,(defined on S^n) from smallest
#to biggest (divided by n^(2k-2))
#For example, try AllAThirdMom(3); 
AllAThirdMom:=proc(k)
local gu,T,i,pi,lu,S,n,mu:
option remember:
gu:=Reps1(k):
mu:={}:

for i from 1 to nops(gu) do
pi:=gu[i]:
lu:=AThirdMom(pi,n)/n^(3*k-2):
T[pi]:=lu:
mu:=mu union {lu}:
od:

for i from 1 to nops(mu) do
S[mu[i]]:={}:
od:

for i from 1 to nops(gu) do
pi:=gu[i]:
S[T[pi]]:=S[T[pi]] union {pi}:
od:
mu:=convert(mu,list):
mu:=sort(mu):
[seq([S[mu[i]],mu[i]],i=1..nops(mu))]:
end:





#######Unique Represative for AllAVariance and AllACorr
#revers(perm): The reverse of a permutation
revers:=proc(perm)
local i:
[seq(perm[nops(perm)-i+1],i=1..nops(perm))]:
end:
 
#REVERS(Perms): The set of reverses of a set of permutations
#Perms
 
REVERS:=proc(Perms)
local i,gu:
 
gu:={}:
 
for i from 1 to nops(Perms) do
gu:=gu  union{revers(op(i,Perms))}:
od:
gu:
end:
 
 
#OPPOS(Perms): The set of oposites of a set of permutations
#Perms
 
OPPOS:=proc(Perms)
local i,gu:
 
gu:={}:
 
for i from 1 to nops(Perms) do
gu:=gu  union{oppos(op(i,Perms))}:
od:
gu:
end:
 
 
#oppos(perm): The reverse of a permutation
oppos:=proc(perm)
local i:
[seq(nops(perm)+1-perm[i],i=1..nops(perm))]:
end:
 
 
#khaverim(perm): all the images of the permutation perm
#under the group generated by reversal and opposite
 
khaverim:=proc(perm):
{perm, oppos(perm),revers(perm),oppos(revers(perm))}:
end:
 

INVERSE:=proc(Perms) local i:{seq(hofkhi(Perms[i]),i=1..nops(Perms))}:end:

#KHAVERIM(Perms): all the images of the set of permutations Perm
#under the group generated by (global) Reversal and Opposite
KHAVERIM:=proc(Perms):
#{Perms, OPPOS(Perms),REVERS(Perms),OPPOS(REVERS(Perms))}:
{Perms, REVERS(Perms),INVERSE(Perms),
REVERS(INVERSE(Perms)),INVERSE(REVERS(Perms)),
INVERSE(REVERS(INVERSE(Perms))),
REVERS(INVERSE(REVERS(Perms))),
INVERSE(REVERS(INVERSE(REVERS(Perms))))}:
end:

 

#hofkhi(tau): the inverse of the permutation tau
hofkhi:=proc(tau)
local i,sig1,sig:
 
for i from 1 to nops(tau) do
sig1[tau[i]]:=i:
od:
 
sig:=[]:
 
for i from 1 to nops(tau) do
sig:=[op(sig),sig1[i]]:
od:
sig:
end:
 

#Reps1(n): Distinct Representatives of permutations of length n
Reps1:=proc(n)
local mu,gu:
with(combinat):
gu:=convert(permute(n),set):
mu:={}:

while gu<>{} do
 mu:=mu union {gu[1]}:
 gu:=gu minus khaverim(gu[1]):
od:
mu:
end:

#Reps2(n): Distinct Representatives of pairs of permutations of length n
Reps2:=proc(n)
local mu,gu,i,j:
with(combinat):
gu:=permute(n):
gu:={seq(seq({gu[i],gu[j]},j=i+1..nops(gu)),i=1..nops(gu))}:

mu:={}:

while gu<>{} do
 mu:=mu union {gu[1]}:
 gu:=gu minus KHAVERIM(gu[1]):
od:
mu:
end:

#Reps2DiffSize(n1,n2): 
#Distinct Representatives of pairs of permutations of lengths n1,n2
#For Example try:  Reps2DiffSize(n1,n2);
Reps2DiffSize:=proc(n1,n2)
local mu,gu,gu1,gu2,i,j:
with(combinat):
if n1=n2 then
 RETURN(Reps2(n1)):
fi:
gu1:=permute(n1):
gu2:=permute(n2):
gu:={seq(seq({gu1[i],gu2[j]},j=1..nops(gu2)),i=1..nops(gu1))}:

mu:={}:

while gu<>{} do
 mu:=mu union {gu[1]}:
 gu:=gu minus KHAVERIM(gu[1]):
od:
mu:
end:

#########End Unique Reprs

#AllkthMoments(a,k,n): The list of all kth Moments
#for permutations of size a
#For example, try AllkthMoments(3,2);
AllkthMoments:=proc(a,k,n)
local gu,i:
option remember:
gu:=Reps1(a):

for i from 1 to nops(gu) do
print(`The exact (not asymptotic) expression for the`, k, `-th moment`):
print(`of the Random Variable on S_n, #of occurences of the pattern`):
print( gu[i], `and its images is`):
print(kthMoment(gu[i],k,n)):
od:
end:



#Kurtosis(Pat1,n): The  Kurtosis of the r.v. #Pat1 defined on S_n
#For example, try: Kurtosis([1,2],n):
Kurtosis:=proc(Pat1,n):kthMoment(Pat1,4,n)/kthMoment(Pat1,2,n)^2:end:


#kthMoment12Vn0(k0,n0): The k0^th moment of the r.v. #12, on S_n0
#(n0 an integer) a faster way
#(using generating functions) for testing subs(n=n0,kthMomentV([1,2],k,n));
#For example, try: kthMoment12Vn0(k,n0) 
kthMoment12Vn0:=proc(k0,n0)
local gu,q,i,j:
gu:=1:
for i from  1 to n0 do
gu:=expand(gu*normal((1-q^i)/(1-q))):
od:

for j from 1 to k0 do
gu:=expand(q*diff(gu,q)):
od:
subs(q=1,gu)/n0!:
end:

#kthMoment12n0(k0,n0): The k^th moment ABOUT THE MEAN
#of the r.v. #12 defined on S_n0 (n0 an integer)
#For example, try kthMoment12n0([1,2],3,5);
kthMoment12n0:=proc(k,n0) local  m,r1: 
m:=expand(binomial(n0,2))/2!:
expand(convert([seq((-1)^(k-r1)*binomial(k,r1)*
m^(k-r1)*kthMoment12Vn0(r1,n0),r1=0..k)],`+`)):
end:

#Ain(i,n): E(binomial(X,i)) for the r.v. #12 over S_n (symbolic n)
Ain:=proc(i,n) local gu,j,N:
option remember :

if i=0 then
 RETURN(1):
fi:
gu:=0:

for j from 1 to i do
gu:=expand(gu+expand(binomial(n,j+1))/n*expand(subs(n=n-1,Ain(i-j,n)))):
od:
factor(expand(subs(N=n,sum(gu,n=1..N)))):
end:


Bin:=proc(i,n) local k:
if i=0 then
 RETURN(1):
fi:
factor(expand(convert(
[seq(stirling2(i, k)* k!* Ain(k,n), k = 0 .. i)], `+`))):
end:

Cin:=proc(i,n)  local r1,m: m:=expand(binomial(n,2))/2!:
factor(expand(convert([seq((-1)^(i-r1)*binomial(i,r1)*
m^(i-r1)*Bin(r1,n),r1=0..i)],`+`))):
end:

kthMoment12:=proc(i,n):Cin(i,n):end:

##begin waste-basket
###This should be deleted eventually it is the special case
#of MishkalT2b with Resh1 of length 2
#MishkalT2b2(Resh1): A quick way to compute 
#the (coeff. of the) third-to-leading contribution to Mishkal(Resh1,n)
#i.e. Mishkal(Resh1,n,sums(nops(Resh1[i]),i=1..nops(Resh1))-2)
#So far Resh1 can only have length 2
#For example, try MishkalT2b2([[1,2],[1,2]]);
MishkalT2b2:=proc(Resh1) local pi,sig,a,b,r1,s1,r2,s2,schum:

if nops(Resh1)<>2 then
ERROR(`nops(Resh1) must be 2, the general case is not yet implemented`):
fi:

pi:=Resh1[1]:
sig:=Resh1[2]:
a:=nops(pi):
b:=nops(sig):
schum:=0:
for r1 from 1 to a do
for r2 from r1+1 to a do
 for s1 from 1 to b do
  for s2 from s1+1 to b do
     if (pi[r1]<pi[r2] and sig[s1]<sig[s2]) then
      schum:=schum+binomial(r1+s1-2,r1-1)*
                   binomial(r2+s2-r1-s1-2,r2-r1-1)*
                   binomial(a+b-r2-s2,a-r2)*
          binomial(pi[r1]+sig[s1]-2,pi[r1]-1)*
          binomial(pi[r2]-pi[r1]+sig[s2]-sig[s1]-2,pi[r2]-pi[r1]-1)*
          binomial(a+b-pi[r2]-sig[s2],a-pi[r2]):
      fi:

     if (pi[r2]<pi[r1] and sig[s2]<sig[s1]) then
      schum:=schum+binomial(r1+s1-2,r1-1)*
                   binomial(r2+s2-r1-s1-2,r2-r1-1)*
                   binomial(a+b-r2-s2,a-r2)*
          binomial(pi[r2]+sig[s2]-2,pi[r2]-1)*
          binomial(pi[r1]-pi[r2]+sig[s1]-sig[s2]-2,pi[r1]-pi[r2]-1)*
          binomial(a+b-pi[r1]-sig[s1],a-pi[r1]):
      fi:

od:
od:         
od:
od:

schum: 

end:

MishkalT2:=proc(Resh1):
MishkalT2a(Resh1)+MishkalT2b(Resh1)+MishkalT2c(Resh1):
end:

#Profile3(Triplet): Given a triplet of sets
#outputs [nops(S1 intersect S2),nops(S1 intersect S3),nops(S2 intersect S3)]
#For example, do Profile3([[1,2],[1,3],[2,3]]);
Profile3:=proc(Tri)
local S1,S2,S3:
S1:=convert(Tri[1],set):
S2:=convert(Tri[2],set):
S3:=convert(Tri[3],set):
[nops(S1 intersect S2),nops(S1 intersect S3),nops(S2 intersect S3)]:
end:

ez:=proc():print(`S110(Resh1), SubProf110(Tri)`):
print(`AllProf110(Resh1), Comp(k,n), PotProf110(Resh), Bdok110(Resh)`):
print(`IsCompat(Resh,Pair1), EligProf110(Resh), Bdok110A`):
print(`TW110(Resh), Bdok110B(Resh), MishkalT2c,BdokMish(Resh1):`):
end:


##New Stuff
#S110(Resh1): the set of triplets with profile 110
#(i.e. |S1 intersect S2|=1, |S1 intersect S3|=1)
#and 12 before 13 for example try:S110([[1,2],[1,2],[1,2]]);
S110:=proc(Resh1)
local gu,i,mu,t,mu1,mu2,S1,S2,S3,ku:
option remember:
t:=convert([seq(nops(Resh1[i]),i=1..nops(Resh1))],`+`)-2:
mu:=SetMishkalT(Resh1,t):

gu:={}:

for i from 1 to nops(mu) do
 if Profile3(mu[i][1])=[1,1,0] then
   ku:=mu[i]:
    S1:=convert(ku[1][1],set):
    S2:=convert(ku[1][2],set):
    S3:=convert(ku[1][3],set):
  mu1:=S1 intersect S2:
  mu2:=S1 intersect S3:
 if  mu1[1] <mu2[1] then
     gu:=gu union {mu[i]}:
  
 fi:
 fi:
od:
gu:
end:

#SubProf110(Tri): Given a pair
#[triplet,perm] [[S1,S2,S3],Per] with [12,13] profile
#outputs [[x,y],[x1,y1]] where x,y,x1,y1, are triplets of integers
#For example, try SubProf([[1,2],[1,3],[2,4]],[1,2,3,4]]);
SubProf110:=proc(Tri)
local S,T,S1S2,S1S3,a,b,A,B,i,x,y,x1,y1,pi,i1,tem:
for i from 1 to 3 do
S[i]:=convert(Tri[1][i],set):
od:

pi:=Tri[2]:

S1S2:=S[1] intersect S[2]:
S1S3:=S[1] intersect S[3]:

if nops(S1S2)<>1 or nops(S1S3)<>1 then
 ERROR(`Wrong input`):
fi:

a:=S1S2[1]:b:=S1S3[1]:

A:=pi[a]:
B:=pi[b]:

if A>B then
tem:=A:
A:=B:
B:=tem:
fi:

if not (a<b) then
 ERROR(`Wrong input`):
fi:

for i from 1 to 3 do
T[i]:={seq(pi[S[i][i1]],i1=1..nops(S[i]))}
od:


for i from 1 to 3 do
x[i]:=0: y[i]:=0: x1[i]:=0: y1[i]:=0:
od:

for i from 1 to a-1 do
 if member(i,S[1]) then
    x[1]:=x[1]+1:
 elif member(i,S[2]) then
    x[2]:=x[2]+1:
 else
    x[3]:=x[3]+1:
 fi:
od:

for i from a+1 to b-1 do
 if member(i,S[1]) then
    y[1]:=y[1]+1:
 elif member(i,S[2]) then
    y[2]:=y[2]+1:
 else
    y[3]:=y[3]+1:
 fi:
od:



for i from 1 to A-1 do
 if member(i,T[1]) then
    x1[1]:=x1[1]+1:
 elif member(i,T[2]) then
    x1[2]:=x1[2]+1:
 else
    x1[3]:=x1[3]+1:
 fi:
od:

for i from A+1 to B-1 do
 if member(i,T[1]) then
    y1[1]:=y1[1]+1:
 elif member(i,T[2]) then
    y1[2]:=y1[2]+1:
 else
    y1[3]:=y1[3]+1:
 fi:
od:

[
 [
   [x[1],x[2],x[3]],
   [y[1],y[2],y[3]]
 ]
,
 [
   [x1[1],x1[2],x1[3]],
   [y1[1],y1[2],y1[3]]
 ]
]:

end:

AllProf110:=proc(Resh1) local lu, i:
lu:=S110(Resh1):
{seq(SubProf110(lu[i])[1],i=1..nops(lu))},
{seq(SubProf110(lu[i])[2],i=1..nops(lu))},
{seq(SubProf110(lu[i]),i=1..nops(lu))}
:
end:

#Comp(k,n): All the set of k-vectors of non-negative integers
#that sum to n. Try Comp(2,4);
Comp:=proc(k,n)
local lu,gu,m,i:
option remember:
if k=0 then
 if n=0 then
  RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:
gu:={}:
for m from 0 to n do
 lu:=Comp(k-1,n-m):
 gu:= gu union {seq([op(lu[i]),m],i=1..nops(lu))}:
od:
gu:
end:

#PotProf110(Resh): All the potential first-components of the subprofile
#of a triplet Resh
PotProf110:=proc(Resh)
local a1,a2,a3,lu1,lu2,lu3,i1,i2,i3,gu:
a1:=nops(Resh[1]):
a2:=nops(Resh[2]):
a3:=nops(Resh[3]):
lu1:=Comp(3,a1-2):
lu2:=Comp(3,a2-1):
lu3:=Comp(3,a3-1):
gu:={}:
for i1 from 1 to nops(lu1) do
for i2 from 1 to nops(lu2) do
for i3 from 1 to nops(lu3) do
 gu:=gu union {[[lu1[i1][1],lu2[i2][1],lu3[i3][1]],
[lu1[i1][2],lu2[i2][2],lu3[i3][2]]]}:
od:
od:
od:
gu:
end:

Bdok110:=proc(Resh1):evalb(PotProf110(Resh1)=AllProf110(Resh1)[1]):end:

rev:=proc(pi) local i:[seq(pi[nops(pi)-i+1],i=1..nops(pi))]:end:
Rev:=proc(Resh) local i:[seq(rev(Resh[i]),i=1..nops(Resh))]:end:

#IsCompat(Resh,Pair1): Given a triplet of patterns, and a pair Pair1
#compositions for the first component, finds all the x1y1
#finds whether is is compatible  with them
#For example, try:Compat([[1,2],[1,2],[1,2]],[[0,1,0],[0,0,1]]);
IsCompat:=proc(Resh,Pair1)
local xy,x1y1,x,y,x1,y1:
xy:=Pair1[1]:
x1y1:=Pair1[2]:
x:=xy[1]:
y:=xy[2]:
x1:=x1y1[1]:
y1:=x1y1[2]:


if Resh[1][x[1]+1]<Resh[1][x[1]+y[1]+2] then
 if (x1[1]=Resh[1][x[1]+1]-1 and 
  y1[1]=Resh[1][x[1]+y[1]+2]-Resh[1][x[1]+1]-1 and
  x1[2]=Resh[2][x[2]+1]-1 and
  x1[3]+y1[3]=Resh[3][x[3]+y[3]+1]-1) then
   RETURN(true):
  else
   RETURN(false):
 fi: 
else
 if (x1[1]=Resh[1][x[1]+y[1]+2]-1 and
   y1[1]=Resh[1][x[1]+1]-Resh[1][x[1]+y[1]+2]-1 and
   x1[2]+y1[2]=Resh[2][x[2]+1]-1 and
  x1[3]=Resh[3][x[3]+y[3]+1]-1) then
   RETURN(true):
  else
   RETURN(false):
  fi:
fi:


end:


#EligProf110(Resh): All the Eligible subprofiles
#of a triplet Resh
EligProf110:=proc(Resh)
local lu,i,j,gu,lu1,lu2,Pair1:

lu:=PotProf110(Resh):
gu:={}:

for i from 1 to nops(lu) do
  lu1:=lu[i]:
 
   for j from 1 to nops(lu) do
    lu2:=lu[j]:
    Pair1:=[lu1,lu2]:
     if IsCompat(Resh,Pair1) then
       gu:=gu union  {Pair1}:
     fi:
   od:
od:
gu:

end:

Bdok110A:=proc(Resh1):evalb(EligProf110(Resh1)=AllProf110(Resh1)[3]):end:


#TW110(Resh): the total number of triplets of sets 
#[S1,S2,S3] with two distinct common points (e.g. S1 and S2 intersect
# at a point and S1 and S3 intersect at another point)
#rresponding to the triplet of patterns Resh 
#For example, try TW110([[1,2],[1,2],[1,2]]);
TW110:=proc(Resh)
local lu,i,x,y,z,x1,y1,z1,a1,a2,a3,gu:
option remember:
lu:=EligProf110(Resh):
a1:=nops(Resh[1]):
a2:=nops(Resh[2]):
a3:=nops(Resh[3]):
gu:=0:

for i from 1 to nops(lu) do
 x:=lu[i][1][1]:y:=lu[i][1][2]:
 z:=[a1-2-x[1]-y[1],a2-1-x[2]-y[2],a3-1-x[3]-y[3]]:
 x1:=lu[i][2][1]:y1:=lu[i][2][2]:
 z1:=[a1-2-x1[1]-y1[1],a2-1-x1[2]-y1[2],a3-1-x1[3]-y1[3]]:
 gu:=gu+MTC(op(x))*MTC(op(y))*MTC(op(z))*MTC(op(x1))*MTC(op(y1))*MTC(op(z1)):
od:
gu:
end:


Bdok110B:=proc(Resh1):evalb(nops(S110(Resh1))=TW110(Resh1)):end:

#MishkalT2c(Resh1): A quick way to compute 
#the (coeff. of the) third-to-leading contribution to Mishkal(Resh1,n)
#(coming from interaction between three patterns)
#where one of the members of Resh1 intersects each of the other two
#For example, try MishkalT2c([[1,2],[1,2],[1,2]]);
MishkalT2c:=proc(Resh1) local a,b,c:
a:=Resh1[1]:b:=Resh1[2]:c:=Resh1[3]:
TW110([a,b,c])+TW110([a,c,b])+
TW110([b,a,c])+TW110([b,c,a])+
TW110([c,a,b])+TW110([c,b,a]):
end:

BdokMish:=proc(Resh1) local i,ku:
ku:=convert([seq(nops(Resh1[i]),i=1..nops(Resh1))],`+`):
evalb(MishkalT2(Resh1)=nops(SetMishkalT(Resh1,ku-2))):
end:




#######Unique Represative for AllAVariance and AllACorr
#revers(perm): The reverse of a permutation
revers:=proc(perm)
local i:
[seq(perm[nops(perm)-i+1],i=1..nops(perm))]:
end:
 
#REVERS(Perms): The set of reverses of a set of permutations
#Perms
 
REVERS:=proc(Perms)
local i,gu:
 
gu:={}:
 
for i from 1 to nops(Perms) do
gu:=gu  union{revers(op(i,Perms))}:
od:
gu:
end:
 
 
#OPPOS(Perms): The set of oposites of a set of permutations
#Perms
 
OPPOS:=proc(Perms)
local i,gu:
 
gu:={}:
 
for i from 1 to nops(Perms) do
gu:=gu  union{oppos(op(i,Perms))}:
od:
gu:
end:
 
 
#oppos(perm): The reverse of a permutation
oppos:=proc(perm)
local i:
[seq(nops(perm)+1-perm[i],i=1..nops(perm))]:
end:
 
 
#khaverim(perm): all the images of the permutation perm
#under the group generated by reversal and opposite
 
khaverim:=proc(perm):
{perm, oppos(perm),revers(perm),oppos(revers(perm))}:
end:
 

INVERSE:=proc(Perms) local i:{seq(hofkhi(Perms[i]),i=1..nops(Perms))}:end:

#KHAVERIM(Perms): all the images of the set of permutations Perm
#under the group generated by (global) Reversal and Opposite
KHAVERIM:=proc(Perms):
#{Perms, OPPOS(Perms),REVERS(Perms),OPPOS(REVERS(Perms))}:
{Perms, REVERS(Perms),INVERSE(Perms),
REVERS(INVERSE(Perms)),INVERSE(REVERS(Perms)),
INVERSE(REVERS(INVERSE(Perms))),
REVERS(INVERSE(REVERS(Perms))),
INVERSE(REVERS(INVERSE(REVERS(Perms))))}:
end:

 

#hofkhi(tau): the inverse of the permutation tau
hofkhi:=proc(tau)
local i,sig1,sig:
 
for i from 1 to nops(tau) do
sig1[tau[i]]:=i:
od:
 
sig:=[]:
 
for i from 1 to nops(tau) do
sig:=[op(sig),sig1[i]]:
od:
sig:
end:
 

#Reps1(n): Distinct Representatives of permutations of length n
Reps1:=proc(n)
local mu,gu:
with(combinat):
gu:=convert(permute(n),set):
mu:={}:

while gu<>{} do
 mu:=mu union {gu[1]}:
 gu:=gu minus khaverim(gu[1]):
od:
mu:
end:

#Reps2(n): Distinct Representatives of pairs of permutations of length n
Reps2:=proc(n)
local mu,gu,i,j:
with(combinat):
gu:=permute(n):
gu:={seq(seq({gu[i],gu[j]},j=i+1..nops(gu)),i=1..nops(gu))}:

mu:={}:

while gu<>{} do
 mu:=mu union {gu[1]}:
 gu:=gu minus KHAVERIM(gu[1]):
od:
mu:
end:

#Reps2DiffSize(n1,n2): 
#Distinct Representatives of pairs of permutations of lengths n1,n2
#For Example try:  Reps2DiffSize(n1,n2);
Reps2DiffSize:=proc(n1,n2)
local mu,gu,gu1,gu2,i,j:
with(combinat):
if n1=n2 then
 RETURN(Reps2(n1)):
fi:
gu1:=permute(n1):
gu2:=permute(n2):
gu:={seq(seq({gu1[i],gu2[j]},j=1..nops(gu2)),i=1..nops(gu1))}:

mu:={}:

while gu<>{} do
 mu:=mu union {gu[1]}:
 gu:=gu minus KHAVERIM(gu[1]):
od:
mu:
end:

#########End Unique Reprs

#AllkthMoments(a,k,n): The list of all kth Moments
#for permutations of size a
#For example, try AllkthMoments(3,2);
AllkthMoments:=proc(a,k,n)
local gu,i:
option remember:
gu:=Reps1(a):

for i from 1 to nops(gu) do
print(`The exact (not asymptotic) expression for the`, k, `-th moment`):
print(`of the Random Variable on S_n, #of occurences of the pattern`):
print( gu[i], `and its images is`):
print(kthMoment(gu[i],k,n)):
od:
end:
 


#kthMomentE(Pat1,k,n): The k^th moment ABOUT THE MEAN
#of the r.v. #Pat1 defined on S_n, using the Empirical 
#(yet rigorous!) method
#For example, try kthMomentE([1,2],2,n);
kthMomentE:=proc(Pat1,k,n):  
EListES([Pat1$k],n):
end:

#AllkthMomentsE(a,k,n): The list of all kth Moments (with empirical method)
#for permutations of size a
#For example, try AllkthMomentsE(3,2);
AllkthMomentsE:=proc(a,k,n)
local gu,i:
option remember:
gu:=Reps1(a):

for i from 1 to nops(gu) do
print(`The exact (not asymptotic) expression for the`, k, `-th moment`):
print(`of the Random Variable on S_n, #of occurences of the pattern`):
print( gu[i], `and its images is`):
print(kthMomentE(gu[i],k,n)):
od:
end:


#AllAThirdMomV(k): A Verbose version of
#AllAThirdMomV(k) (q.v.)
AllAThirdMomV:=proc(k)
local gu,i:
gu:=AllAThirdMom(k):

print(`There are`, nops(gu), `different asymptotic third moments`):
print(` (about the mean) for`):
print(`the random variable #Occurrences of the pattern pi, of length`, k):
print(`defined over all permutations of length n`):
print(`The ratio between the biggest and the smallest is`):
print(evalf(gu[nops(gu)][2]/gu[1][2])):

for i from 1 to nops(gu) do
print(`The `, i, `-th smallest Asymptotic Third Moment belongs to the `):
print(`permutations in the set (and their obvious images)`):
print(gu[i][1]):
print(`whose asymptotic third moment is  n to the power`, (2*k-2), ` times ` ):
print(gu[i][2], `which is roughly`, evalf(gu[i][2])):
od:
end: