######################################################################
##CompositionsPlus.txt: Save this file as                            #
#CompositionsPlus.txt                                                #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read `CompositionsPlus.txt`<Enter>                 #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Jan. 11, 2019
 
print(`Created:  Jan. 11, 2019`):
print(` This is CompositionsPlus.txt `):
print(`It is one of the Maple packages that accompany the article `):
print(`  The "Monkey Typing Shakespeare" Problem for Compositions`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
print(`published in the Personal Journal of Shalosh B. Ekhad and Doron Zeilberger and `):
print(` arxiv.org `):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://sites.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/comp.html  .`):


print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Story procedures for the ORIGINAL case type ezraStory();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Story procedures for the GENERALIZED case type ezraXStory();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the  MAIN procedures about  for the GENERALIZED case type ezraX(); `):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the  procedures about  the MAIN GENERALIZED procedures, type ezraX(); `):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the  procedures about Checking the GENERALIZED procedures, type ezraXch(); `):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


print(`---------------------------------------`):
 print(`For a list of the Supporting procedures of the GENERALIZED case type ezraX1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Multi-Composition procedures type ezraM();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


print(`---------------------------------------`):
 print(`For a list of the Cluster procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):


print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezraXch:=proc()

if args=NULL then
 print(` The Checking procedures for the generalized procedure GFsetX are`):
 print(` CheckGFsetX, CheckGFsetX1, CheckGFsetX11,  Mishkal, NuSin, WtCd`):
else
ezra(args):
fi:

end:


ezra1X:=proc()

if args=NULL then
 print(` The Supporint GENERALIZED procedures (keeping track of the number of offenders) are: AbaSetX,  InitialWtX, InitialWtStX, Kesem, KesemA  `):

else
ezra(args):
fi:

end:

ezraX:=proc()

if args=NULL then
 print(` The MAIN GENERALIZED procedures (keeping track of the number of offenders) are: GFsetX,  `):
 print(`MomsK, MomsKA, MomsKA2, MomsKA2st, MomsKAst, MomsK2 `):

else
ezra(args):
fi:

end:

ezraM:=proc()

if args=NULL then
 print(` The multi-compositions procedures are:  AbaSet, AvotSet, GFset, GoodCseqD, InitialStatesS, InitialWt, InitialWtSt, Reduce1, Reduce11`):
 print(`StatesCset  `):

else
ezra(args):
fi:

end:

ezraC:=proc()

if args=NULL then
 print(` The the cluster procedures are: Aba, Avot, Clus, Clus1, MaxC, StatesC `):
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AsyGF, Av,  CheckGFone, CheckGFset, CheckGFset1, CheckGFset11, Comps, Comps1, IsUni, `):
 print(` GFone, GFoneU, GoodC, GoodC1, GoodCnaive, GoodC1naive, GoodCseqD,  IsCon, IsLess1, ND2, ND2m, ND2mTable `):
 print(` Ove, Ove1, Skyline, Taur, Wt `):
 print(``):

else
ezra(args):
fi:

end:


ezraXStory:=proc()

if args=NULL then
 print(` The story procedure for the GENERALIZED case are`):
print(`The main procedures are:   InfoX, InfoX2V,InfoX2Vth, InfoXV, SuperInfo1X, SuperInfo1XV, SuperInfoX2Vone , SuperInfoX2Vtwo  `):
 print(``):

else
ezra(args):
fi:

end:

ezraStory:=proc()

if args=NULL then
 print(` The story procedure for the original case are`):
print(`The main procedures are:   Info,  InfoV,  SuperInfo1, SuperInfo1V,  SuperInfo2kr, SuperInfo2krV   `):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The MAIN (simple) procedure is:    GFset    `):
 print(` `):


elif nops([args])=1 and op(1,[args])=Aba then
print(`Aba(C,S,r,x,t): inputs a composition C, a state S, and r between 2 and nops(C), outputs the new state and the`):
print(`change in weight in terms of x and t. Try: `):
print(`Aba([2,1,2,1,2],[2,1,2,2],4,x,t);`):


elif nops([args])=1 and op(1,[args])=AbaSet then
print(`AbaSet(SetC1,S,r,x,t): inputs a set of compositions SetC1, a state S, and r between 2 and nops(Skyline(SetC)), outputs the new state and the`):
print(`change in weight in terms of x and t. Try: `):
print(`AbaSet({[2,1,2,1,2]},[2,1,2,2],4,x,t);`):

elif nops([args])=1 and op(1,[args])=AbaSetX then
print(`AbaSetX(SetC1,S,r,x,t,X): inputs a set of compositions SetC1, a state S, and r between 2 and nops(Skyline(SetC)), outputs the new state and the`):
print(`change in weight in terms of x and t for the GENERALIZED case wher we keep track of occurences of offenders. Try: `):
print(`AbaSetX{[2,1,2,1,2]},[2,1,2,2],4,x,t,X);`):


elif nops([args])=1 and op(1,[args])=AsyGF then
print(`AsyGF(f,x): inputs a rational function f and finds the pair [alpha,A] such that the coefficient of x^n in the Taylor expansion of f`):
print(`is A*alpha^n. Try`):
print(`AsyGF(1/(1-x-x^2),x);`):

elif nops([args])=1 and op(1,[args])=Av then
print(`Av(C,SetC): does the composition C avoid the set of compositions SetC. Try:`):
print(`Av([3,4,2],{[2,5],[4,1]});`):

elif nops([args])=1 and op(1,[args])=Avot then
print(`Avot(C,S): all the fathers of the state S with composition C. Try:`):
print(`Avot([2,1,2,1,2],[2,1,2,2]);`):

elif nops([args])=1 and op(1,[args])=AvotSet then
print(`AvotSet(SetC,S): all the fathers of the state S with a set of compositions SetC. Try:`):
print(`AvotSet({[2,1,2,1,2]},[2,1,2,2]);`):

elif nops([args])=1 and op(1,[args])=CheckGFone then
print(`CheckGFone(C,N): checks GFone(C,x) up to N terms. Try:`):
print(`CheckGFone([2,1,2],12);`):

elif nops([args])=1 and op(1,[args])=CheckGFset then
print(`CheckGFset(R,N): checks GFset(SetC,x) up to N terms for all subsets SetC of the set of compositions of <=r with the same number of  parts. Try:`):
print(`CheckGFset(3,10);`):

elif nops([args])=1 and op(1,[args])=CheckGFset1 then
print(`CheckGFset1(r,k,N): checks GFset(SetC,x) up to N terms for all subsets SetC of the set of compositions of r with k parts. Try:`):
print(`CheckGFset1(3,2,10);`):

elif nops([args])=1 and op(1,[args])=CheckGFset11 then
print(`CheckGFset11(SetC,N): checks GFset(SetC,x) up to N terms. Try:`):
print(`CheckGFset11({[2,1,2]},12);`):

elif nops([args])=1 and op(1,[args])=CheckGFsetX then
print(`CheckGFsetX(R,N): checks GFsetX(SetC,x) up to N terms for all subsets SetC of the set of compositions of <=r with the same number of  parts. `):
print(`Warning; very slow. Try:`):
print(`CheckGFsetX(3,10);`):

elif nops([args])=1 and op(1,[args])=CheckGFsetX1 then
print(`CheckGFsetX1(r,k,N): checks GFsetX(SetC,x) up to N terms for all subsets SetC of the set of compositions of r with k parts. Try:`):
print(`CheckGFsetX1(3,2,10);`):



elif nops([args])=1 and op(1,[args])=CheckGFsetX11 then
print(`CheckGFsetX11(SetC,N): checks procedure GFsetX(SetC,x,X) up to the coefficient of x^N. Try:`):
print(`CheckGFsetX11({[2,3,2],[3,2,3]},10);`):

elif nops([args])=1 and op(1,[args])=Clus then
print(`Clus(C,k): The set of clusters with one composition C of height k. Try: `):
print(`Clus([2,1,2],3);`):

elif nops([args])=1 and op(1,[args])=Clus1 then
print(`Clus1(C,k,m): The set of k by m clusters with one composition C. Try: `):
print(`Clus1([2,1,2],3,6);`):



elif nops([args])=1 and op(1,[args])=Comps then
print(`Comps(n,k): the set of compositions of n Try:`):
print(`Comps(5);`):

elif nops([args])=1 and op(1,[args])=Comps1 then
print(`Comps1(n,k): the set of compositions of n into exactly k parts. Try:`):
print(`Comps1(5,2);`):

elif nops([args])=1 and op(1,[args])=GFset then
print(`GFset(S,x): inputs a set of compositions S and outputs the generating function enumerating compositions that avoid`):
print(`the members of S as sub-compositions. Try:`):
print(`GFset({[2,3,2],[4,1,4]},x);`):

elif nops([args])=1 and op(1,[args])=GFsetX then
print(`GFsetX(SetC,x,X): inputs a set of compositions SetC, abd variable names x, and X,`):
print(` and outputs the generating function enumerating compositions accodring to the weight`):
print(` Product(X[C]^(Number of occurrences of C as a subcomposition over all C in SetC`):
print(` Try: `):
print(` GFsetX({[2,3,2],[4,1,4]},x,X); `):

elif nops([args])=1 and op(1,[args])=GoodC1 then
print(`GoodC1(n,k,SetC): The set of compositions of n into k partsavooding the set of compositions SetC as subcompositions. Try:`):
print(`GoodC1(6,3,{[1,2,1]});`):

elif nops([args])=1 and op(1,[args])=GoodC then
print(`GoodC(n,SetC): The set of compositions of n avooding the set of compositions SetC as subcompositions.Try:`):
print(`GoodC(6,{[1,2,1]});`):

elif nops([args])=1 and op(1,[args])=GoodCnaive then
print(`GoodCnaive(n,SetC): The set of compositions of n avooding the set of compositions SetC as subcompositions. Naive approach. Try:`):
print(`GoodCnaive(6,{[1,2,1]});`):

elif nops([args])=1 and op(1,[args])=GoodC1naive then
print(`GoodC1naive(n,k,SetC): The set of compositions of n into k partsavooding the set of compositions SetC as subcompositions. Naive approach. Try:`):
print(`GoodC1naive(6,3,{[1,2,1]});`):

elif nops([args])=1 and op(1,[args])=GFone then
print(`GFone(C,x): the generating function of the sequence enumerating compositions of n avoiding the sub-composition C. Try`):
print(`GFone([2,1,2],x);`):

elif nops([args])=1 and op(1,[args])=GFoneU then
print(`GFoneU(C,x): inputs a UNIMODAL composition  C, and a variable name, x, and outputs the generating function`):
print(`of the enumerating sequence of compositions of n avooding C as a sub-composition. Try:`):
print(`GFoneU([2,3,2],x);`):

elif nops([args])=1 and op(1,[args])=GoodCseqD then
print(`GoodCseqD(N,SetC): the first N+1 terms, starting at n=0 and ending with n=N,`):
print(` of the enumerating sequence for compositions of n avoiding as sub-compositions the`):
print(`members of the set SetC, done DIRECTLY (by counting). For checking only! Try:`):
print(`GoodCseqD(10,{[2,3,2]});`):


elif nops([args])=1 and op(1,[args])=Info then
print(`Info(SetC,x,N): inputs a set of compositions of the same length, SetC, and a variable x, and a positive integer N, returns`):
print(`(i) The input set , SetC`):
print(`(ii) The generating function of the seqeunce enumerating compositions that avoid the compositions in SetC. let's call it a(n)`):
print(`(iii): The first N terms of the enumerating sequence.`):
print(`(iv) the limit of a(n+1)/a(n), as n goes to infinity, let's call it alpha`):
print(`(v) The limit of a(n)/alpha^n as n goes to infinity. Try:`):
print(`Info({[2,3,2]},x,30);`):

elif nops([args])=1 and op(1,[args])=InfoX then
print(`InfoX(SetC,x,X,N): inputs a set of compositions of the same length, SetC, and variables x,X, and a positive integer N, returns`):
print(`(i) The input set , SetC`):
print(`(ii) The generating function of the seqeunce enumerating compositions that avoid the compositions in SetC. let's call it a(n)`):
print(`(ii') The FULL generating function of the seqeunce enumerating compositions according to occurrences of the members of SetC.`):
print(`(iii): The first N terms of the enumerating sequence.`):
print(`(iv) the limit of a(n+1)/a(n), as n goes to infinity, let's call it alpha`):
print(`(v) The limit of a(n)/alpha^n as n goes to infinity. Try:`):
print(`InfoX({[2,3,2]},x,X,30);`):



elif nops([args])=1 and op(1,[args])=InfoX2V then
print(`InfoX2V(C1,C2,x,X1,X2,n,r): inputs two compositions C1 and C2 none a sub-composition of the other, and outputs`):
print(`an article about `):
print(`(i) the generating function, R0(x) enumerating compositions of the sequence enumerating compositions of`):
print(`n, AVOIDING both C1 and C2 (gotten from GFset({C1,C2},x))`):
print(`(ii) the asymptoic growth rate of that enumerating sequence enumerating `):
print(`(iii) The Full generating function, the rational function in three variables R(x,X1,X2) where the coefficient of x^n*X1^a1*x2^a2`):
print(`is the number of compositions of n with exactly a1 occurrences (as containments) of C1 and a2 occurrences (as containments) of C2`):
print(`The expectation and  variance, in terms of n, of the random variables "number of occurrences" of C1 and C2 respectively`):
print(`defined on the set of compositions of n`):
print(`(iv) The correlation of these two random variables`):
print(`(v) A confirmation that they are indeed joint asymptotically normal with that correlation. For example, try:`):
print(`InfoX2V([2,3,2],[3,2,3],x,X,Y,n,4); `):

elif nops([args])=1 and op(1,[args])=InfoX2Vth then
print(`InfoX2Vth(C1,C2,x,X1,X2,n,r,NUM): Like InfoX2V(C1,C2,x,X1,X2,n,r) (q.v.) but with the theorem numbered NUM. Try:`):
print(`InfoX2Vth([2,3,2],[3,2,3],x,X,Y,n,4,11); `):

elif nops([args])=1 and op(1,[args])=InfoXV then
print(`InfoXV(SetC,x,X,N): verbose form of InfoX(SetC,x,X,N) (q.v.). Try: `):
print(`InfoXV({[2,3,2]},x,X,30);`):

elif nops([args])=1 and op(1,[args])=Info1X then
print(`Info1X(C,x,X,n,r): inputs a  composition C, and a variable x, and a variable name X, a variable name n, and a positive`):
print(`integer r. Outputs`):
print(`(i) The generating function ,  a rational function of x and X[op(C)] whose coefficient of x^n is`):
print(`the weight enumerator of the set of compositions of n according to the weight X[op(C)]^NumberOfTimesCisContainedInIt`):
print(`(ii) the  list consisting of the asympotic average, variance, and higher moments (about the mean) modulo peanuts.`):
print(`Try:`):
print(`Info1X([2,3,2],x,X,n,4);`):

elif nops([args])=1 and op(1,[args])=InfoV then
print(`InfoV(SetC,x,N): verbose form of Info(SetC,x,N) (q.v.). Try: `):
print(`InfoV({[2,3,2]},x,30);`):

elif nops([args])=1 and op(1,[args])=InitialStatesS then
print(`InitialStatesS(SetC): the set of initial states for a set of compositions SetC. Try:`):
print(`InitialStatesS({[1,2,1],[4,1]});`):

elif nops([args])=1 and op(1,[args])=InitialWt then
print(`InitialWt(SetC1,x,t): The  initial weight of the set SetC1. Try:`):
print(`InitialWt({[1,3,1],[3,1,3]},x,t);`):

elif nops([args])=1 and op(1,[args])=InitialWtSt then
print(`InitialWtSt(SetC,S,x,t): The  initial weight of the state S with set of compositions SetC. Try:`):
print(`InitialWtSt({[2,1,2]},[2,1,2],x,t);`):

elif nops([args])=1 and op(1,[args])=InitialWtStX then
print(`InitialWtStX(SetC,S,x,t,X): The  initial weight of the state S with set of compositions SetC for the GENRALIZED case. Try:`):
print(`InitialWtStX({[2,1,2]},[2,1,2],x,t,X);`):

elif nops([args])=1 and op(1,[args])=InitialWtX then
print(`InitialWtX(SetC1,x,t,X): The  initial GENERALIZED weight of the set SetC1. Try:`):
print(`InitialWtX({[1,3,1],[3,1,3]},x,t,X);`):

elif nops([args])=1 and op(1,[args])=IsCon then
print(`IsCon(C1,C2): is composition C1 contained in composition C2? Try: `):
print(` IsCon([1,2,1],[3,4,3,2]); `):

elif nops([args])=1 and op(1,[args])=IsLess1 then
print(`IsLess1(L1,L2) is the list L1 pointwise less than L2. try:`):
print(`IsLess1([3,4],[4,5]);`):

elif nops([args])=1 and op(1,[args])=IsUni then
print(`IsUni(C): Is the composition C unimodal? Try:`):
print(`IsUni([1,2,1]); IsUni([2,1,2]);`):


elif nops([args])=1 and op(1,[args])=Kesem then
print(`Kesem(R,x,n): inputs an expression of the form Polynomial in x + Sum(C_i/(1-x)^i, i=1..r)+ Sum(C_i/(1-2x)^j, i=1..s) out`):
print(`outputs the expression of the form P(x)+ Q(n)2^n describing the coefficient of x^n, as well as the first few coefficients`):
print(`Try: `):
print(`Kesem((1+x+x^5)/(1-x)^2/(1-2*x)^3, x, n);`):

elif nops([args])=1 and op(1,[args])=KesemA then
print(`KesemA(R,x,n): inputs an expression of the form Polynomial in x + Sum(C_i/(1-x)^i, i=1..r)+ Sum(C_i/(1-2x)^j, i=1..s) out`):
print(`outputs the leading expression Q(n), of Kesem(R,x,n), w/o the 2^n `):
print(`Try: `):
print(`KesemA((1+x+x^5)/(1-x)^2/(1-2*x)^3, x, n);`):


elif nops([args])=1 and op(1,[args])=MaxC then
print(`MaxC(M): given a rectangular matrix, outputs the row vector that is the max of each column. Try:`):
print(`MaxC([[2,1,2],[1,2,1]]);`):


elif nops([args])=1 and op(1,[args])=Mishkal then
print(`Mishkal(C,SetC,X): Given a composition C, a set of compositions SetC, and a variable name X`):
print(`outputs the weight arrording to SetC of C, i.e. `):
print(`mul( X[op(c)]^NuSin(c,C), c in SetC). Try:`):
print(`Mishkal([4,2,4,3],{[1,2,3],[2,3,1]},X);`):

elif nops([args])=1 and op(1,[args])=MomsK then
print(`MomsK(C,n,K): Inputs a composition C, and a variable n, as well as a positive integer K, outputs explicit`):
print(`expressions for the first K moments of the r.v. "number of times C is contained in a composition of n)"`):
print(`starting with the average. Try:`):
print(`MomsK([2,3,2],n,4);`):


elif nops([args])=1 and op(1,[args])=MomsK2 then
print(`MomsK2(C1,C2,n,K): Inputs  compositions C1 and C2 of the same length, and a variable n, as well as a positive integer K, outputs `):
print(` a K+1 by K+1 matrix (presented as a list of lists), let's call it L, such the L[i][j] is the (i-1,j-1) mixed moment of `):
print(` the random variables "Number of occurrences of C1 (C2) in a composition of length n". `):
print(` In particular, L[2][2] is the (1,1) mixed moment. Try: `):
print(` MomsK2([2,3,2],[3,2,3],n,4); `):

elif nops([args])=1 and op(1,[args])=MomsKA then
print(`MomsKA(C,n,K): The leading asymptotic of MomsK(C,n,K) (q.v.) ABOUT THE MEAN. The first entry is the expectation, the second the variance,etc. `):
print(` Try: `):
print(`MomsKA([2,3,2],n,4);`):

elif nops([args])=1 and op(1,[args])=MomsKA2 then
print(`MomsKA2(C1,C2,n,K): Inputs  compositions C1 and C2 of the same length, and a variable n, as well as a positive integer K, outputs `):
print(` a K+1 by K+1 matrix (presented as a list of lists), let's call it L, such the L[i][j] is the (i-1,j-1) mixed moment ABOUT THE MEAN `):
print(` for the random variables "Number of occurrences of C1 (C2) in a composition of length n".  Ignoring peanuts`):
print(` In particular, L[2][2] is the covariance. Try: `):
print(` MomsKA2([2,3,2],[3,2,3],n,4); `):


elif nops([args])=1 and op(1,[args])=MomsKA2st then
print(`MomsKA2stL(C1,C2,n,K): Inputs  compositions C1 and C2 of the same length, and a variable n, as well as a positive integer K, outputs `):
print(`(i): the [ave,variance] of C1`):
print(`(ii): the [ave,variance] of C2`):
print(`(iii) a K by K matrix (presented as a list of lists), let's call it L, such the L[i][j] is the `):
print(`LIMIT of the (i-1,j-1) standardized  mixed moment of as n goes to infinity.`):
print(`In particular, L[2][2] is the limiting correlation. Try:`):
print(`MomsKA2st([2,3,2],[3,2,3],n,4);`):


elif nops([args])=1 and op(1,[args])=MomsKAst then
print(`MomsKAst(C,n,K): Like MomsKA(C,n,K) (q.v.) but starting with the the third moments, gives the limits of`):
print(`the scaled moments confirming asymtotic normality. Try:`):
print(`MomsKAst([2,3,2],n,4);`):

elif nops([args])=1 and op(1,[args])=ND2 then
print(`ND2(x,y,a): the bi-variate normal distribution . Try:`):
print(`ND2(x,y,1/2)`):

elif nops([args])=1 and op(1,[args])=ND2m then
print(`ND2m(a,r,s): the (r,s) mixed moment of ND2(x,y,a). Try:`):
print(`ND2m(1/3,1,1);`):

elif nops([args])=1 and op(1,[args])=ND2mTable then
print(`ND2mTable(a,K): the moment table of ND2(x,y,a) `):


elif nops([args])=1 and op(1,[args])=NuSin then
print(`NuSin(C1,C2): How many times does the composition C1 contained in C2. Try: `):
print(`NuSin([3,2,3], [5,4,2,4,4,2]); `):

elif nops([args])=1 and op(1,[args])=Ove then
print(`Ove(S1,S2,x,t): the overlaps between sets of compositions S1 and S2 try:`):
print(`Ove({[2,3,2]},{[4,3,2]},x,t);`):

elif nops([args])=1 and op(1,[args])=Ove1 then
print(`Ove1(C,x,t): the self overlap of a composition C . Try:`):
print(`Ove1([2,3,4,3,2],x,t);`):


elif nops([args])=1 and op(1,[args])=Reduce1 then
print(`Reduce1(Cset): given a set of compositions finds a minimal subset. Try:`):
print(`Reduce1({[1,2,1],[3,3,3],[6,4,3]});`):

elif nops([args])=1 and op(1,[args])=Reduce11 then
print(`Reduce11(Cset): given a set of compositions performs one step in kicking out superflous compositions. Try:`):
print(`Reduce11({[1,2,1],[3,3,3]});`):

elif nops([args])=1 and op(1,[args])=Skyline then
print(`Skyline(S): inputs a set of compositions, outputs the smallest composition that  includes them all. Try: `):
print(` Skyline({[1,2,1],[2,1,1],[4]}); `):

elif nops([args])=1 and op(1,[args])=StatesC then
print(`StatesC(C): The set of states for the clusters of the composition C. Try:`):
print(`StatesC([2,1,2]);`):

elif nops([args])=1 and op(1,[args])=StatesCset then
print(`StatesCset(SetC): The set of states for the clusters of the set of compositions SetC. Try:`):
print(`StatesCset({[2,1,2]});`):

elif nops([args])=1 and op(1,[args])=SuperInfo1 then
print(`SuperInfo1(k,x,N): Inputs a positive integer k, variable x, and positive integer N, gathers Info({C},x,N) for all`):
print(`compositions C of k. Try:`):
print(`SuperInfo1(5,x,30);`):

elif nops([args])=1 and op(1,[args])=SuperInfo1X then
print(`SuperInfo1X(k,x,N,X,n,r): Like SuperInfo1(k,x,N) but with additional information about the generalized case and the moments. Try:`):
print(`SuperInfo1X(5,x,20,X,n,4);`):

elif nops([args])=1 and op(1,[args])=SuperInfo1XV then
print(`SuperInfo1XV(k,x,N,X,n,r): verbose verson of SuperInfo1X(k,x,N,X,n,r) (q.v.). Try:`):
print(`SuperInfo1XV(5,x,30,X,n,4);`):

elif nops([args])=1 and op(1,[args])=SuperInfo1V then
print(`SuperInfo1V(k,x,N): Verbose form of SuperInfo1(k,x,N) (q.v.). Try: `):
print(`SuperInfo1V(5,x,30);`):

elif nops([args])=1 and op(1,[args])=SuperInfo2kr then
print(`SuperInfo2kr(k,r,x,N): Inputs  positive integers k and r, variable x, and positive integer N, gathers Info(S,x,N) `):
print(`for all subsets of the set of compositions of k with r parts. Try:`):
print(`SuperInfo2kr(5,2,x,30);`):

elif nops([args])=1 and op(1,[args])=SuperInfo2krV then
print(`SuperInfo2krV(k,r,x,N): verbose form of SuperInfo2kr(r,k,x,N) (q.v.). Try: `):
print(`SuperInfo2krV(5,2,x,30);`):


elif nops([args])=1 and op(1,[args])=SuperInfoX2Vone then
print(`SuperInfoX2Vone(m1,m2,x,X1,X2,n,r):inputs  positive integers m1, m2, (m2<m1)`):
print(` continuous variable names x,X1,X2, discrete variable name n, and`):
print(`a positive integer r >=2, outputs an article about the joint generating functions of all pairs of`):
print(`of compositions of m1 into m2 parts. Try:`):
print(`SuperInfoX2Vone(5,2,x,X1,X2,n,4);`):


elif nops([args])=1 and op(1,[args])=SuperInfoX2Vtwo then
print(`SuperInfoX2Vtwo(m1a,m1b, m2,x,X1,X2,n,r):inputs  positive integers m1a, ma1bm, m2, (m2<min(m1a,m1b))`):
print(` continuous variable names x,X1,X2, discrete variable name n, and`):
print(`a positive integer r >=2, outputs an article about the joint generating functions of all pairs of`):
print(`of compositions of m1a into m2 parts and compositions of m1b into m2 parts that are not included in each other. Try:`):
print(`SuperInfoX2Vtwo(5,6,2,x,X1,X2,n,2);`):

elif nops([args])=1 and op(1,[args])=Taur then
print(`Taur(r): the composition r(r+1)...(2*r-2)(2*r-1)...(r+1)r. For example, try:`):
print(`Taur(4);`):

elif nops([args])=1 and op(1,[args])=Wt then
print(`Wt(S,x,t): the weight of a set of compositions S. Try:`):
print(`Wt({[1,2,1],[2,1,1],[4]},x,t);`):


elif nops([args])=1 and op(1,[args])=WtCd then
print(`WtCd(n,SetC,X): The weight-enumerator of the set of compositions of length n according to X[SetC], done directly.`):
print(`For checking purposes only. Try: `):
print(` WtCd(6,{[1,2,1],[2,1,2]},X); `):

else
print(`There is no ezra for`,args):
fi:
 
end:


#Comps1(n,k): the set of compositions of n into exactly k parts. Try:
#Comps1(5,2);

Comps1:=proc(n,k) local gu,mu1,mu,i:
option remember:
if k=0 then
 if n=0 then
  RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:

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

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

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

gu:={}:

for i from 1 to n do
 mu:=Comps1(n-i,k-1):
 gu:=gu union {seq([i,op(mu1)],mu1 in mu)}:
od:

gu:

end:


#Comps(n): the set of compositions of n. Try:
#Comps(5);
Comps:=proc(n) local k:
option remember:
if n=0 then
 RETURN({[]}):
fi:
{seq(op(Comps1(n,k)),k=1..n)}:
end:


#IsLess1(L1,L2) is the list L1 pointwise less than L2. try:
#IsLess1([3,4],[4,5]);
IsLess1:=proc(L1,L2) local i:
if nops(L1)<>nops(L2) then
 RETURN(FAIL):
fi:

for i from 1 to nops(L1) do
 if L1[i]>L2[i] then
  RETURN(false):
 fi:
od:
true:
end:


#IsCon(C1,C2): is composition C1 contained in composition C2? Try
#IsCon([1,2,1],[3,4,3,2]);
IsCon:=proc(C1,C2) local i:

if nops(C1)>nops(C2) then
 RETURN(false):
fi:

for i from 1 to nops(C2)-nops(C1)+1 do
if IsLess1(C1,[op(i..i+nops(C1)-1,C2)]) then
 RETURN(true):
fi:
od:

false:

end:



#Av(C,SetC): does the composition C avoid the set of compositions SetC. Try
#Av([3,4,2],{[2,5],[4,1]});
Av:=proc(C,SetC) local C1:

for C1 in SetC do
 if IsCon(C1,C) then
   RETURN(false):
 fi:
od:

true:

end:

#GoodCnaive(n,SetC): The set of compositions of n avooding the set of compositions SetC as subcompositions. Naive approach. Try:
#GoodCnaive(6,{[1,2,1]});

GoodCnaive:=proc(n,SetC) local mu,mu1,gu:
option remember:
mu:=Comps(n):
gu:={}:

 for mu1 in mu do
 if Av(mu1,SetC) then
   gu:=gu union {mu1}:
 fi:
od:
gu:
end:


#GoodC1naive(n,k,SetC): The set of compositions of n into k partsavooding the set of compositions SetC as subcompositions. Naive approach. Try:
#GoodC1naive(6,3,{[1,2,1]});

GoodC1naive:=proc(n,k,SetC) local mu,mu1,gu:
option remember:
mu:=Comps1(n,k):
gu:={}:

 for mu1 in mu do
 if Av(mu1,SetC) then
   gu:=gu union {mu1}:
 fi:
od:
gu:
end:


#GoodC1(n,k,SetC): the set of compositions of n into exactly k parts avoiding the set of compositions SetC . Try:
#GoodC1(5,2,{[1,2,1]});

GoodC1:=proc(n,k,SetC) local gu,mu1,mu,i,cand:
option remember:
if k=0 then
 if n=0 then
  RETURN({[]}):
 else
  RETURN({}):
 fi:
fi:

if  (n<k or n=k or k=1) then
 RETURN(GoodC1naive(n,k,SetC)):
fi:

gu:={}:

for i from 1 to n do
 mu:=GoodC1(n-i,k-1,SetC):

  for mu1 in mu do
   cand:=[i,op(mu1)]:
    if Av(cand,SetC)  then
     gu:=gu union {cand}:
    fi:
  od:
 
od:

gu:

end:

#GoodC(n,SetC): The set of compositions of n avooding the set of compositions SetC as subcompositions. Try:
#GoodC(6,{[1,2,1]});
GoodC:=proc(n,SetC) local k:
option remember:
{seq(op(GoodC1(n,k,SetC)),k=1..n)}:
end:

#GoodCseqD(N,SetC): the first N+1 terms, starting at n=0 and ending with n=N,
# of the enumerating sequence for compositions of n avoiding as sub-compositions the
#members of the set SetC, done DIRECTLY (by counting). For checking only! Try:
#GoodCseqD(10,{[2,3,2]});
GoodCseqD:=proc(N,SetC) local n1:
[1,seq(nops(GoodC(n1,SetC)),n1=1..N)]:
end:

#GFoneU(C,x inputs a variable name, t, and a composition C, and outputs the generating function
#of the enumerating sequence of compositions of n avooding C as a sub-composition. Try:
#GFoneU([2,3,2],x
GFoneU:=proc(C,x) local L,t:
option remember:

if C=[] then
 RETURN(0):
fi:

if not IsUni(C) then
 print(C, `is not unimodal `):
 RETURN(FAIL):
fi:

L:=normal(Wt({C},x,t) /(1+Ove1(C,x,t))):

L:=subs(t=1/(1-x),L):
normal(1/(1-x/(1-x)-L)):

end:


#Ove1(C,x,t): the self overlap of a composition C . Try
#Ove1([2,3,4,3,2],x,t);
Ove1:=proc(C,x,t) local k,j,i:
k:=nops(C):

if not (type(C,list) and type(x,symbol) and type(t,symbol) ) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

add( t^(i-1)*x^(add(C[j],j=1..i-1)+add( max(C[j]-C[j-i+1],0)  , j=i..k)),i=2..k):
end:



#Taur(r): the composition r(r+1)...(2*r-2)(2*r-1)...(r+1)r. For example, try:
#Taur(4);
Taur:=proc(r) local i,gu:
gu:=[seq(i,i=r..2*r-2)]:
[op(gu),2*r-1,seq(gu[nops(gu)+1-i],i=1..nops(gu))]:
end:



#Skyline(S): inputs a set of compositions, outputs the smallest composition that  includes them all. Try
#Skyline({[1,2,1],[2,1,1],[4]});
Skyline:=proc(S) local s,k,S1,i:

k:=max(seq(nops(s), s in S)):

S1:={seq([op(s),0$(k-nops(s))],s in S)}:
[seq(max( seq(s[i],s in S1)),i=1..k)]:
end:


#Wt(S,x,t): the weight of a set of compositions S. Try:
#Wt({[1,2,1],[2,1,1],[4]},x,t);
Wt:=proc(S,x,t) local gu:
gu:=Skyline(S):

(-1)^nops(S)*t^nops(gu)*x^convert(gu,`+`):
end:





#Ove(S1,S2,x,t): the overlaps between sets of compositions S1 and S2 try:
#Ove({[2,3,2]},{[4,3,2]},x,t);
Ove:=proc(S1,S2,x,t) local C1,C2,k,i,j,lu:
C1:=Skyline(S1):
C2:=Skyline(S2):
k:=nops(C1):

lu:=0:

for i from 2 to k do

 if k-i+1<=nops(C2) then

  lu:=lu+t^(i-1)*x^(add(C1[j],j=1..i-1)+add( max(C1[j]-C2[j-i+1],0)  , j=i..k) ):
 else
   lu:=lu+t^(i-1)*x^(add(C1[j],j=1..i-1)+add( max(C1[j]-C2[j-i+1],0)  , j=i..k-i+1)+add( C1[j]  , j=k-i+2..k))
 fi:

od:
end:





#IsUni(C): Is the composition C unimodal? Try:
#IsUni([1,2,1]); IsUni([2,1,2]);
IsUni:=proc(C) local i:
if nops(C)<=1 then
 RETURN(true):
fi:
if C[2]<C[1] then
   RETURN(false):
fi:

for i from 2 to nops(C)-1 do
 if C[i-1]>C[i] and C[i]<C[i+1] then
   RETURN(false):
 fi:
od:
true:
end:

#Clus1(C,k,m): The set of k by m clusters with one composition. Try
#Clus1([2,1,2],3,6);
Clus1:=proc(C,k,m) local gu,r,gu1,gu11,i1:
option remember:

if not (type(C,list) and type(k, integer) and type(m,integer) and k>0 and m>0) then
 RETURN({}):
fi:

if k=1 then
 if m=nops(C) then
  RETURN({[C]}):
  else
  RETURN({}):
 fi:
fi:

gu:={}:

for r from 1 to nops(C)-1 do
 gu1:=Clus1(C,k-1,m-r):
 gu1:={seq([[op(C),0$(m-nops(C))], seq([0$(r),op(gu11[i1])],i1=1..k-1)],gu11 in gu1)}:
 gu:=gu union gu1:
od:
gu:


end:

#Clus(C,k): The set of k by m clusters with one composition with ANY m. Try
#Clus([2,1,2],3);
Clus:=proc(C,k) local m,m1,gu:
option remember:

for m from 1 while Clus1(C,k,m)={} do od:
m1:=m:

gu:={}:

for m from m1 while Clus1(C,k,m)<>{} do
 gu:=gu union Clus1(C,k,m):
od:
gu:

end:


#MaxC(M): given a rectangular matrix, outputs the row vector that is the max of each column. Try
#MaxC([[2,1,2],[1,2,1]]);
MaxC:=proc(M) local i,i1:

[seq(max(seq(M[i][i1],i=1..nops(M))),i1=1..nops(M[1]))]:

end:


#Aba(C,S,r,x,t): inputs a composition C, a state S, and r between 2 and nops(C), outputs the new state and the
#change in weight in terms of x and t. Try
#Aba([2,1,2,1,2],[2,1,2,2],4,x,t);
Aba:=proc(C,S,r,x,t) local gu,mu,i:
option remember:
if (nops(C)<>nops(S)  or r<=1 or r>nops(C)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=[[op(C),0$(r-1)],[0$(r-1),op(S)]]:
mu:=[seq(max(gu[1][i],gu[2][i]),i=1..nops(C)+r-1)]:

[[op(1..nops(S),mu)],-t^(r-1)*x^(convert(mu,`+`)-convert(S,`+`))]:
end:

#Avot(C,S): all the fathers of the state S with composition C. Try:
#Avot([2,1,2,1,2],[2,1,2,2]);
Avot:=proc(C,S) local r,x,t:
{seq(Aba(C,S,r,x,t)[1],r=2..nops(C))}:
end:



#StatesC(C): The set of states for the clusters of the composition C. Try:
#StatesC([2,1,2]);
StatesC:=proc(C) local S,gu,gu1,gu2, gua:

S:=C:

gu:={S}:

gu1:=Avot(C,S) union gu:

while gu<>gu1 do
 gu2:=gu1 union {seq(op(Avot(C,gua)),gua in gu1)}:
 gu:=gu1:
 gu1:=gu2:
od:
gu1:
end:




#GFone(C,x): the generating function of the sequence enumerating compositions of n avoiding the sub-composition C. Try
#GFone([2,1,2],x);
GFone:=proc(C,x) local gu,eq,var, T,S,gu1,gu2,r,lu,t,X,eq1,L:
option remember:
gu:=StatesC(C):

S:=C:

for  gu1 in gu do
 for gu2 in gu do
  T[gu1,gu2]:=0:
 od:
od:

for gu1 in gu do
 for r from 2 to nops(C) do
  lu:=Aba(C,gu1,r,x,t):
   T[lu[1],gu1]:=T[lu[1],gu1]+lu[2]:
 od:
od:




var:={seq(X[gu1],gu1 in gu)}:

eq1:=X[S]+x^convert(C,`+`)*t^nops(C)-add(T[S,gu1]*X[gu1],gu1 in gu):
eq:={eq1}:


for gu1 in gu minus {S} do
 eq1:=X[gu1]-add(T[gu1,gu2]*X[gu2],gu2 in gu):
 eq:=eq union {eq1}:
od:



var:=solve(eq,var):

lu:=normal(add(subs(var,X[gu1]),gu1 in gu)):

L:=normal(subs(t=1/(1-x),lu)):

factor(normal(1/(1-x/(1-x)-L))):
end:

  

#CheckGFone(C,N): checks GFone(C,x) up to N terms. Try:
#CheckGFone([2,1,2],12);
CheckGFone:=proc(C,N) local gu,i,x,mu:
gu:=GFone(C,x):
gu:=taylor(gu,x=0,N+1):
gu:=[seq(coeff(gu,x,i),i=0..N)]:
mu:=GoodCseqD(N,{C}):

if gu=mu then
 RETURN(true):
else
 RETURN(mu,gu):
fi:


end:

#Reduce11(Cset): given a set of compositions performs one step in kicking out superflous compositions. Try:
#Reduce11({[1,2,1],[3,3,3]});
Reduce11:=proc(Cset) local C1,C2:

if nops(Cset)=1 then
 RETURN(Cset):
fi:

for C1 in Cset do
 for C2 in Cset do
  if C1<>C2 and IsCon(C1,C2) then
    RETURN(Cset minus {C2}):
 fi:
 od:
od:
Cset:
end:


#Reduce1(Cset): given a set of compositions performs finds a minimal subset. Try
#Reduce1({[1,2,1],[3,3,3]});
Reduce1:=proc(Cset) local gu,gu1,gu2:
gu:=Cset:
gu1:=Reduce11(gu):

while gu<>gu1 do
gu2:=Reduce11(gu1):
gu:=gu1:
gu1:=gu2:
od:
gu:
end:




#InitialStatesS(SetC): the set of initial states for a set of compositions SetC. Try
#InitialStatesS({[1,2,1],[4,1]});
InitialStatesS:=proc(SetC) local gu,gu1:

if Reduce1(SetC)<>SetC then
 print(`Not reduced`):
 RETURN(FAIL):
fi:


gu:=powerset(SetC) minus {{}}:

{seq(Skyline(gu1),gu1 in gu)}:
end:




#InitialWt(SetC1,x,t): The  initial weight of the set SetC1. Try:
#InitialWt({[1,3,1],[3,1,3]},x,t);
InitialWt:=proc(SetC1,x,t) local gu:
option remember:
 gu:=Skyline(SetC1):

 (-1)^nops(SetC1)*t^nops(gu)*x^convert(gu,`+`):
end:

#InitialWtSt(SetC,S,x,t): The  initial weight of the state S with set of compositions SetC. Try:
#InitialWtSt({[2,1,2]},[2,1,2],x,t);
InitialWtSt:=proc(SetC,S,x,t) local lu1,gu,lu:
lu:=powerset(SetC):

if not member(S,InitialStatesS(SetC)) then
  RETURN(FAIL):
fi:

gu:=0:

for lu1 in lu do
 if Skyline(lu1)=S then
  gu:=gu+InitialWt(lu1,x,t):
fi:
od:
gu:
end:

 




#AbaSet(SetC1,S,r,x,t): inputs a set of compositions SetC1, a state S, and r between 2 and nops(Skyline(SetC)), outputs the new state and the
#change in weight in terms of x and t. Try
#AbaSet({[2,1,2,1,2]},[2,1,2,2],4,x,t);
AbaSet:=proc(SetC1,S,r,x,t) local C1,gu,mu,i:
option remember:

if Reduce1(SetC1)<>SetC1 then
 print(`Not reduced`):
 RETURN(FAIL):
fi:

C1:=Skyline(SetC1):

if ( r<=1 or r>nops(C1)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=[[op(C1),0$(r-1+nops(S)-nops(C1))],[0$(r-1),op(S)]]:
mu:=[seq(max(gu[1][i],gu[2][i]),i=1..nops(S)+r-1)]:

[[op(1..nops(S),mu)],(-1)^nops(SetC1)*t^(r-1)*x^(convert(mu,`+`)-convert(S,`+`))]:
end:




#AvotSet(SetC,S): all the fathers of the state S with a set of compositions SetC. Try:
#AvotSet({[2,1,2,1,2]},[2,1,2,2]);
AvotSet:=proc(SetC,S) local SetC1,C1,r,x,t,gu,lu:

if Reduce1(SetC)<>SetC then
 print(`Not reduced`):
 RETURN(FAIL):
fi:

gu:={}:

lu:=powerset(SetC) minus {{}}:

for SetC1 in lu do
 C1:=Skyline(SetC1):

gu:=gu union {seq(AbaSet(SetC1,S,r,x,t)[1],r=2..nops(C1))}:

od:

gu:

end:

#StatesCset(SetC): The set of states for the clusters of the set of compositions SetC. Try:
#StatesCset({[2,1,2]});
StatesCset:=proc(SetC) local lu,lu1,gu,gu1,gu2, gua:

if Reduce1(SetC)<>SetC then
 print(`Not reduced`):
 RETURN(FAIL):
fi:

lu:=InitialStatesS(SetC):
gu:=InitialStatesS(SetC):
gu1:=gu union {seq(op(AvotSet(SetC,lu1)),lu1 in lu)}:

while gu<>gu1 do
 gu2:=gu1 union {seq(op(AvotSet(SetC,gua)),gua in gu1)}:
 gu:=gu1:
 gu1:=gu2:
od:
gu1:
end:



#GFset(SetC,x): inputs a set of compositions S and outputs the generating function enumerating compositions that avoid
#the members of S as sub-compositions. Try:
#GFset({[2,3,2],[4,1,4]},x);
GFset:=proc(S,x) local SetC,gu,ku,ku1,lu,gu1,gu2,T,t,eq,var,vu,L,X,r,eq1:
option remember:
if not type(S,set) then
 print(`bad input`):
 RETURN(FAIL):
fi:

if nops({seq(nops(SetC),SetC in S)})<>1 then
 print(`The members of`, S, `are not of the same length`):
 RETURN(FAIL):
fi:


SetC:=Reduce1(S):

gu:=StatesCset(SetC):

vu:=InitialStatesS(SetC):


ku:=powerset(SetC) minus {{}}:

for  gu1 in gu do
 for gu2 in gu do
  T[gu1,gu2]:=0:
 od:
od:

for gu1 in gu do
 for ku1 in ku do
 for r from 2 to nops(Skyline(ku1)) do
  lu:=AbaSet(ku1,gu1,r,x,t):
   T[lu[1],gu1]:=T[lu[1],gu1]+lu[2]:
 od:
od:
od:



var:={seq(X[gu1],gu1 in gu)}:

eq:={}:

for gu1 in gu  do


 eq1:=X[gu1]-add(T[gu1,gu2]*X[gu2],gu2 in gu):


if member(gu1,vu) then

 eq1:=eq1-InitialWtSt(SetC,gu1,x,t):

fi:

 eq:=eq union {eq1}:
od:




var:=solve(eq,var):

lu:=normal(add(subs(var,X[gu1]),gu1 in gu)):

L:=normal(subs(t=1/(1-x),lu)):

factor(normal(1/(1-x/(1-x)-L))):
end:

#CheckGFset11(SetC,N): checks GFset(SetC,x) up to N terms. Try:
#CheckGFset11({[2,1,2]},12);
CheckGFset11:=proc(SetC,N) local gu,i,x,mu:
option remember:
gu:=GFset(SetC,x):
gu:=taylor(gu,x=0,N+1):
gu:=[seq(coeff(gu,x,i),i=0..N)]:
mu:=GoodCseqD(N,SetC):

if gu=mu then
 RETURN(true):
else
print(mu,gu):
 RETURN(false):
fi:


end:


#CheckGFset1(r,k,N): checks GFset(SetC,x) up to N terms for all non-empty subsets of the set of compositions of r with k parts. Try:
#CheckGFset1(5,2,12);
CheckGFset1:=proc(r,k,N) local SetC,gu:
gu:=powerset(Comps1(r,k)) minus {{}}:
 for SetC in gu do
  if not CheckGFset11(Reduce1(SetC),N) then
   print(`Wrong for the set`, SetC):
   RETURN(false):
 fi:
 od:
true:
end:


#CheckGFset(R,N): checks GFset(SetC,x) up to N terms for all non-empty subsets of the set of compositions of <=R with the same number of parts. Try:
#CheckGFset(3,12);
CheckGFset:=proc(R,N) local r,k:

for r from 1 to R do
 for k from 1 to r do

 if not CheckGFset1(r,k,N) then
  print(r,k):
   RETURN(false):
 fi:
od:
od:

true:

end:


#Info(S,x,N): inputs a set of compositions of the same length, SetC, and a variable x, and a positive integer N, returns
#(i) The generating function of the seqeunce enumerating compositions that avoid the compositions in SetC. let's call it a(n)
#(ii) the limit of a(n+1)/a(n), as n goes to infinity, let's call it alpha
#(iii) The limit of a(n)/alpha^n as n goes to infinity
#(iv): The first N terms of the enumerating sequence. Try:
#Info({[2,3,2]},x,30);
Info:=proc(S,x,N) local C,gu,lu,i,hal,aluf,si,alpha,gu1,Sid:
option remember:
Digits:=20:
if not type(S,set) then
 print(`bad input`):
 RETURN(FAIL):
fi:

if S={} then
 RETURN([S,normal(1+x/(1-2*x)),[1,seq(2^(i-1),i=1..N)],2,1/2]):
fi:

if nops({seq(nops(C),C in S)})<>1 then
 print(`The members of`, S, `are not of the same length`):
 RETURN(FAIL):
fi:


gu:=GFset(S,x):

gu1:=taylor(gu,x=0,max(N,201)):

Sid:=[seq(coeff(gu1,x,i),i=0..N)]:

lu:=[fsolve(denom(gu))]:

if lu=[] then
 RETURN([S,gu,Sid]):
fi:



aluf:=lu[1]:
si:=abs(lu[1]):

for i from 2 to nops(lu) do
 hal:=abs(lu[i]):
 if hal<si then
   aluf:=lu[i]:
   si:=hal:
 fi:
od:

if aluf<0 then
 RETURN([S,gu,Sid]):
fi:

alpha:=1/si:



if abs(coeff(gu1,x,200)/coeff(gu1,x,199)-alpha)>10^(-12) then
 RETURN([S,gu,Sid]):
fi:


C:=-alpha/subs(x=si,diff(1/gu,x)):

if abs(coeff(gu1,x,200)/alpha^200-C)>10^(-12) then
 RETURN([S,gu,Sid,alpha]):
fi:

[S,gu,Sid,alpha,C]:
end:




#InfoV(S,x,N): Verbose form of Info(S,x,N) (q.v.). Try:
#InfoV({[2,3,2]},x,30);
InfoV:=proc(S,x,N) local gu,a,n:
gu:=Info(S,x,N):

if gu=FAIL then
 RETURN(FAIL):
fi:

print(``):
if nops(S)=1 then
print(`Theorem: Let a(n) be the number of compositions of n avoiding, as a subcomposition`, S[1]):
else
print(`Theorem: Let a(n) be the number of compositions of n avoiding as subcompositions, the members of the set`, S):
fi:
print(``):
print(`We have the following rational function for the (ordinary) generating function of a(n)`):
print(``):
print(Sum(a(n)*x^n,n=0..infinity)=gu[2]):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu[2]):
print(``):
print(`The first`, N+1, `terms of a(n),  starting at n=0, are`):
print(``):
print(gu[3]):
print(``):

if nops(gu)>=4 then
print(`The limit of a(n+1)/a(n) as n goes to infinity is`):
print(``):
lprint(evalf(gu[4],12)):
print(``):
fi:

if nops(gu)=5 then
print(`a(n) is asymptotic to`):
print(``):
lprint(evalf(gu[5],12)*evalf(gu[4],12)^n):
print(``):
fi:
print(``):
print(`----------------------------------------------------------------------------`):
print(``):
end:

#InfoVth(S,x,N,co): Verbose form of Info(S,x,N) (q.v.), with theorem numbered. Try:
#InfoVth({[2,3,2]},x,30,3);
InfoVth:=proc(S,x,N,co) local gu,a,n:
gu:=Info(S,x,N):

if gu=FAIL then
 RETURN(FAIL):
fi:

print(``):
if nops(S)=1 then
print(`Theorem Number`, co, `: Let a(n) be the number of compositions of n avoiding, as a subcomposition`, S[1]):
else
print(`Theorem Number`,co, `: Let a(n) be the number of compositions of n avoiding as subcompositions, the members of the set`, S):
fi:
print(``):
print(`We have the following rational function for the (ordinary) generating function of a(n)`):
print(``):
print(Sum(a(n)*x^n,n=0..infinity)=gu[2]):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu[2]):
print(``):
print(`The first`, N+1, `terms of a(n),  starting at n=0, are`):
print(``):
print(gu[3]):
print(``):

if nops(gu)>=4 then
print(`The limit of a(n+1)/a(n) as n goes to infinity is`):
print(``):
lprint(evalf(gu[4],12)):
print(``):
fi:

if nops(gu)=5 then
print(`a(n) is asymptotic to`):
print(``):
lprint(evalf(gu[5],12)*evalf(gu[4],12)^n):
print(``):
fi:
print(``):
print(`----------------------------------------------------------------------------`):
print(``):
end:


#SuperInfo1(k,x,N): Inputs a positive integer k, variable x, and positive integer N, gathers Info({C},x,N) for all
#compositions C of k. Try:
#SuperInfo1(5,x,30);
SuperInfo1:=proc(k,x,N) local lu,mu,mu1,T,lu1,ku,ru,K1,K2,i,gu:
lu:=Comps(k):
mu:={seq(GFone(lu[i],x),i=1..nops(lu))}:

for mu1 in mu do
 T[mu1]:={}:
od:

for lu1 in lu do
 ku:=GFset({lu1},x):
 T[ku]:=T[ku] union {lu1}:
od:

gu:={}:

for mu1 in mu do
 ru:=Info({T[mu1][1]},x,N):
 if nops(ru)<=3 then
 gu:=gu union {1}:
  K1[mu1]:=1:
  K2[1]:=mu1:
 else
 gu:=gu union {ru[4]}:   

   
 K1[mu1]:=ru[4]:
  K2[ru[4]]:=mu1:
   
 fi:
od:

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



[seq([ gu[i],  T[K2[gu[i]]], Info({T[K2[gu[i]]][1]},x,N)] ,i=1..nops(gu))]:

end:




#SuperInfo1V(k,x,N): verbose verson of SuperInfo1(k,x,N) (q.v.). Try:
#SuperInfo1V(5,x,30);
SuperInfo1V:=proc(k,x,N) local gu,co,t0:
t0:=time():
gu:=SuperInfo1(k,x,N):



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

print(`Generating functions and Growth rates for the Enumerating Sequences for Comositions Avoiding Each of the compositions of`, k):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
co:=1:
print(`The compositions of`, k, `that yield polynomial growth (hence growth ratio 1, that is the lowest) are`, op(gu[1][2])):
print(``):
InfoVth({gu[1][2][1]},x,N,co):
print(``):

for co from 2 to nops(gu) do
print(`The compositions of`, k, `that yield the`,co, `-th largest growth, that is`, gu[co][1], ` are `,  op(gu[co][2]) ):
print(``):
InfoVth({gu[co][2][1]},x,N,co):
print(``):

od:
print(``):
print(`This ends this article, that took`, time()-t0, `seconds to generate. `):
print(``):
print(`----------------------------------------------------------`):
print(``):
end:



#SuperInfo2kr(k,r,x,N): Inputs  positive integers k and r, variable x, and positive integer N, gathers Info(S,x,N) 
#for all subsets of the set of compositions of k with r parts. Try:
#SuperInfo2kr(5,2,x,30);
SuperInfo2kr:=proc(k,r,x,N) local lu,mu,mu1,T,lu1,ku,ru,K1,K2,i,gu:

lu:=Comps1(k,r):
lu:=powerset(lu) minus {{}}:

mu:={seq(GFset(lu[i],x),i=1..nops(lu))}:


for mu1 in mu do
 T[mu1]:={}:
od:

for lu1 in lu do
 ku:=GFset(lu1,x):
 T[ku]:=T[ku] union {lu1}:
od:


gu:={}:

for mu1 in mu do
 ru:=Info(T[mu1][1],x,N):
 if nops(ru)<=3 then
 gu:=gu union {1}:
  K1[mu1]:=1:
  K2[1]:=mu1:
 else
 gu:=gu union {ru[4]}:   

   
 K1[mu1]:=ru[4]:
  K2[ru[4]]:=mu1:
   
 fi:
od:

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



#[seq([gu[i],K2[gu[i]],T[K2[gu[i]]] ] ,i=1..nops(gu))]:

[seq([ gu[i],  T[K2[gu[i]]], Info(T[K2[gu[i]]][1],x,N)] ,i=1..nops(gu))]:

end:


#SuperInfo2krV(k,r,x,N): verbose verson of SuperInfo2kr(k,r,x,N) (q.v.). Try:
#SuperInfo2krV(5,2,x,30);
SuperInfo2krV:=proc(k,r,x,N) local gu,co,t0:
t0:=time():
gu:=SuperInfo2kr(k,r,x,N):



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

print(`Generating functions and Growth rates for the Enumerating Sequences for Comositions Avoiding Each of the Possible subsets of`):
print(`the set of compositions of`, k, `into exactly`, r, `parts `):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
for co from 1 to nops(gu) do
print(`The subsets of the set of compositions of`, k, `with exactly`, r, `parts, that yield the`):
print(co, `-th largest growth, that is`, gu[co][1], ` are `,  op(gu[co][2]) ):
print(``):
InfoVth(gu[co][2][1],x,N,co):
print(``):

od:
print(``):
print(`This ends this article, that took`, time()-t0, `seconds to generate. `):
print(``):
print(`----------------------------------------------------------`):
print(``):
end:



#################Start the statistical part where we keep track of the number of occurrences of each offender

#InitialWtX(SetC1,x,t,X): The  initial GENERALIZED weight of the set SetC1. Try:
#InitialWtX({[1,3,1],[3,1,3]},x,t,X);
InitialWtX:=proc(SetC1,x,t,X) local gu,C:
option remember:
 gu:=Skyline(SetC1):

mul( (X[op(C)]-1), C in SetC1)*t^nops(gu)*x^convert(gu,`+`):
end:



#AbaSetX(SetC1,S,r,x,t,X): inputs a set of compositions SetC1, a state S, and r between 2 and nops(Skyline(SetC)), outputs the new state and the
#change in weight in terms of x and t for the GENERALIZED case keeping track of offenders. Try
#AbaSetX({[2,1,2,1,2]},[2,1,2,2],4,x,t,X);
AbaSetX:=proc(SetC1,S,r,x,t,X) local C1,gu,mu,i,C:
option remember:

if Reduce1(SetC1)<>SetC1 then
 print(`Not reduced`):
 RETURN(FAIL):
fi:

C1:=Skyline(SetC1):

if ( r<=1 or r>nops(C1)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=[[op(C1),0$(r-1+nops(S)-nops(C1))],[0$(r-1),op(S)]]:
mu:=[seq(max(gu[1][i],gu[2][i]),i=1..nops(S)+r-1)]:

[[op(1..nops(S),mu)],mul( X[op(C)]-1, C in SetC1)*t^(r-1)*x^(convert(mu,`+`)-convert(S,`+`))]:
end:



#InitialWtStX(SetC,S,x,t,X): The  initial weight of the state S with set of compositions SetC for the GENRALIZED case. Try:
#InitialWtStX({[2,1,2]},[2,1,2],x,t,X);
InitialWtStX:=proc(SetC,S,x,t,X) local lu1,gu,lu:
lu:=powerset(SetC):

if not member(S,InitialStatesS(SetC)) then
  RETURN(FAIL):
fi:

gu:=0:

for lu1 in lu do
 if Skyline(lu1)=S then
  gu:=gu+InitialWtX(lu1,x,t,X):
fi:
od:
gu:
end:



#GFsetX(SetC,x,X): inputs a set of compositions SetC, abd variable names x, and T,
# and outputs the generating function enumerating compositions accodring to the weight
# Product(X[C]^(Number of occurrences of C as a subcomposition over all C in SetC
#Try:
#GFsetX({[2,3,2],[4,1,4]},x,X);
GFsetX:=proc(S,x,X) local SetC,gu,ku,ku1,lu,gu1,gu2,T,t,eq,var,vu,L,X1,r,eq1:
option remember:
if not type(S,set) then
 print(`bad input`):
 RETURN(FAIL):
fi:

if nops({seq(nops(SetC),SetC in S)})<>1 then
 print(`The members of`, S, `are not of the same length`):
 RETURN(FAIL):
fi:


SetC:=Reduce1(S):

gu:=StatesCset(SetC):

vu:=InitialStatesS(SetC):


ku:=powerset(SetC) minus {{}}:

for  gu1 in gu do
 for gu2 in gu do
  T[gu1,gu2]:=0:
 od:
od:

for gu1 in gu do
 for ku1 in ku do
 for r from 2 to nops(Skyline(ku1)) do
  lu:=AbaSetX(ku1,gu1,r,x,t,X):
   T[lu[1],gu1]:=T[lu[1],gu1]+lu[2]:
 od:
od:
od:



var:={seq(X1[gu1],gu1 in gu)}:

eq:={}:

for gu1 in gu  do


 eq1:=X1[gu1]-add(T[gu1,gu2]*X1[gu2],gu2 in gu):


if member(gu1,vu) then

 eq1:=eq1-InitialWtStX(SetC,gu1,x,t,X):

fi:

 eq:=eq union {eq1}:
od:




var:=solve(eq,var):

lu:=normal(add(subs(var,X1[gu1]),gu1 in gu)):

L:=normal(subs(t=1/(1-x),lu)):

factor(normal(1/(1-x/(1-x)-L))):

end:



#NuSin(C1,C2): How many times does the composition C1 contained in C2. Try
#NuSin([3,2,3],[5,4,2,4,4,2]);
NuSin:=proc(C1,C2) local i,gu:

if nops(C1)>nops(C2) then
 RETURN(0):
fi:

gu:=0:
for i from 1 to nops(C2)-nops(C1)+1 do
if IsLess1(C1,[op(i..i+nops(C1)-1,C2)]) then
 gu:=gu+1:
fi:
od:

gu:

end:


#Mishkal(C,SetC,X): Given a composition C, a set of compositions SetC, and a variable name X
#outputs the weight arrording to SetC of C, i.e. 
#mul( X[op(c)]^NuSin(c,C), c in SetC). Try:
#Mishkal([4,2,4,3],{[1,2,3],[2,3,1]},X);
Mishkal:=proc(C,SetC,X) local c:
mul(X[op(c)]^NuSin(c,C),c in SetC):
end:



#WtCd(n,SetC,X): The weight-enumerator of the set of compositions of length n according to X[SetC], done directly
#For checking purposes only. Try:
#WtCd(6,{[1,2,1],[2,1,2]},X);
WtCd:=proc(n,SetC,X) local gu,gu1:
option remember:

gu:=Comps(n):

add(Mishkal(gu1,SetC,X),gu1 in gu):

end:

#CheckGFsetX11(SetC,N): checks procedure GFsetX(SetC,x,X) up to the coefficient of x^N. Try:
#CheckGFsetX11({[2,3,2],[3,2,3]},10);
CheckGFsetX11:=proc(SetC,N) local gu,x,X,i:
gu:=GFsetX(SetC,x,X):
gu:=taylor(gu,x=0,N+1):
evalb({seq(expand(coeff(gu,x,i))-WtCd(i,SetC,X),i=0..N)}={0}):
end:



#CheckGFsetX1(r,k,N): checks GFsetX(SetC,x,X) up to N terms for all non-empty subsets of the set of compositions of r with k parts. Try:
#CheckGFsetX1(5,2,12);
CheckGFsetX1:=proc(r,k,N) local SetC,gu:
gu:=powerset(Comps1(r,k)) minus {{}}:
 for SetC in gu do
  if not CheckGFsetX11(Reduce1(SetC),N) then
   print(`Wrong for the set`, SetC):
   RETURN(false):
 fi:
 od:
true:
end:


#CheckGFsetX(R,N): checks GFsetX(SetC,x) up to N terms for all non-empty subsets of the set of compositions of <=R with the same number of parts. 
#Warning: very slow
#Try:
#CheckGFsetX(3,12);
CheckGFsetX:=proc(R,N) local r,k:

for r from 1 to R do
 for k from 1 to r do

 if not CheckGFsetX1(r,k,N) then
  print(r,k):
   RETURN(false):
 fi:
od:
od:

true:

end:


#Kesem(R,x,n): inputs an expression of the form Polynomial in x + Sum(C_i/(1-x)^i, i=1..r)+ Sum(C_i/(1-2x)^j, i=1..s) out
#outputs the expression of the form P(x)+ Q(n)2^n describing the coefficient of x^n, as well as the first few coefficients
#Try
#Kesem((1+x+x^5)/(1-x)^2/(1-2*x)^3,x, n);
Kesem:=proc(R,x,n) local R1,R11,T,kv,P,Q,d,co,lu,i,Gadol,gu,mu,ku:
R1:=convert(R,parfrac):

if not type(R1,`+`) then
 RETURN(FAIL):
fi:

kv:={}:
P:=0:
Q:=0:

for i from 1 to nops(R1) do
 R11:=op(i,R1):


 if degree(denom(R11),x)=0 then
   d:=degree(R11,x):
    co:= normal(R11/x^d):
    if not type(co,numeric) then
      RETURN(FAIL):
    else
    T[d]:=co:
     kv:=kv union {d}:
    fi:


 else
  d:=degree(denom(R11),x):
    if degree(denom(normal(R11*(1-x)^d)),x)=0 then
      lu:=normal(R11*(1-x)^d):
      if not type(lu,numeric) then
         RETURN(FAIL):
      fi:
    P:=P+lu*expand(binomial(d-1+n,d-1)):       

    elif degree(denom(normal(R11*(1-2*x)^d)),x)=0 then

            lu:=normal(R11*(1-2*x)^d):
      if not type(lu,numeric) then
         RETURN(FAIL):
      fi:
    Q:=Q+lu*expand(binomial(d-1+n,d-1))*2^n:       

    else
    RETURN(FAIL):
  fi:

 fi:
od:

P:=factor(P):
Q:=factor(Q):

if kv<>{} then
Gadol:=max(op(kv)):
else
 Gadol:=0:
fi:

for i from 0 to Gadol do
 if not member(i,kv) then
   T[i]:=0:
 fi:
od:

gu:=P+Q:

mu:=[seq(T[i]+subs(n=i,gu),i=0..Gadol)]:


ku:=taylor(R,x=0,101):

if {seq(mu[i+1]-coeff(ku,x,i),i=0..Gadol),seq(subs(n=i,gu)-coeff(ku,x,i),i=Gadol+1..100)}<>{0} then
  print(mu,gu, `did not work out`):
  RETURN(FAIL):
fi:


[mu,gu]:

end:

#KesemA(R,x,n): inputs an expression of the form Polynomial in x + Sum(C_i/(1-x)^i, i=1..r)+ Sum(C_i/(1-2x)^j, i=1..s) out
#outputs the expression Q(n)2^n describing the leading term of the coefficient of x^n, as well as the first few coefficients
#Try
#KesemA((1+x+x^5)/(1-x)^2/(1-2*x)^3,x, n);
KesemA:=proc(R,x,n) local R1,R11,T,kv,P,Q,d,co,lu,i,Gadol,gu,mu,ku,Q1:
R1:=convert(R,parfrac):

if not type(R1,`+`) then
 RETURN(FAIL):
fi:

kv:={}:
P:=0:
Q:=0:
Q1:=0:
for i from 1 to nops(R1) do
 R11:=op(i,R1):


 if degree(denom(R11),x)=0 then
   d:=degree(R11,x):
    co:= normal(R11/x^d):
    if not type(co,numeric) then
      RETURN(FAIL):
    else
    T[d]:=co:
     kv:=kv union {d}:
    fi:


 else
  d:=degree(denom(R11),x):
    if degree(denom(normal(R11*(1-x)^d)),x)=0 then
      lu:=normal(R11*(1-x)^d):
      if not type(lu,numeric) then
         RETURN(FAIL):
      fi:
    P:=P+lu*expand(binomial(d-1+n,d-1)):       

    elif degree(denom(normal(R11*(1-2*x)^d)),x)=0 then

            lu:=normal(R11*(1-2*x)^d):
      if not type(lu,numeric) then
         RETURN(FAIL):
      fi:
    Q:=Q+lu*expand(binomial(d-1+n,d-1))*2^n:       
    Q1:=Q1+lu*expand(binomial(d-1+n,d-1)):       

    else
    RETURN(FAIL):
  fi:

 fi:
od:

P:=factor(P):
Q:=factor(Q):

if kv<>{} then
Gadol:=max(op(kv)):
else
 Gadol:=0:
fi:

for i from 0 to Gadol do
 if not member(i,kv) then
   T[i]:=0:
 fi:
od:

gu:=P+Q:

mu:=[seq(T[i]+subs(n=i,gu),i=0..Gadol)]:


ku:=taylor(R,x=0,101):

if {seq(mu[i+1]-coeff(ku,x,i),i=0..Gadol),seq(subs(n=i,gu)-coeff(ku,x,i),i=Gadol+1..100)}<>{0} then
#  print(mu,gu, `did not work out`):
  RETURN(FAIL):
fi:


Q1:

end:




#MomsK(C,n,K): Inputs a composition C, and a variable n, as well as a positive integer K, outputs explicit
#expressions for the first K moments of the r.v. "number of times C is contained in a composition of n)"
#starting with the average. Try:
#MomsK([2,3,2],n,4);
MomsK:=proc(C,n,K) local gu,mu,lu,R,x,X,i:
option remember:
R:=GFsetX({C},x,X):

lu:=X[op(C)]*diff(R,X[op(C)]):

mu:=Kesem(normal(subs(X[op(C)]=1,lu)),x,n):
if mu=FAIL then
  RETURN(FAIL):
else
 mu:=mu[2]/2^(n-1):
gu:=[mu]:
fi:

for i from 2 to K do
lu:=X[op(C)]*diff(lu,X[op(C)]):
mu:=Kesem(normal(subs(X[op(C)]=1,lu)),x,n):

if mu=FAIL then
  RETURN(FAIL):
else
 mu:=mu[2]/2^(n-1):
gu:=[op(gu),mu]:
fi:
od:

simplify(gu):

end:


#MomsKA(C,n,K): Inputs a composition C, and a variable n, as well as a positive integer K, outputs explicit
#expressions for the Asympotics of the first K moments ABOUT THE MEAN of the r.v. "number of times C is contained in a composition of n)"
#starting with the average. Try:
#MomsKA([2,3,2],n,4);
MomsKA:=proc(C,n,K) local gu,mu,lu,R,x,X,i,ave,r:
option remember:
R:=GFsetX({C},x,X):

lu:=X[op(C)]*diff(R,X[op(C)]):

mu:=2*KesemA(normal(subs(X[op(C)]=1,lu)),x,n):
if mu=FAIL then
  RETURN(FAIL):
else
gu:=[mu]:
fi:

for i from 2 to K do
lu:=X[op(C)]*diff(lu,X[op(C)]):
mu:=2*KesemA(normal(subs(X[op(C)]=1,lu)),x,n):

if mu=FAIL then
  RETURN(FAIL):
else
gu:=[op(gu),mu]:
fi:
od:

ave:=gu[1]:

factor([gu[1],seq((-ave)^r+add(gu[r-i]*(-ave)^(i)*binomial(r,i),i=0..r-1),r=2..K)]):

end:


#MomsKAnaked(C,n,K): Inputs a composition C, and a variable n, as well as a positive integer K, outputs explicit
#expressions for the Asympotics of the first K moments ABOUT THE MEAN of the r.v. "number of times C is contained in a composition of n)"
#starting with the average. Try:
#MomsKAnaked([2,3,2],n,4);
MomsKAnaked:=proc(C,n,K) local gu,mu,lu,R,x,X,i,r:
option remember:
R:=GFsetX({C},x,X):

lu:=X[op(C)]*diff(R,X[op(C)]):

mu:=2*KesemA(normal(subs(X[op(C)]=1,lu)),x,n):
if mu=FAIL then
  RETURN(FAIL):
else
gu:=[mu]:
fi:

for i from 2 to K do
lu:=X[op(C)]*diff(lu,X[op(C)]):
mu:=2*KesemA(normal(subs(X[op(C)]=1,lu)),x,n):

if mu=FAIL then
  RETURN(FAIL):
else
gu:=[op(gu),mu]:
fi:
od:

gu:



end:


#MomsKAst(C,n,K): Like MomsKA(C,n,K) (q.v.) but starting with the the third moments, gives the limits of
#the scaled moments confirming asymtotic normality. Try:
#MomsKAst([2,3,2],n,4);
MomsKAst:=proc(C,n,K) local gu,mu,i:

if not (type(C,list) and type(n,symbol) and type(K, integer) and K>=2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

gu:=MomsKA(C,n,K):


mu:=[gu[1],gu[2]]:

for i from 3 to K do

 if i mod 2=1 then
  if degree(gu[i],n)>=i/2 then
   print(`The `, i, `-th moment is`, gu[i] ,` and hence not asymptotically normal `):
   RETURN(mu):
     else
   mu:=[op(mu),0]:
  fi:
else
  if degree(gu[i],n)> i/2 then
  print(`The `, i, `-th moment is`, gu[i] ,` and hence not asymptotically normal `):
  RETURN(mu):
  else
  mu:=[op(mu), coeff(gu[i],n,i/2)/coeff(gu[2],n,1)^(i/2)]:
 fi:
fi:

od:

mu:

end:



#MomsK2(C1,C2,n,K): Inputs  compositions C1 and C2 of the same length, and a variable n, as well as a positive integer K, outputs 
#a K+1 by K+1 matrix (presented as a list of lists), let's call it L, such the L[i][j] is the (i-1,j-1) mixed moment of
#the random variables "Number of occurrences of C1 (C2) in a composition of length n".
#In particular, L[2][2] is the (1,1)-moment. Try:
#MomsK2([2,3,2],[3,2,3],n,4);
MomsK2:=proc(C1,C2,n,K) local gu1,lu,R,x,X,i,T,R1,R2,j:

if not (type(C1,list) and type(C2,list) and nops(C1)=nops(C2)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

T[0,0]:=1:

gu1:=MomsK(C1,n,K):

for i from 1 to K do
 T[i,0]:=gu1[i]:
od:



gu1:=MomsK(C2,n,K):

for i from 1 to K do
 T[0,i]:=gu1[i]:
od:


R:=GFsetX({C1,C2},x,X):

R1:=R:
for i from 1 to K do
 R1:=X[op(C1)]*diff(R1,X[op(C1)]):

 R2:=R1:
 for j from 1 to K do
  R2:=X[op(C2)]*diff(R2,X[op(C2)]):
 lu:=Kesem(subs({X[op(C1)]=1, X[op(C2)]=1},R2),x,n):
   
  T[i,j]:=[lu[1],simplify(lu[2]/2^(n-1)) ]:
 od:

od:

[seq( [seq(T[i,j],j=0..K)],i=0..K)]:

end:




#MomsKA2(C1,C2,n,K): Inputs  compositions C1 and C2 of the same length, and a variable n, as well as a positive integer K, outputs 
#a K+1 by K+1 matrix (presented as a list of lists), let's call it L, such the L[i][j] is the (i-1,j-1) mixed moment of
#ABOUT the mean only retaining  the 2^n part
#In particular, L[2][2] is the covariance (ignoring peanuts). Try:
#MomsKA2([2,3,2],[3,2,3],n,4);
MomsKA2:=proc(C1,C2,n,K) local gu1,lu,R,x,X,i,T,R1,R2,j,T1,r,s,ave1,ave2:

if not (type(C1,list) and type(C2,list) and nops(C1)=nops(C2)) then
 print(`Bad input`):
 RETURN(FAIL):
fi:

T[0,0]:=1:

gu1:=MomsKAnaked(C1,n,K):

for i from 1 to K do
 T[i,0]:=gu1[i]:
od:



gu1:=MomsKAnaked(C2,n,K):

for i from 1 to K do
 T[0,i]:=gu1[i]:
od:


R:=GFsetX({C1,C2},x,X):

R1:=R:
for i from 1 to K do
 R1:=X[op(C1)]*diff(R1,X[op(C1)]):

 R2:=R1:
 for j from 1 to K do
  R2:=X[op(C2)]*diff(R2,X[op(C2)]):
 lu:=KesemA(subs({X[op(C1)]=1, X[op(C2)]=1},R2),x,n):
   
 if lu=FAIL then
   RETURN(FAIL):
fi:
  T[i,j]:=2*lu:
 od:

od:

ave1:=T[1,0]:
ave2:=T[0,1]:

for r from 0 to K do
for s from 0 to K do
if (r=0 and s=0)  then
 T1[r,s]:=1:
elif r=1 and s=0 then
 T1[r,s]:=ave1:
elif r=0 and s=1 then
 T1[r,s]:=ave2:

else

lu:=add(add(binomial(r,i)*binomial(s,j)*(-ave1)^(r-i)*(-ave2)^(s-j)*T[i,j],j=0..s),i=0..r):
lu:=factor(lu):
T1[r,s]:=lu:
fi:

od:
od:

[seq( [seq(T1[i,j],j=0..K)],i=0..K)]:


end:




#MomsKA2stL(C1,C2,n,K): Inputs  compositions C1 and C2 of the same length, and a variable n, as well as a positive integer K, outputs 
#(i): the [ave,variance] of C1
#(ii): the [ave,variance] of C2
#(iii) a K by K matrix (presented as a list of lists), let's call it L, such the L[i][j] is the 
#LIMIT of the (i-1,j-1) standardized mixed moment of as n goes to infinity.
#In particular, L[2][2] is the limiting correlation. Try:
#MomsKA2st([2,3,2],[3,2,3],n,4);
MomsKA2st:=proc(C1,C2,n,K) local gu, ave1,ave2,var1,var2,T,var1a,var2a,i,j,lu:
gu:=MomsKA2(C1,C2,n,K):

if gu=FAIL then
 RETURN(FAIL):
fi:

ave1:=gu[2][1]:
var1:=gu[3][1]:

ave2:=gu[1][2]:
var2:=gu[1][3]:

var1a:=coeff(var1,n,1):
var2a:=coeff(var2,n,1):

for i from 1 to K do
 for j from 1 to K do
   lu:=gu[i+1][j+1]:
   if degree(lu,n)<(i+j)/2 then
     T[i,j]:=0:
    elif degree(lu,n)=(i+j)/2 then
        T[i,j]:=coeff(lu,n,(i+j)/2)/var1a^(i/2)/var2a^(j/2):
   else
    RETURN(FAIL):
   fi:
 od:
od:

[[ave1,var1],[ave2,var2],simplify([seq([seq(T[i,j],j=1..K)],i=1..K)])]:

end:





#ND2(x,y,a): the bi-variate normal distribution . Try:
#ND2(x,y,1/2)
ND2:=proc(x,y,a)
sqrt(1-a^2)*exp(-x^2/2-y^2/2+a*x*y)/(2*Pi)
end:


#ND2m(a,r,s): the (r,s) mixed moment of ND2(x,y,a). Try:
#ND2m(1/3,1,1);
ND2m:=proc(a,r,s) local x,y:

int(int(x^r*y^s*ND2(x,y,a),x=-infinity..infinity),y=-infinity..infinity):

end:


#ND2mTable(a,K): the moment table of ND2(x,y,a) 
ND2mTable:=proc(a,K) local i,j,lu:
if evalf(a)>=1 or evalf(a)<=0 then
 RETURN(FAIL):
fi:

lu:=ND2m(a,2,0):
[seq([seq(ND2m(a,i,j)/lu^((i+j)/2),j=1..K)],i=1..K)]:
end:




#Info1X(C,x,X,n,r): inputs a  composition C, and a variable x, and a variable name X, a variable name n, and a positive
#integer r. Outputs
#(i) The generating function ,  a rational function of x and X[op(C)] whose coefficient of x^n is
#the weight enumerator of the set of compositions of n according to the weight X[op(C)]^NumberOfTimesCisContainedInIt
#(ii) the  list consisting of the asympotic average, variance, and higher moments (about the mean) modulo peanuts.
#Try:
#Info1X([2,3,2],x,X,n,4);
Info1X:=proc(C,x,X,n,r) local gu,mu:
option remember:

gu:=GFsetX({C},x,X):

mu:=MomsKA(C,n,r):

[gu,mu]:
end:




#InfoX(S,x,X,N): an extension of Info(S,x,N) with the additional item of the complete weight enumerator given using X[op(C)]
#to indicate number of occurrences of the offending pattern. In other words it outputs
#inputs a set of compositions of the same length, SetC, and a variable x, and a positive integer N, returns
#(i) The generating function of the seqeunce enumerating compositions that avoid the compositions in SetC. let's call it a(n)
#(ii) the limit of a(n+1)/a(n), as n goes to infinity, let's call it alpha
#(iii) The limit of a(n)/alpha^n as n goes to infinity
#(iv): The first N terms of the enumerating sequence. 
#(v): The full generating function in X[op(C)]
#Try:
#InfoX({[2,3,2]},x,X,30);
InfoX:=proc(S,x,X,N) local C,gu,lu,i,hal,aluf,si,alpha,gu1,Sid,fu:
option remember:
Digits:=20:
if not type(S,set) then
 print(`bad input`):
 RETURN(FAIL):
fi:

if S={} then
 RETURN([S,normal(1+x/(1-2*x)),[1,seq(2^(i-1),i=1..N)],2,1/2]):
fi:

if nops({seq(nops(C),C in S)})<>1 then
 print(`The members of`, S, `are not of the same length`):
 RETURN(FAIL):
fi:


gu:=GFset(S,x):
fu:=GFsetX(S,x,X):

gu1:=taylor(gu,x=0,max(N,201)):

Sid:=[seq(coeff(gu1,x,i),i=0..N)]:

lu:=[fsolve(denom(gu))]:

if lu=[] then
 RETURN([S,gu,fu,Sid]):
fi:



aluf:=lu[1]:
si:=abs(lu[1]):

for i from 2 to nops(lu) do
 hal:=abs(lu[i]):
 if hal<si then
   aluf:=lu[i]:
   si:=hal:
 fi:
od:

if aluf<0 then
 RETURN([S,gu,fu,Sid]):
fi:

alpha:=1/si:



if abs(coeff(gu1,x,200)/coeff(gu1,x,199)-alpha)>10^(-12) then
 RETURN([S,gu,fu,Sid]):
fi:


C:=-alpha/subs(x=si,diff(1/gu,x)):

if abs(coeff(gu1,x,200)/alpha^200-C)>10^(-12) then
 RETURN([S,gu,fu,Sid,alpha]):
fi:

[S,gu,fu,Sid,alpha,C]:
end:


#InfoXV(S,x,N): Verbose form of InfoX(S,x,X,N) (q.v.). Try:
#InfoXV({[2,3,2]},x,30);
InfoXV:=proc(S,x,X,N) local gu,a,n,i:
gu:=InfoX(S,x,X,N):

if gu=FAIL then
 RETURN(FAIL):
fi:

print(``):
if nops(S)=1 then
print(`Theorem: Let a(n) be the number of compositions of n avoiding, as a subcomposition`, S[1]):
else
print(`Theorem: Let a(n) be the number of compositions of n avoiding as subcompositions, the members of the set`, S):
fi:
print(``):
print(`We have the following rational function for the (ordinary) generating function of a(n)`):
print(``):
print(Sum(a(n)*x^n,n=0..infinity)=gu[2]):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu[2]):
print(``):
print(`The first`, N+1, `terms of a(n),  starting at n=0, are`):
print(``):
print(gu[4]):
print(``):

print(`More generally, if we keep track of the number of occurrences, using the indexed variables`):
print(seq(X[op(S[i])],i=1..nops(S))):
print(` this more general generating function is`):
print(``):
print(gu[3]):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu[3]):
print(``):

if normal(subs({seq(X[op(S[i])]=0,i=1..nops(S))},gu[3])-gu[2])<>0 then
 print(`Something is terribly wrong`):
  RETURN(FAIL):
else
print(`Note that when we set all the variables`, seq(X[op(S[i])],i=1..nops(S)), ` to zero, we get the former generating function in`, x):
fi:

if normal(subs({seq(X[op(S[i])]=1,i=1..nops(S))},gu[3])-(1-x)/(1-2*x))<>0 then
 print(`Something is terribly wrong`):
  RETURN(FAIL):
else
print(`Also note that when we set all the variables`, seq(X[op(S[i])],i=1..nops(S)), ` to 1, we get the generating function for compositions`, (1-x)/(1-2*x)):
fi:




if nops(gu)>=5 then
print(`The limit of a(n+1)/a(n) as n goes to infinity is`):
print(``):
lprint(evalf(gu[5],12)):
print(``):
fi:

if nops(gu)=6 then
print(`a(n) is asymptotic to`):
print(``):
lprint(evalf(gu[6],12)*evalf(gu[5],12)^n):
print(``):
fi:
print(``):
print(`----------------------------------------------------------------------------`):
print(``):
end:



#SuperInfo1X(k,x,N,X,n,r): Like SuperInfo1(k,x,N) but with additional information about the generalized case and the moments. Try:
#SuperInfo1X(5,x,20,n,4);
SuperInfo1X:=proc(k,x,N,X,n,r)  local gu,i:
gu:=SuperInfo1(k,x,N):

[seq([op(gu[i]), Info1X(gu[i][3][1][1],x,X,n,r)],i=1..nops(gu))]:

end:





#SuperInfo1XV(k,x,N,X,n,r): verbose verson of SuperInfo1X(k,x,N,X,n,r) (q.v.). Try:
#SuperInfo1XV(5,x,30,X,n,4);
SuperInfo1XV:=proc(k,x,N,X,n,r) local gu,co,t0,lu,r1:
t0:=time():
gu:=SuperInfo1X(k,x,N,X,n,r):



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

print(`Generating functions and Growth rates for Weight-Enumerating Sequences According to the Number of occurrences of compositions of`, k):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
co:=1:
print(`The compositions of`, k, `that yield polynomial growth (hence growth ratio 1, that is the lowest) are`, op(gu[1][2])):
print(``):
InfoVth({gu[1][2][1]},x,N,co):
print(``):

lu:=gu[1][4]:

for co from 2 to nops(gu) do
print(`The compositions of`, k, `that yield the`,co, `-th largest growth, that is`, gu[co][1], ` are `,  op(gu[co][2]) ):
print(``):
InfoVth({gu[co][2][1]},x,N,co):
print(``):
lu:=gu[co][4]:
#print(`Also lu is`, lu):
print(``):
print(`The generating function keeping track of the number of occurrences of`, gu[co][2][1], `denoted by the variable`, X[op(gu[co][2][1])], `is `):
print(``):
print(lu[1]):
print(``):
print(`and in Maple format`):
print(``):
lprint(lu[1]):
print(``):
print(`Furthermore, the expectation of the random variable number of occurences (by containment) of `,  gu[co][2][1], ` equals `, lu[2][1] ):
print(``):
print(`The variance  equals `, lu[2][2] ):
for r1 from 3 to r do
 print(`The `, r1, `-th moment about the mean is `, lu[2][r1] ):
od:

od:
print(``):
print(`This ends this article, that took`, time()-t0, `seconds to generate. `):
print(``):
print(`----------------------------------------------------------`):
print(``):
end:




#AsyGF(f,x): inputs a rational function f and finds the pair [alpha,A] such that the coefficient of x^n in the Taylor expansion of f
#is A*alpha^n. Try
#AsyGF(1/(1-x-x^2),x);
AsyGF:=proc(f,x) local gu1,lu,aluf,si,C,alpha,hal,i:
Digits:=20:
gu1:=taylor(f,x=0,201):

lu:=[fsolve(denom(f))]:

if lu=[] then
 RETURN(FAIL):
fi:



aluf:=lu[1]:
si:=abs(lu[1]):

for i from 2 to nops(lu) do
 hal:=abs(lu[i]):
 if hal<si then
   aluf:=lu[i]:
   si:=hal:
 fi:
od:

if aluf<0 then
 RETURN(FAIL):
fi:

alpha:=1/si:



if abs(coeff(gu1,x,200)/coeff(gu1,x,199)-alpha)>10^(-12) then
 RETURN(FAIL):
fi:


C:=-alpha/subs(x=si,diff(1/f,x)):

if abs(coeff(gu1,x,200)/alpha^200-C)>10^(-12) then
 RETURN(FAIL):
fi:

[alpha,C]:
end:

#InfoX2V(C1,C2,x,X1,X2,n,r): inputs two compositions C1 and C2 none a sub-composition of the other, and outputs
#an article about 
#(i) the generating function, R0(x) enumerating compositions of the sequence enumerating compositions of
#n, AVOIDING both C1 and C2 (gotten from GFset({C1,C2},x))
#(ii) the asymptoic growth rate of that enumerating sequence enumerating 
#(iii) The Full generating function, the rational function in three variables R(x,X1,X2) where the coefficient of x^n*X1^a1*x2^a2
#is the number of compositions of n with exactly a1 occurrences (as containments) of C1 and a2 occurrences (as containments) of C2
#The expectation and  variance, in terms of n, of the random variables "number of occurrences" of C1 and C2 respectively
#defined on the set of compositions of n
#(iv) The correlation of these two random variables
#(v) A confirmation that they are indeed joint asymptotically normal with that correlation. For example, try:
#InfoX2V([2,3,2],[4,2,1],x,X,Y,n,6);
InfoX2V:=proc(C1,C2,x,X1,X2,n,r) local X,gu,a,A,c,d,lu,mu,i,j:
print(``):
if nops(Reduce1({C1,C2}))<>2 then
 print(`one of them is contained in the other`):
 RETURN(FAIL):
fi:

gu:=GFset({C1,C2},x):

print(``):
print(`Theorem: The following statements are true`):
print(``):
print(`Let a(n) be the number of compositions of n that contain neither the composition`, C1, `nor the composition`, C2):
print(``):
print(` Then `):
print(``):
print(Sum(a(n)*x^n,n=0..infinity)=gu):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu):
print(``):

lu:=AsyGF(gu,x):
if lu<>FAIL then
 print(`The asymptotic  expression for a(n) is`, lu[2]*lu[1]^n ):
fi:
print(``):
gu:=GFsetX({C1,C2},x,X):
gu:=subs({X[op(C1)]=X1, X[op(C2)]=X2},gu):
print(``):
print(` Let `, A(n,c,d), ` be the number of compositions of n with c occurrences (as containment) of the composition `, C1):
print(`and d occurrences (as containment) of the composition`, C2, ` (Note that A(n,0,0)=a(n)), Then `):
print(``):
print(Sum(Sum(Sum(A(n,c,d)*x^n*X1^c*X2^d,c=0..infinity),d=0..infinity),n=0..infinity) =gu):
print(``):
print(`and in Maple format`):
print(``):
print(``):
lprint(gu):
print(``):

gu:=MomsKA2st(C1,C2,n,r):
if gu<>FAIL then
print(``):
print(`Furthermore the average and variance of the random variable: Number of occurrences of`, C1, `are `):
print(``):
print(gu[1][1] , `and `, gu[1][2] , `respectively, while `):
print(``):
print(`Furthermore the average and variance of the random variable: Number of occurrences of`, C2, `are `):
print(``):
print(gu[2][1] , `and `, gu[2][2] , `respectively, while  the asymptotic correlation between these two random variables is`):
print(``):
print(gu[3][1][1]):
print(``):
print(`and in floating point`):
print(``):
print(evalf(gu[3][1][1])):
print(``):

mu:=ND2mTable(gu[3][1][1],r):

if mu<>FAIL then
if {seq(seq(simplify(mu[i][j]-gu[3][i][j]),j=1..r),i=1..r)}<>{0} then
 print(`Something amazing happened, the two random variables are not asymptotically normal`):
 RETURN(gu[3],mu):
else
print(`The asymptotic standardized mixed moments are indeed as they are for bi-variate normal pair with correlation`, mu[1][1]):
print(``):
print(`i.e. `, mu):
fi:

fi:
fi:

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


end:






#InfoX2Vth(C1,C2,x,X1,X2,n,r,Num): Like InfoX2Vth(C1,C2,x,X1,X2,n,r) but with the theorem Numbered. Try:
#InfoX2Vth([2,3,2],[4,2,1],x,X,Y,n,4,11);
InfoX2Vth:=proc(C1,C2,x,X1,X2,n,r,Num) local X,gu,a,A,c,d,lu,mu,i,j:
print(``):
if nops(Reduce1({C1,C2}))<>2 then
 print(`one of them is contained in the other`):
 RETURN(FAIL):
fi:

gu:=GFset({C1,C2},x):

print(``):
print(`Theorem Number`, Num, `: The following statements are true`):
print(``):
print(`Let a(n) be the number of compositions of n that contain neither the composition`, C1, `nor the composition`, C2):
print(``):
print(` Then `):
print(``):
print(Sum(a(n)*x^n,n=0..infinity)=gu):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu):
print(``):

lu:=AsyGF(gu,x):
if lu<>FAIL then
 print(`The asymptotic  expression for a(n) is`, lu[2]*lu[1]^n ):
fi:
print(``):
gu:=GFsetX({C1,C2},x,X):
gu:=subs({X[op(C1)]=X1, X[op(C2)]=X2},gu):
print(``):
print(` Let `, A(n,c,d), ` be the number of compositions of n with c occurrences (as containment) of the composition `, C1):
print(`and d occurrences (as containment) of the composition`, C2, ` (Note that A(n,0,0)=a(n)), Then `):
print(``):
print(Sum(Sum(Sum(A(n,c,d)*x^n*X1^c*X2^d,c=0..infinity),d=0..infinity),n=0..infinity) =gu):
print(``):
print(`and in Maple format`):
print(``):
print(``):
lprint(gu):
print(``):

gu:=MomsKA2st(C1,C2,n,r):
if gu<>FAIL then
print(``):
print(`Furthermore the average and variance of the random variable: Number of occurrences of`, C1, `are `):
print(``):
print(gu[1][1] , `and `, gu[1][2] , `respectively, while `):
print(``):
print(`Furthermore the average and variance of the random variable: Number of occurrences of`, C2, `are `):
print(``):
print(gu[2][1] , `and `, gu[2][2] , `respectively, while  the asymptotic correlation between these two random variables is`):
print(``):

print(``):
print(gu[3][1][1]):
print(``):
print(`and in floating point`):
print(``):
print(evalf(gu[3][1][1])):
print(``):

mu:=ND2mTable(gu[3][1][1],r):

if mu<>FAIL then
if {seq(seq(simplify(mu[i][j]-gu[3][i][j]),j=1..r),i=1..r)}<>{0} then
 print(`Something amazing happened, the two random variables are not asymptotically normal`):
 RETURN(gu[3],mu):
else
print(`The asymptotic standardized mixed moments are indeed as they are for bi-variate normal pair with correlation`, mu[1][1]):
print(``):
print(`i.e. `, mu):
fi:

fi:
fi:

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


end:



#SuperInfoX2Vone(m1,m2,x,X1,X2,n,r):inputs  positive integers m1, m2, (m2<m1)
# continuous variable names x,X1,X2, discrete variable name n, and
#a positive integer r >=2, outputs an article about the joint generating functions of all pairs of
#of compositions of m1 into m2 parts. Try:
#SuperInfoX2Vone(5,2,x,X1,X2,n,4);
SuperInfoX2Vone:=proc(m1,m2,x,X1,X2,n,r) local co,gu,i,j,t0:

t0:=time():
if not (type(m1, integer) and m1>=2 and m2>=1 and m2<=m1 and type(x,symbol) and type(X1,symbol) and type(X2,symbol) and type(n,symbol) and 
type(r,integer)  and r>=2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:
print(``):
print(`----------------------------------`):
print(``):  
print(`All the generating functions (and statistical information) Enumerating compositions by number of occurrences, as containments of all possible pairs`):
print(`of offending compositions of`, m1, `with `, m2, `parts `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
co:=0:

gu:=Comps1(m1,m2):

for i from 1 to nops(gu) do
 for j from i+1 to nops(gu) do
 co:=co+1:
  InfoX2Vth(gu[i] , gu[j],x,X1,X2,n,r,co):
 od:
od:

print(``):
print(`------------------------`):
print(``):
print(`This ends this article, that took`, time()-t0, ` to generate. `):
print(``):
end:


#SuperInfoX2Vtwo(m1a,m1b, m2,x,X1,X2,n,r):inputs  positive integers m1a, ma1bm, m2, (m2<min(m1a,m1b))
# continuous variable names x,X1,X2, discrete variable name n, and
#a positive integer r >=2, outputs an article about the joint generating functions of all pairs of
#of compositions of m1a into m2 parts and compositions of m1b into m2 parts that are not included in each other. Try:
#SuperInfoX2Vtwo(5,6,2,x,X1,X2,n,2);
SuperInfoX2Vtwo:=proc(m1a,m1b, m2,x,X1,X2,n,r) local co,gu1,gu2,i,j,t0:

t0:=time():
if not (type(m1a, integer) and type(m1a, integer) 
and m1a>=2 and m1b>=2 and m2>=1 and m2<=min(m1a,m1b) and type(x,symbol) and type(X1,symbol) and type(X2,symbol) and type(n,symbol) and 
type(r,integer)  and r>=2) then
 print(`Bad input`):
 RETURN(FAIL):
fi:
print(``):
print(`----------------------------------`):
print(``):  
print(`All the generating functions (and statistical information) Enumerating compositions by number of occurrences, as containments of all possible pairs`):
print(`of offending compositions of`, m1a, `with `, m2, `parts and of `, m1b, `also with `, m2, ` parts `  ):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
co:=0:

gu1:=Comps1(m1a,m2):
gu2:=Comps1(m1b,m2):

for i from 1 to nops(gu1) do
 for j from 1 to nops(gu2) do
 if nops(Reduce1({gu1[i],gu2[j]}))=2 then

 co:=co+1:
  InfoX2Vth(gu1[i] , gu2[j],x,X1,X2,n,r,co):
 fi:
 od:
od:

print(``):
print(`------------------------`):
print(``):
print(`This ends this article, that took`, time()-t0, ` to generate. `):
print(``):
end: