######################################################################
##SERGI: Save this file as  SERGI                                    #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read SERGI<Enter>                                  #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Sept. 2010
#This version: Jan. 13 , 2011 
print(`Created:  Sept. 2010`):
print(`This version:  Jan. 13, 2011`):
print(` This is SERGI `):
print(`It is one of two Maple packages accompanying `):
print(`the article by Andrew Baxter, Brian Nakamura and Doron Zeilberger`):
print(`"Automatic Generation of Theorems and Proofs on `):
print(`Enumerating Consecutive-Wilf classes "`):
print(`available from the authors' websites`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 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 procedure, type ezra(procedure_name);   .`):
 print(``):
 print(`For a list of the ad-hoc procedures for computing, for help with`):
 print(`special families using formulas of Elizalde-Noy, Elizalde,`):
 print(` and Dotsenko-Khoroshkin  type ezraS() `):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(` AnkABrute, AnkAx, AnkAxBrute,AnkBrute, `):
 print(`ApplyUmbra, ApplyUmbralMatrix, CleanUpU, ECGF, findrec, Findrec`):
 print(` Init, `):
 print(` NiceSU, PUAt, PUAtV, PU1, PU1At, , redu , `):
 print(` Sader, Reps, SeqFromRec, SimpU`):
 print(`  SuperApplyUmbra, ToUmbra, U1At, U1AtV, Wt, Wtx`):

else
ezra(args):
fi:

end:


ezraS:=proc()

if args=NULL then
 print(` The MAIN special-purpose procedures are: `):
 print(` AnPam , AnPDK, AnPincr, AnP1324, AnP1423, AnP2413 `):

else
ezra(args):
fi:

end:


ezraS1:=proc()

if args=NULL then
 print(` The supporting special-purpose procedures are: `):
 print(` DK1324cn, DK1324cnl, DKfk, DKfkp `):
print(`DK1423cn, DK1423cnl,`):
else
ezra(args):
fi:

end:


ezraOld:=proc()

if args=NULL then
 print(`The main procedures are: `):
 print(`   UA, UAt, UAtV, UAtVsimp,  `):
else
 ezra(args):
fi:
end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: `):
 print(` AllSeqs, Ank, AnP, AnPt, AnPlus, AnPW, Theorem, Theoremt`):
 print(`Um, Umt, UmNice, UmNicet, UmS  `)
:


elif nops([args])=1 and op(1,[args])=AllSeqs then
print(`AllSeqs(k,r,N): All possible sequences enumerating`):
print(` r-sets of patterns each of length N for each`):
print(`representative class, sorted according to the last (N-th element)`):
print(`For example, try:`):
print(`AllSeqs(3,1,30);`):



elif nops([args])=1 and op(1,[args])=Ank then
print(`Ank(n,k,t): the weight-enumerator of all permutations of`):
print(`length n with weight being product of t[redu([p[i],..,p[i+k-1]])]`):
print(`for i=1..n-k+1 by the umbral method. For example try:`):
print(`Ank(4,3,t);`):

elif nops([args])=1 and op(1,[args])=AnkAx then
print(`AnkAx(n,k,t,A,x): the weight-enumerator of all permutations of`):
print(`length n with weight(p) and indicated states`):
print(`product of t[redu([p[i],..,p[i+k-1]])]`):
print(`for i=1..n-k+1 by the umbral appraoch. For example try:`):
print(`AnkAx(4,3,t,A,x);`):



elif nops([args])=1 and op(1,[args])=AnkMx then
print(`AnkMx(n,k,t,x,z): the weight-enumerator of all permutations of`):
print(`length n with weight(p) via matrix approach`):
print(`product of t[redu([p[i],..,p[i+k-1]])]`):
print(`for i=1..n-k+1 by the umbral matrix matrix. For example try:`):
print(`AnkMx(4,3,t,x,z);`):


elif nops([args])=1 and op(1,[args])=AnkM then
print(`AnkM(n,k,t): the weight-enumerator of all permutations of`):
print(`length n with weight(p) via matrix approach`):
print(`product of t[redu([p[i],..,p[i+k-1]])]`):
print(`for i=1..n-k+1 by the umbral matrix matrix. For example try:`):
print(`AnkM(4,3,t);`):

elif nops([args])=1 and op(1,[args])=AnkBrute then
print(`AnkBrute(n,k,t): the weight-enumerator of all permutations of`):
print(`length n with weight(p) product of t[redu([p[i],..,p[i+k-1]])]`):
print(`for i=1..n-k+1 by complete bruce-force. For example try:`):
print(`AnkBrute(4,3,t);`):

elif nops([args])=1 and op(1,[args])=AnkAxBrute then
print(`AnkAxBrute(n,k,t,A,x): the weight-enumerator of all permutations of`):
print(`length n with weight(p) product of t[redu([p[i],..,p[i+k-1]])]`):
print(`for i=1..n-k+1  and product of x[1]^i[1]*x[2]^i2*...`):
print(`where i1, i2 are the last entries of the permutation`):
print(`by complete bruce-force. For example try:`):
print(`AnkAxBrute(4,3,t,A,x);`):

elif nops([args])=1 and op(1,[args])=AnkABrute then
print(`AnkABrute(n,k,t,A): the weight-enumerator of all permutations of`):
print(`length n with weight(p) product of t[redu([p[i],..,p[i+k-1]])]`):
print(`for i=1..n-k+1 , and idicated states given by the indexed variable`):
print(` A, by complete bruce-force. For example try:`):
print(`AnkABrute(4,3,t,A);`):



elif nops([args])=1 and op(1,[args])=AnPincr then
print(`AnPincr(a,N): The first N terms of the sequence`):
print(`enumerating the weight-enumerators of permutations of length n`):
print(`according to the weight enumerator for`):
print(`occurrences of the  consecutive pattern [1,2,3,.., a]`):
print(`using the formula of Sergi Elizalde and Marc Noy. For exaple, try:`):
print(`AnPincr(3,20);`):

elif nops([args])=1 and op(1,[args])=AnPold then
print(`AnPold(N,FP): The first N terms of the sequence`):
print(`enumerating the number of permutations of length n`):
print(`avoiding the set of patterns FP`):
print(`For example, try:`):
print(`AnPold(10,{[1,2,3]});`):

elif nops([args])=1 and op(1,[args])=AnP then
print(`AnP(FP,N): The first N terms of the sequence`):
print(`enumerating the number of permutations of length n`):
print(`avoiding the set of patterns FP`):
print(`For example, try:`):
print(`AnP({[1,2,3]},10);`):


elif nops([args])=1 and op(1,[args])=AnPam then
print(`AnPam(a,m,N,t): The first N terms of the sequence`):
print(`enumerating the weight-enumerators of permutations of length n`):
print(`according to the number of`):
print(`occurrences of the  consecutive pattern 1...a tau (a+1) in S_{m+1}`):
print(`for any tau.`):
print(`For example, try:`):
print(`AnPam(1,2,10,t);`):

elif nops([args])=1 and op(1,[args])=AnPDK then
print(`AnPDK(a,b,m,N,t): The first N terms of the sequence`):
print(`enumerating the weight-enumerators of permutations of length n`):
print(`according to the weight enumerator for`):
print(`occurrences of the  consecutive pattern tau of length m+1`):
print(`with exactly one overlap and`):
print(`tau[1]=a+1 and tau[m+1]=b+1 following the formula`):
print(`in Vladimir Dotsenko and Anton Khoroshkin's `):
print(`paper "Anick-type resolutions and consecutive pattern`):
print(`avoidance."`):
print(`For example, try:`):
print(`AnPDK(1,2,3,10,t);`):

elif nops([args])=1 and op(1,[args])=AnP1324 then
print(`AnP1324(N,t): The first N terms of the sequence`):
print(`of weight-enumerators of permutations of length n`):
print(`according to the weight enumerator for`):
print(`occurrences of the  consecutive pattern 1324`):
print(`taken from Theorem 4 in`):
print(`Vladimir Dotsenko and Anton Khoroshkin's `);
print(`paper "Anick-type resolutions and consecutive pattern`):
print(`avoidance." For examle, try:`):
print(`AnP1324(10,t);`):

elif nops([args])=1 and op(1,[args])=AnP1423 then
print(`AnP1423(N,t): The first N terms of the sequence`):
print(`of weight-enumerators of permutations of length n`):
print(`according to the weight enumerator for`):
print(`occurrences of the  consecutive pattern 1423`):
print(`taken from Theorem 5 in`):
print(`Vladimir Dotsenko and Anton Khoroshkin's `);
print(`paper "Anick-type resolutions and consecutive pattern`):
print(`avoidance." For examle, try:`):
print(`AnP1423(10,t);`):


elif nops([args])=1 and op(1,[args])=AnP2413 then
print(`AnP2413(N,t): The first N terms of the sequence`):
print(`of weight-enumerators of permutations of length n`):
print(`according to the weight enumerator for`):
print(`occurrences of the  consecutive pattern 2413`):
print(`taken from Theorem 7 in`):
print(`Vladimir Dotsenko and Anton Khoroshkin's `);
print(`paper "Anick-type resolutions and consecutive pattern`):
print(`avoidance." For examle, try:`):
print(`AnP2413(10,t);`):

elif nops([args])=1 and op(1,[args])=AnPlus then
print(`AnPlus(N,FP,N1,C,n,E): The first N1 terms of the sequence`):
print(`enumerating the number of permutations of length n`):
print(`avoiding the set of patterns FP using the first N terms`):
print(`and conjecturing a linear recurrence of complexity <=C`):
print(`satisfied by`):
print(`the coefficients of the cluster generating function`):
print(`For example, try:`):
print(`AnPlus(20,{[1,2,3]},40,6,n,E);`):

elif nops([args])=1 and op(1,[args])=AnPlusBrian then
print(`AnPlusBrian(N,FP,N1,C,n,E): Like AnPlus but using`):
print(`Brian Nakamura's faster program`):
print(`The first N1 terms of the sequence`):
print(`enumerating the number of permutations of length n`):
print(`avoiding the set of patterns FP using the first N terms`):
print(`and conjecturing a linear recurrence of complexity <=C`):
print(`satisfied by`):
print(`the coefficients of the cluster generating function`):
print(`For example, try:`):
print(`AnPlusBrian(20,{[1,2,3]},40,6,n,E);`):

elif nops([args])=1 and op(1,[args])=AnPt then
print(`AnPt(FP,N,t): The first N terms of the sequence`):
print(`of weight-enumerators of permutations of length n`):
print(`according to the weight product of t[p]#occurrences of p in FP.`):
print(`For example, try:`):
print(`AnPt({[1,2,3]},10,t);`):

elif nops([args])=1 and op(1,[args])=AnPW then
print(`AnPW(N,FP,accur1): The first N terms of the sequence`):
print(`enumerating the number of permutations of length n`):
print(`avoiding the set of patterns FP plus the`):
print(`Warlimont constant (if it seems to exist) rho`):
print(`and the const. A such that the n-th term`):
print(`of the sequence is A*rho^n . accur1 is the required accuracy `):
print(`(number of digits).`):
print(`For example, try:`):
print(`AnPW(40,{[1,2,3]},10);`):

elif nops([args])=1 and op(1,[args])=ApplyUmbra then
print(`ApplyUmbra(Umbra1,f,ylist): given an umbra, Umbra1, applies`):
print(` it to the expression f in the variables in ylist`):


elif nops([args])=1 and op(1,[args])=ApplyUmbralMatrix then
print(`ApplyUmbralMatrix(UM,f,ylist): applies the umbral matrix UM`):
print(`to the vector of functions in the variables ylist.`):
print(`For example, try:`):
print(`ApplyUmbralMatrix(UAtV(2,x,z,t),Init(2,x,z),[x[1],x[2],z]);`):

elif nops([args])=1 and op(1,[args])=CleanUpU then
print(` CleanUpU(U): cleans up the Umbra U `):

elif nops([args])=1 and op(1,[args])=DKfk then
print(`DKfk(k,m,a,b): The number of k-clusters of length m*k+1`):
print(`for a pattern tau of length m+1 that starts with`):
print(`a+1 and ends with b+1 in`):
print(`Anton Khoroshkin's paper "Anick-type resolutions and consecutive pattern`):
print(`avoidance. `):
print(`For example, try:`):
print(`DKfk(2,3,1,2);`):

elif nops([args])=1 and op(1,[args])=DKfkp then
print(`DKfkp(k,p,m,a,b): `):
print(`implementing the recursion (5) of Vladimir Dotsenko and`):
print(`Anton Khoroshkin's paper "Anick-type resolutions and `):
print(`consecutive pattern avoidance".`):
print(`It is the number of k-clusters for any consecutive pattern`):
print(`that starts with a+1 and ends at b+1 such that the cluster`):
print(`starts at p+1. `):
print(` For example, try:`):
print(`DKfkp(2,6,3,1,2);`):


elif nops([args])=1 and op(1,[args])=DK1324cnl then
print(`DK1324cnl(n,l): Formula (6) `):
print(`in Vladimir Dotsenko and Anton Khoroshkin's `):
print(`paper "Anick-type resolutions and consecutive pattern`):
print(`avoidance." For example, try:`):
print(`DK1324cnl(10,3);`):

elif nops([args])=1 and op(1,[args])=DK1324cn then
print(`DK1324cn(n,t): DK1324cnl(n,l)*(t-1)^l for l from 1 to n`):
print(`For example, try:`):
print(`DK1324cn(10,t);`):


elif nops([args])=1 and op(1,[args])=ECGF then
print(`ECGF(FP,K): The first K terms of the empirical `):
print(`cluster exponential generating function. For example, try:`):
print(`ECGF({[1,3,2]},20);`):

elif nargs=1 and args[1]=findrec then
print(`findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating`):
print(`the sequence f of degree DEGREE and order ORDER.`):
print(`For example, try: findrec([seq(i,i=1..10)],0,2,n,N);`):

elif nargs=1 and args[1]=Findrec then
print(`Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with`):
print(`poly coffs. ope(n,N), where n is the discrete variable, and N is the shift operator `):
print(`of maximum DEGREE+ORDER<=MaxC`):
print(`e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);`):

elif nargs=1 and args[1]=FindrecF then
print(`FindrecF(f,n,N): Given a function f of a single variable tries to find a linear recurrence equation with`):
print(`poly coffs. .g. try FindrecF(i->i!,n,N);`):


elif nargs=1 and args[1]=SeqFromRec then
print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`):
print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`):
print(`satisfied by the recurrence ope(n,N)f(n)=0`):
print(`extends it to the first K values`):
print(`For example, try:`):
print(`SeqFromRec(N-n-1,n,N,[1],10);`):

elif nops([args])=1 and op(1,[args])=Init then
print(`Init(k,x,z): The initial generating function for permutations`):
print(`of length x. For example, try:`):
print(`Init(2,x,z);`):

elif nops([args])=1 and op(1,[args])=NiceSU then
print(`NiceSU(U,F): A nice rendition of the umbral operator U`):
print(`using the function F. For example, try:`):
print(`NiceSU(UA(x,z,A,{[1,2,3]}),F);`):

elif nops([args])=1 and op(1,[args])=PUAt then
print(`PUAt(pi,x,i,n,A,t): Given a permutation pi (of size k-1, say), `):
print(`and sympols x,i,n,A,t, gives the pre-umbra`):
print(`that transforms a permutation in final state pi`):
print(`to any other feasible state `):
print(`For example, try:`):
print(`PUAt([2,3,1],x,i,n,A,t);`):

elif nops([args])=1 and op(1,[args])=PUAtV then
print(`PUAtV(pi,x,i,n,A,t): The vector form of  PUAt. For example, try:`):
print(`PUAtV([2,3,1],x,i,n,t);`):


elif nops([args])=1 and op(1,[args])=PU1 then
print(`PU1(pi,r,x,i,n): Given a permutation pi (of size k-1, say), `):
print(`a half integer r (1/2, ...,k/2) outputs the pre-umbra`):
print(`that transforms a permutation in final state pi`):
print(`to a permutation of final state redu([pi[2], ..., pi[k],r])`):
print(`by appending (and making room) at the very end`):
print(`an entry between i[r-1/2]+1 and i[r+1/2]`):
print(`with i[0]=0 and i[k]=n+1`):
print(`For example, try:`):
print(`PU1([2,3,1],3/2,x,i,n);`):

elif nops([args])=1 and op(1,[args])=PU1At then
print(`PU1At(pi,r,x,i,n,A,t): Given a permutation pi (of size k-1, say), `):
print(`a half integer r (1/2, ...,k/2) outputs the pre-umbra`):
print(`that transforms a permutation in final state pi`):
print(`to a permutation of final state redu([pi[2], ..., pi[k],r])`):
print(`by appending (and making room) at the very end`):
print(`an entry between i[r-1/2]+1 and i[r+1/2]`):
print(`with i[0]=0 and i[k]=n+1`):
print(`Here it includes the new state and the earned brownie point`):
print(`For example, try:`):
print(`PU1At([2,3,1],3/2,x,i,n,A,t);`):

elif nops([args])=1 and op(1,[args])=redu then
print(`redu(perm): Given a permutation on a set of positive integers`):
print(`finds its reduced form.`):
print(`For example, try:`):
print(` redu([16,24,3]); `):

elif nops([args])=1 and op(1,[args])=Reps then
print(`Reps(k,r): a set of representatives for the `):
print(`sets of sets of r patterns all of length r`):
print(`For example, try: `):
print(`Reps(3,2);`):

elif nops([args])=1 and op(1,[args])=Sader then
print(`Sader(L): inputs a list L and outputs the list whose`):
print(`i-th item is the i-th largest element`):
print(`For example, try:`):
print(`Sader([3,1,2,3,2]);`):

elif nops([args])=1 and op(1,[args])=SimpU then
print(`SimpU(U): simplifying an umbral operator that has no differentiation`):
print(`to simpler format {[coeff,argument]}.`):
print(`For example, try:`):
print(`SimpU({[1,[0,0,0],[x,y,z]]});`):

elif nops([args])=1 and op(1,[args])=SuperApplyUmbra then
print(`SuperApplyUmbra(U,f,x,z,k,A): Given a super-umbra U with states`):
print(`parametrized by k-permutation and a polynomial f in`):
print(`x[1],..., x[k],z and a linear expression in A[pi] with`):
print(`pi k-perms, applies U to f. For example, try:`):
print(`SuperApplyUmbra(UAt(1,x,z,A,t),x[1]*z*A[1],x,z,1,A);`):


elif nops([args])=1 and op(1,[args])=Theorem then
print(`Theorem(FP): The theorem that enables fast enumeration`):
print(`of permutations avoiding the set of patterns FP`):
print(`of the same length.`):
print(`For example, try:`):
print(`Theorem({[1,2,3]});`):

elif nops([args])=1 and op(1,[args])=Theoremt then

print(`Theoremt(k,t): The theorem that enables fast `):
print(`weighted enumeration according to the weight`):
print(`keep track of all consecutive patterns of length k`):
print(`t is the variable for indexing the consecutive patters.`):
print(`For example, try:`):
print(`Theoremt(3,t);`):


elif nops([args])=1 and op(1,[args])=ToUmbra then
print(` ToUmbra(poly1,xlist,alist): given a polynomial, poly1, in the`):
print(` variables xlist with exponents that are affine-linear`):
print(` expressions in the discrete variables in the`):
print(` list alist, and in the variables of alist themselves`):
print(` outputs the corresponding Umbra, such that`):
print(` poly1 is the image, under that umbra, of the`):
print(` generic monomial y[1]^alist[1]*...*y[k]^alist[k],`):
print(` (where k=nops(alist), and y[1], ..., y[k] are generic`):
print(` continuous variables that correspond to the disrete variables`):
print(` in alist`):
print(` The format of the output is a set, each element of which`):
print(` is a list of 3-elements`):
print(` [FRONT, Diffs,SubsList], where the FRONT`):
print(` is a rational function in whatever, Diffs is a list of`):
print(` integers of the same length as alist, and SubsList is`):
print(` the list of substitutions that the continuous variables`):
print(` that correspond to alist have to be substituted by`):
print(` For example ToUmbra(a*x^b+b*y^a,[x,y],[a,b]) should yield`):
print(`  {[1, [1, 0], [1, x]], [1, [0, 1], [y, 1]]}`):


elif nops([args])=1 and op(1,[args])=UA then
print(`UA(x,z,A,FP): The global`):
print(`umbral evolution operator, for enumerating`):
print(`permutations avoiding the patterns (of the same length) in FP`):
print(`phrased in terms of the indexed variables`):
print(`x[1],x[2],..x[k] and  the variable z (for the size of the permutation`):
print(`For example, try:`):
print(`UA(x,z,A,{[1,2,3]});`):

elif nops([args])=1 and op(1,[args])=UAt then
print(`UAt(k,x,z,A,t): The global`):
print(`umbral evolution operator, for all consecutive patterns`):
print(`of length k+1`):
print(`phrased in terms of the indexed variables`):
print(`x[1],x[2],..x[k] and  the variable z (for the size of the permutation`):
print(`and the index variable A keeping track of the states.`):
print(`For example, try:`):
print(` UAt(2,x,z,A,t); `):

elif nops([args])=1 and op(1,[args])=UAtV then
print(`UAtV(k,x,z,t): The global`):
print(`umbral evolution matrix operator, for all consecutive patterns`):
print(`of length k+1`):
print(`phrased in terms of the indexed variables`):
print(`x[1],x[2],..x[k] and  the variable z (for the size of the permutation`):
print(`For example, try:`):
print(` UAtV(2,x,z,t); `):

elif nops([args])=1 and op(1,[args])=UAtVsimp then
print(`UAtVsimp(k,x,z,t): The simplifed form of the global`):
print(`umbral evolution matrix operator, for all consecutive patterns`):
print(`of length k+1`):
print(`phrased in terms of the indexed variables`):
print(`x[1],x[2],..x[k] and  the variable z (for the size of the permutation`):
print(`For example, try:`):
print(` UAtVsimp(2,x,z,t); `):

elif nops([args])=1 and op(1,[args])=U1At then
print(`U1At(pi,x,z,A,t): Given a permutation pi (of size k, say), `):
print(`outputs the umbra, phrased in terms of the indexed variables`):
print(`x[1],x[2],..x[k] and  the variable z (for the size of the permutation`):
print(`that transforms a permutation in final state pi`):
print(`by appending (and making room) at the very end`):
print(`any entry`):
print(` Here it includes the new state and the`):
print(`earned brownie point. A is a variable for indexing the states`):
print(`For example, try:`):
print(`U1At([2,3,1],x,z,A,t);`):


elif nops([args])=1 and op(1,[args])=U1AtV then
print(`U1AtV(pi,x,z,t): Given a permutation pi (of size k, say), `):
print(`outputs the umbral , vector of size nops(pi)! `):
print(`phrased in terms of the indexed variables`):
print(`x[1],x[2],..x[k] and  the variable z (for the size of the permutation`):
print(`that transforms a permutation in final state pi`):
print(`For example, try:`):
print(`U1AtV([2,3,1],x,z,t);`):


elif nops([args])=1 and op(1,[args])=Um then
print(`Um(x,z,FP): The global`):
print(`umbral evolution matrix operator, for enumerating`):
print(`permutations avoiding the patterns (of the same length) in FP`):
print(`in terms of an umbral matrix and catalytic variables`):
print(`x[1],x[2],..x[k] and  the variable z `):
print(`(for the size of the permutation)`):
print(`For example, try:`):
print(`Um(x,z,{[1,2,3]});`):

elif nops([args])=1 and op(1,[args])=Umt then
print(`Umt(x,z,FP,t): The global`):
print(`umbral evolution matrix operator, for enumerating`):
print(`permutations avoiding the patterns (of the same length) in FP`):
print(`in terms of an umbral matrix and catalytic variables`):
print(`x[1],x[2],..x[k] and  the variable z `):
print(`where the weight is the prouduct of t[p]#(occurences of p) `):
print(`over all partterns in FP`):
print(`(for the size of the permutation)`):
print(`For example, try:`):
print(`Umt(x,z,{[1,2,3]},t);`):

elif nops([args])=1 and op(1,[args])=UmNice then
print(`UmNice(x,z,FP,f,g): A nice rendition of `):
print(`an Umbral matrix operator with the input function`):
print(`being called f[i] and the output function called g[i]`):
print(`For example, try:`):
print(`UmNice(x,z,{[1,2,3]},f,g);`):

elif nops([args])=1 and op(1,[args])=UmNicet then
print(`UmNicet(x,z,k,t,f,g): A nice rendition of `):
print(`an Umbral matrix operator with the input function`):
print(`being called f[i] and the output function called g[i]`):
print(`for the general problem taking care of all patterns of`):
print(`length k`):
print(`For example, try:`):
print(`UmNicet(x,z,3,t,f,g);`):

elif nops([args])=1 and op(1,[args])=UmS then
print(`UmS(x,z,FP): The simplified form (without differentiations `):
print(`that don't show up here of the global`):
print(`umbral evolution matrix operator, for enumerating`):
print(`permutations avoiding the patterns (of the same length) in FP`):
print(`in terms of an umbral matrix and catalytic variables`):
print(`x[1],x[2],..x[k] and  the variable z `):
print(`(for the size of the permutation)`):
print(`For example, try:`):
print(`UmS(x,z,{[1,2,3]});`):

elif nops([args])=1 and op(1,[args])=Wt then
print(`Wt(pi,k,t): the weight of a permutation pi of level-k`):
print(`the product of t[] of all consecutive patterns.`):
print(`For example, try:`):
print(`Wt([2,4,1,3,5],3,t);`):

elif nops([args])=1 and op(1,[args])=Wtx then
print(`Wtx(pi,k,t,x): the weight of a permutation pi of level-k`):
print(`the product of t[] of all consecutive patterns`):
print(`and the x's of the last entries.`):
print(`For example, try:`):
print(`Wtx([2,4,1,3,5],3,t,x);`):

else
print(`There is no ezra for`,args):
fi:
 
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:

with(combinat):

#Wt(pi,k,t): the weight of a permutation pi of level-k
#the product of t[] of all consecutive patterns
#For example, try:
#Wt([2,4,1,3,5],3,t);
Wt:=proc(pi,k,t) local i:
mul(t[op(redu([op(i..i+k-1,pi)]))],i=1..nops(pi)-k+1):
end:

#Wtx(pi,k,t,x): the weight of a permutation pi of level-k
#the product of t[] of all consecutive patterns
#and the x's of the last entries.
#For example, try:
#Wtx([2,4,1,3,5],3,t,x);
Wtx:=proc(pi,k,t,x) local i,sig,n:
n:=nops(pi):
sig:=sort([op(n-k+1..n,pi)]):
mul(t[op(redu([op(i..i+k,pi)]))],i=1..nops(pi)-k)
*mul(x[i]^sig[i],i=1..nops(sig)):
end:

#AnkBrute(n,k,t): the weight-enumerator of all permutations of
#length n with weight(p) product of t[redu([p[i],..,p[i+k-1]])]
#for i=1..n-k+1 by complete bruce-force. For example try:
#AnkBrute(4,3,t);
AnkBrute:=proc(n,k,t) local gu,pi:
gu:=permute(n):
add(Wt(pi,k,t),pi in gu):
end:




#AnkABrute(n,k,t,A): the weight-enumerator of all permutations of
#length n with weight(p) and indicated states
#product of t[redu([p[i],..,p[i+k-1]])]
#for i=1..n-k+1 by complete bruce-force. For example try:
#AnkABrute(4,3,t,A);
AnkABrute:=proc(n,k,t,A) local gu,pi:
gu:=permute(n):
add(Wt(pi,k,t)*A[op(redu([op(n-k+2..n,pi)]))],pi in gu):
end:

#AnkAxBrute(n,k,t,A,x): the weight-enumerator of all permutations of
#length n with weight(p) and indicated states
#product of t[redu([p[i],..,p[i+k-1]])] and prouct of x[last_entries]
#for i=1..n-k+1 by complete bruce-force. For example try:
#AnkAxBrute(4,3,t,A,x);
AnkAxBrute:=proc(n,k,t,A,x) local gu,pi:
gu:=permute(n):
add(Wtx(pi,k,t,x)*A[op(redu([op(n-k+1..n,pi)]))],pi in gu):
end:


#PU1(pi,r,x,i,n): Given a permutation pi (of size k-1, say), 
#a half integer r (1/2, ...,k/2) outputs the pre-umbra
#that transforms a permutation in final state pi
#to a permutation of final state redu([pi[2], ..., pi[k],r])
#by appending (and making room) at the very end
#an entry between i[r-1/2]+1 and i[r+1/2]
#with i[0]=0 and i[k]=n+1
#For example, try:
#PU1([2,3,1],3/2,x,i,n);
PU1:=proc(pi,r,x,i,n) local lu,i1,j,k,npi,nnpi:
k:=nops(pi):
lu:=mul(x[pi[i1]]^i[pi[i1]],i1=1..k):

for i1 from 1 to k do
if pi[i1]>r then
 lu:=x[pi[i1]]*lu:
fi:
od:

lu:=lu*sum(x[r]^j,j=i[r-1/2]+1..i[r+1/2]):
lu:=subs(i[0]=0,lu):
lu:=subs(i[k+1]=n+1,lu):


npi:=redu([op(pi),r]):

lu:=subs({seq(x[pi[i1]]=x[npi[i1]],i1=1..nops(pi)), x[r]=x[r+1/2]},lu):

nnpi:=redu([op(2..nops(npi),npi)]):

lu:=subs({x[npi[1]]=1, seq(x[npi[i1]]=x[nnpi[i1-1]],i1=2..nops(npi))},lu):
expand(lu):
end:


#PU1At(pi,r,x,i,n,A,t): Given a permutation pi (of size k-1, say), 
#a half integer r (1/2, ...,k/2) outputs the pre-umbra
#that transforms a permutation in final state pi
#to a permutation of final state redu([pi[2], ..., pi[k],r])
#by appending (and making room) at the very end
#an entry between i[r-1/2]+1 and i[r+1/2]
#with i[0]=0 and i[k]=n+1. Here it includes the new state and the
#earned brownie point
#For example, try:
#PU1At([2,3,1],3/2,x,i,n,A,t);
PU1At:=proc(pi,r,x,i,n,A,t) local gu,pi1,pi2:
gu:=PU1(pi,r,x,i,n):
pi1:=redu([op(pi),r]):
pi2:=redu([op(2..nops(pi1),pi1)]):
gu:=expand(gu*A[op(redu(pi2))]*t[op(pi1)]):
end:




#PUAt(pi,x,i,n,A,t): Given a permutation pi (of size k-1, say), 
#outputs the pre-umbra
#that transforms a permutation in final state pi
#by appending (and making room) at the very end
#any entry
# Here it includes the new state and the
#earned brownie point
#For example, try:
#PUAt([2,3,1],x,i,n,A,t);
PUAt:=proc(pi,x,i,n,A,t) local r:
expand(add(PU1At(pi,r-1/2,x,i,n,A,t),r=1..nops(pi)+1)):
end:

##############Stuff from ROTA################
# ToUmbra(poly1,xlist,alist): given a polynomial, poly1, in the
# variables xlist with exponents that are affine-linear
# expressions in the discrete variables in the
# list alist, and in the variables of alist themselves
# outputs the corresponding Umbra, such that
# poly1 is the image, under that umbra, of the
# generic monomial y[1]^alist[1]*...*y[k]^alist[k],
# (where k=nops(alist), and y[1], ..., y[k] are generic
# continuous variables that correspond to the disrete variables
# in alist
# The format of the output is a set, each element of whcih
# is a list of 3-elements
# [FRONT, Diffs,SubsList], where the FRONT
# is a rational function in whatever, Diffs is a list of
# integers of the same length as alist, and SubsList is
# the list of substitutions that the continuous variables
# that correspond to alist have to be substituted by
# For example ToUmbra(a*x^b+b*y^a,[x,y],[a,b]) should yield
#  {[1, [1, 0], [1, x]], [1, [0, 1], [y, 1]]}
ToUmbra:=proc(poly1,xlist,alist)
local gu,i,poly2:
if expand(poly1)=0 then
 RETURN(0):
fi:

poly2:=expand(poly1):
 
if not type(poly2,`+`) then
RETURN({UmbralTerm(poly2,xlist,alist)}):
fi:
 
gu:={}:
 
for i from 1 to nops(poly2) do
 gu:=gu union {UmbralTerm(op(i,poly2),xlist,alist)}:
od:
 
gu:
end:
 

# UmbralTerm(mono,xlist,alist): given a monomial in the
# variables xlist with exponents that are affine-linear
# expressions in the discrete variables in the
# list alist, and in the variables of alist themselves
# outputs the corresponding term in the Umbra
# The format of the output is a list of 3-elements
# [FRONT, Diffs,SubsList], where the FRONT
# is a rational function in whatever, Diffs is a list of
# integers of the same length as alist, and SubsList is
# the list of substitutions that the continuous variables
# that correspond to alist have to be substituted by
#
UmbralTerm:=proc(mono,xlist,alist)
local mono1,FRONT, Diffs, SubsList,i,d1,gu:
mono1:=expand(mono):
gu:=hafokh(mono1,alist):
FRONT:=gu[2]:
SubsList:=gu[1]:
 
mono1:=expand(FRONT):
Diffs:=[]:
 
for i from 1 to nops(alist) do
d1:=degree(mono1,alist[i]):
mono1:=expand(normal(expand(mono1/alist[i]^d1))):
Diffs:=[op(Diffs),d1]:
od:
 
[mono1,Diffs,SubsList]:
end:
 
#hafokh(mono,alist): given a monomial mono, in the form
#product(x[i]^a[j]) outputs the list [p[1],...p[k]] and the
#left-over monomial, lu, such that
#such that it equals lu*p[1]^a[1]*...*p[k]^a[k] times
hafokh:=proc(mono,alist)
local mono1,i,resh,mu:
mono1:=expand(mono):
resh:=[]:
for i from 1 to nops(alist) do
 mu:=grab(mono1,alist[i]):
 resh:=[op(resh),mu[1]]:
 mono1:=mu[2]:
od:
 
resh,mono1:
end:

#grab(mono,a): given a monomial, mono, in the variables 
#
# exponent (discrete) variable, a, returns
#the largest monomial, lu, such that lu^a is a factor of mono
grab:=proc(mono,a)
local mono1,mono2,i,lu,mu,mu1,mu2,khe,mu11,mu12:
mono1:=expand(mono):
mono2:=expand(mono1):
lu:=1:
 
if type(mono1,`*`)  then
 
for i from 1 to nops(mono1) do
mu:=op(i,mono1):
if type(mu,`^`) then
 mu1:=op(1,mu):
 mu2:=op(2,mu):
  if type(mu1,`^`) then
    mu11:=op(1,mu1):
     if type(mu11, `^`) then
      ERROR(`I give up`):
     fi:
    mu12:=op(2,mu1):
   mu1:=mu11:
   mu2:=mu2*mu12:
  fi:
 
 
 khe:=coeff(expand(mu2),a,1):
 lu:=lu*mu1^khe:
 mono2:=simplify(mono2/mu1^(a*khe)):
fi:
od:
 
elif type(mono1,`^`)  then
 mu1:=op(1,mu):
 mu2:=op(2,mu):
  if type(mu1,`^`) then
    mu11:=op(1,mu1):
     if type(mu11, `^`) then
      ERROR(`I give up`):
     fi:
    mu12:=op(2,mu1):
   mu1:=mu11:
   mu2:=mu2*mu12:
  fi:
 
 khe:=coeff(mu2,a,1):
 lu:=lu*mu1^khe:
 mono2:=simplify(mono2/mu1^(a*khe)):
fi:
lu,simplify(expand(mono2),symbolic):
end:


# ApplyUmbra(Umbra1,f,ylist): given an umbra, Umbra1, applies
# it to the expression f in the variables in ylist
ApplyUmbra:=proc(Umbra1,f,ylist)
local i,gu:
if Umbra1=0 then
 RETURN(0):
fi:
gu:=0:
 
for i from 1 to nops(Umbra1) do
 gu:=normal(gu+ApplyUmbraTerm(Umbra1[i],f,ylist)):
od:
expand(gu):
end:

# ApplyUmbraTerm(Umbra1,f,ylist): given an umbral term, Umbra1, applies
# it to the expression f in the variables in ylist
ApplyUmbraTerm:=proc(Umbra1,f,ylist)
local i, FRONT,Diffs,SubsList,gu,bu:
 
FRONT:=Umbra1[1]:
Diffs:=Umbra1[2]:
SubsList:=Umbra1[3]:
 
if nops(Diffs)<>nops(ylist) then
 ERROR(`Diffs and ylist should have the same length`):
fi:
 
 
if nops(SubsList)<>nops(ylist) then
 ERROR(`SubsList and ylist should have the same length`):
fi:
 
gu:=f:
 
for i from 1 to nops(ylist) do
gu:=normal(xdiff1(gu,ylist[i],Diffs[i])):
od:
 
bu:={}:
for i from 1 to nops(ylist) do
bu:=bu union {ylist[i]=SubsList[i]}:
od:
gu:=subs(bu,gu):
 
normal(FRONT*gu):
 
end:

# xdiff1(f,x,i): given an expression f, and a variable x,
# and an integer i, finds (xD_x)^i f
# 
xdiff1:=proc(f,x,i)
local gu,j:
 
if i=0 then
 RETURN(f):
fi:
 
gu:=f:
 
for j from 1 to i do
gu:=x*diff(gu,x):
od:
 
gu:
 
end:

 
##############End Stuff from ROTA################



#U1At(pi,x,z,A,t): Given a permutation pi (of size k, say), 
#outputs the umbra, phrased in terms of the indexed variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation
#that transforms a permutation in final state pi
#by appending (and making room) at the very end
#any entry
# Here it includes the new state and the
#earned brownie point. A is a variable for indexing the states
#For example, try:
#U1At([2,3,1],x,z,A,t);
U1At:=proc(pi,x,z, A,t) local   i,n,k,i1,gu:
k:=nops(pi):
gu:=PUAt(pi,x,i,n,A,t):
expand(ToUmbra(gu,[seq(x[i1],i1=1..k),z],[seq(i[i1],i1=1..k),n])):
end:



#UAt(k,x,z,A,t): The global
#umbral evolution operator, for all consecutive patterns
#of length k+1
#phrased in terms of the indexed variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation
#and the index variable A keeping track of the states
#For example, try:
#UAt(2,x,z,A,t);
UAt:=proc(k,x,z, A,t) local   gu,pi:
option remember:
gu:=permute(k):
{seq([pi,U1At(pi,x,z, A,t)],pi in gu)}:
end:

#SuperApplyUmbra(U,f,x,z,k,A): Given a super-umbra U with states
#parametrized by k-permutation and a polynomial f in
#x[1],..., x[k],z and a linear expression in A[pi] with
#pi k-perms, applies U to f. For example, try:
#SuperApplyUmbra(UAt(1,x,z,A,t),x[1]*z*A[1],x,z,1,A);
SuperApplyUmbra:=proc(U,f,x,z,k,A) local u,i1:

add(
expand(ApplyUmbra(u[2],coeff(f,A[op(u[1])],1),[seq(x[i1],i1=1..k),z])),
u in U):
end:



#AnkAx(n,k,t,A,x): the weight-enumerator of all permutations of
#length n with weight(p) and indicated states
#product of t[redu([p[i],..,p[i+k-1]])]
#for i=1..n-k+1 by umbral approach
#AnkAx(4,3,t,A,x);
AnkAx:=proc(n,k,t,A,x) local z,gu,U:
option remember:

if n<=k then
RETURN(AnkAxBrute(n,k,t,A,x)):
fi:

U:=UAt(k,x,z, A,t):

gu:=expand(z^(n-1)*AnkAx(n-1,k,t,A,x)):

SuperApplyUmbra(U,gu,x,z,k,A):
end:


#Ank(n,k,t): the weight-enumerator of all permutations of
#length n with weight(p) and indicated states
#product of t[redu([p[i],..,p[i+k-1]])]
#for i=1..n-k+1 by complete bruce-force. For example try:
#Ank(4,3,t);
Ank:=proc(n,k,t) local gu,A,x,mu:
gu:=AnkAx(n,k-1,t,A,x):
mu:=permute(k-1):
subs({seq(x[i]=1,i=1..k-1), seq(A[op(pi)]=1,pi in mu)},gu):
end:


#CleanUpU(U): cleans up the Umbra U
CleanUpU:=proc(U) local U1,u:
U1:={}:

for u in U do
if expand(u[1])<>0 then
 U1:=U1 union {u}:
fi:
od:
U1:
end:

CleanUpU2:=proc(L) local i1,j1:
[seq([seq(CleanUpU(L[i1][j1]),j1=1..nops(L[i1]))],i1=1..nops(L))]:
end:

#UA(x,z,A,FP): The global
#umbral evolution operator, for enumerating
#permutations avoiding the patterns (of the same length) in FP
#phrased in terms of the indexed variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation
#For example, try:
#UA(x,z,A,{[1,2,3]});
UA:=proc(x,z,A,FP) local k,pi,S,t,lu,gu,i:
option remember:

if FP={} then
 print(`Set of patterns must be non-empty`):
fi:

S:={seq(nops(pi), pi in FP)}:


if nops(S)>1 then
 print(`All patterns must have the same length, not implemented`):
 RETURN(FAIL):
fi:

k:=S[1]:

gu:=UAt(k-1,x,z, A,t):

lu:=convert(permute(k),set) minus FP:

gu:=subs({seq(t[op(pi)]=1,pi in lu), seq(t[op(pi)]=0,pi in FP)},gu):

gu:={seq([gu[i][1],CleanUpU(gu[i][2])],i=1..nops(gu))}:
end:



#SimpU(U): simplifying an umbral operator that has no differentiation
#to simpler format {[coeff,argument]}
#For example, try:
#SimpU({[1,[0,0,0],[x,y,z]]});
SimpU:=proc(U) local u:
if U={} then
 RETURN({}):
fi:
if {seq(evalb(u[2]=[0$nops(u[2])]), u in U)}<>{true} then
 RETURN(FAIL):
fi:

{seq([u[1],u[3]], u in U)}:
end:




#NiceSU(U,F): A nice rendition of the umbral operator U
#using the function F. For example, try:
#NiceSU(UA(x,z,A,{[1,2,3]}),F);
NiceSU:=proc(U,F) local u:
add(u[1]*F(op(u[2])),u in U):
end:




#AnPold(N,FP): The first N terms of the sequence
#enumerating the number of permutations of length n
#avoiding the set of patterns FP
#For example, try:
#AnPold(10,{[1,2,3]});
AnPold:=proc(N,FP) local x,z,A,gu,pol,i,lu,pi,t,S,k,i1,mu,n,mu1,gu1:
lu:=UA(x,z,A,FP):

if FP={} then
 print(`Set of patterns must be non-empty`):
fi:

S:={seq(nops(pi), pi in FP)}:

if nops(S)>1 then
 print(`All patterns must have the same length, not implemented`):
 RETURN(FAIL):
fi:

k:=S[1]:

gu:=[seq(i!,i=1..k-1)]:
pol:=AnkAx(k-1,k-1,t,A,x):

mu:=convert(permute(k),set) minus FP:
mu1:=permute(k-1):


pol:=expand(
z^(k-1)*subs({seq(t[op(pi)]=1,pi in mu), seq(t[op(pi)]=0,pi in FP)},pol)):



for n from k to N do
 pol:=expand(z^n*SuperApplyUmbra(lu,pol,x,z,k-1,A)):
 gu1:=subs({z=1,seq(x[i1]=1,i1=1..k-1)},pol):
 gu1:=subs({seq(A[op(pi)]=1,pi in mu1)},gu1):
 gu:=[op(gu),gu1]:
od:

gu:

end:


############Start Trivial Equivalence###########

#Reps(k,r): all the representatives of single patterns of length k
#For example, try:
#Reps(4,1);

Reps:=proc(k,r) local gu,mu:
mu:=convert(permute(k),set):
mu:=choose(mu,r):
gu:={}:
while mu<>{} do
gu:=gu union {mu[1]}:
mu:=mu minus Khaverim(mu[1]):
od:
gu:
end:




############End Trivial Equivalence###########


#AnPW(N,FP,accur1): The first N terms of the sequence
#enumerating the number of permutations of length n
#avoiding the set of patterns FP plus the
#Warlimont constant (if it seems to exist) rho
#and the const. A such that the n-th term
#of the sequence is A*rho^n . accur1 is the required accuracy 
#(number of digits)
#For example, try:
#AnPW(40,{[1,2,3]},10);
AnPW:=proc(N,FP,accur1)  local gu,mu1,mu2,mu3,A1,A2,A3,rho,A:
Digits:=3*accur1:
gu:=AnP(FP,N):
mu1:=gu[N]/gu[N-1]/N:
mu2:=gu[N-1]/gu[N-2]/(N-1):
mu3:=gu[N-2]/gu[N-3]/(N-2):
 if abs(mu1-mu2)<10^(-accur1-3) and abs(mu2-mu3)<10^(-accur1-3) then
   rho:=evalf(mu1,accur1):
  else
  RETURN(gu,FAIL,FAIL):
 fi:


A1:=gu[N]/N!/rho^N:
A2:=gu[N-1]/(N-1)!/rho^(N-1):
A3:=gu[N-2]/(N-2)!/rho^(N-2):
print(A1,A2,A3):
 if abs(A1-A2)<10^(-accur1) and abs(A2-A3)<10^(-accur1) then
   A:=evalf(A1,accur1-2):
   RETURN(gu,rho,A):
  else
  RETURN(gu,rho,FAIL):
 fi:
end:


#Sader(L): inputs a list L and outputs the list whose
#i-th item is the i-th largest element
#For example, try:
#Sader([3,1,2,3,2]);
Sader:=proc(L) local L1,gu,gu1,i,j:
L1:=sort(convert(convert(L,set),list)):
gu:=[]:

for i from 1 to nops(L1) do
gu1:={}:
for j from 1 to nops(L) do
 if L[j]=L1[i] then
    gu1:=gu1 union {j}:
 fi:
od:
gu:=[op(gu),gu1]:
od:
gu:
end:

#AllSeqs(k,r,N): All possible sequences enumerating
# r-sets of patterns each of length N for each
#representative class, sorted according to the last (N-th element)
#For example, try:
#AllSeqs(3,1,30);
AllSeqs:=proc(k,r,N) local gu,i,lu,mu,lu1,i1:

gu:=convert(Reps(k,r),list):

lu:=[]:

for i from 1 to nops(gu) do
lu:=[op(lu),[gu[i],AnP(gu[i],N)]]:
od:

lu1:=[seq(lu[i][2][N],i=1..nops(lu))]:
lu1:=Sader(lu1):
lu1:=[seq(op(lu1[i1]),i1=1..nops(lu1))]:

mu:=[]:

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








###Start matrix form of Umbral operator########


#PUAtV(pi,x,i,n,A,t): The vector form of  PUAt. For example, try:
#PUAtV([2,3,1],x,i,n,t);
PUAtV:=proc(pi,x,i,n,t) local gu,A,mu,i1:
 gu:=PUAt(pi,x,i,n,A,t):
 mu:=permute(nops(pi)):
 [seq( coeff(gu,A[op(mu[i1])],1),i1=1..nops(mu))]:
end:



#U1AtV(pi,x,z,t): Given a permutation pi (of size k, say), 
#outputs the umbra, phrased in terms of the indexed variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation
#that transforms a permutation in final state pi
#by appending (and making room) at the very end
#any entry
# Here it includes the new state and the
#For example, try:
#U1AtV([2,3,1],x,z,t);
U1AtV:=proc(pi,x,z,t) local   gu,i1,k,j1,i,n:
k:=nops(pi):
gu:=PUAtV(pi,x,i,n,t):
[seq(
expand(ToUmbra(gu[j1],[seq(x[i1],i1=1..k),z],[seq(i[i1],i1=1..k),n])),
j1=1..nops(gu))]:
end:


#UAtV(k,x,z,t): The global
#umbral matrix operator, for all consecutive patterns
#of length k+1
#phrased in terms of the indexed variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation
#For example, try:
#UAtV(2,x,z,t);
UAtV:=proc(k,x,z,t) local   gu,i1,mu:
option remember:
gu:=permute(k):
mu:=[seq(U1AtV(gu[i1],x,z,t),i1=1..nops(gu))]:
[seq([seq(mu[i1][j1],i1=1..nops(mu))],j1=1..nops(mu))]:
end:

#UAtVsimp(k,x,z,t): The global
#umbral matrix operator, for all consecutive patterns
#of length k+1
#phrased in terms of the indexed variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation
#For example, try:
#UAtVsimp(2,x,z,t);
UAtVsimp:=proc(k,x,z,t) local   mu,i1,j1:
option remember:
mu:=UAtV(k,x,z,t):
[seq([seq(SimpU(mu[i1][j1]),j1=1..nops(mu[i1]))],i1=1..nops(mu))]:
end:





#Init(k,x,z): The initial generating function for permutations
#of length x. For example, try:
#Init(k,x,z)
Init:=proc(k,x,z) local gu,i:
gu:=z^k*mul(x[i]^i,i=1..k):
[gu$k!]:
end:

#ApplyUmbralMatrix(UM,f,ylist): applies the umbral matrix UM
#to the vector of functions in the variables ylist
#For example, try:
#ApplyUmbralMatrix(UAtV(2,x,z,t),Init(2,x,z),[x[1],x[2],z]);
ApplyUmbralMatrix:=proc(UM,f,ylist) local i1,j1:
if nops({seq(nops(UM[i1]),i1=1..nops(UM))})<>1 then
 RETURN(FAIL):
fi:
if nops(UM[1])<>nops(f) then
 RETURN(FAIL):
fi:

[seq(expand(add(ApplyUmbra(UM[i1][j1],f[j1],ylist),j1=1..nops(UM[i1]))),
i1=1..nops(UM))]:

end:




#AnkMx(n,k,t,x,z): the weight-enumerator of all permutations of
#length n with weight(p) via matrix approach
#product of t[redu([p[i],..,p[i+k-1]])]
#for i=1..n-k+1 by the umbral matrix matrix. For example try:
#AnkMx(4,3,t,x,z);
AnkMx:=proc(n,k,t,x,z) local gu,U,i1:
option remember:

if n<k then
RETURN(FAIL):
fi:

if n=k then 
RETURN(Init(k,x,z)):
fi:

U:=UAtV(k,x,z,t):

gu:=AnkMx(n-1,k,t,x,z):

expand(z^n*ApplyUmbralMatrix(U,gu,[seq(x[i1],i1=1..k),z])):

end:




#AnkM(n,k,t,x,z): the weight-enumerator of all permutations of
#length n  via matrix approach
#product of t[redu([p[i],..,p[i+k-1]])]
#for i=1..n-k+1 by the umbral matrix matrix. For example try:
#AnkM(4,3,t,x,z);
AnkM:=proc(n,k,t) local x,z,i1,gu:
option remember:
gu:=AnkMx(n,k-1,t,x,z):
convert(subs({z=1,seq(x[i1]=1,i1=1..k-1)},gu),`+`):
end:





#Um(x,z,FP): The global
#umbral evolution matrix operator, for enumerating
#permutations avoiding the patterns (of the same length) in FP
#in terms of an umbral matrix and catalytic variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation)
#For example, try:
#Um(x,z,{[1,2,3]});
Um:=proc(x,z,FP) local k,pi,S,t,lu,gu:
option remember:

if FP={} then
 print(`Set of patterns must be non-empty`):
fi:

S:={seq(nops(pi), pi in FP)}:


if nops(S)>1 then
 print(`All patterns must have the same length, not implemented`):
 RETURN(FAIL):
fi:

k:=S[1]:

gu:=UAtV(k-1,x,z,t):



lu:=convert(permute(k),set) minus FP:

gu:=subs({seq(t[op(pi)]=1,pi in lu), seq(t[op(pi)]=0,pi in FP)},gu):



CleanUpU2(gu):

end:


#Umt(x,z,FP,t): The global
#umbral evolution matrix operator, for enumerating
#permutations according to the number
#patterns (of the same length) in FP
#in terms of an umbral matrix and catalytic variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation)
#For example, try:
#Umt(x,z,{[1,2,3]},t);Um
Umt:=proc(x,z,FP,t) local k,pi,S,lu,gu:
option remember:

if FP={} then
 print(`Set of patterns must be non-empty`):
fi:

S:={seq(nops(pi), pi in FP)}:


if nops(S)>1 then
 print(`All patterns must have the same length, not implemented`):
 RETURN(FAIL):
fi:

k:=S[1]:

gu:=UAtV(k-1,x,z,t):



lu:=convert(permute(k),set) minus FP:

gu:=subs({seq(t[op(pi)]=1,pi in lu)},gu):



CleanUpU2(gu):

end:



#AnP(FP,N): The first N terms of the sequence
#enumerating the number of permutations of length n
#avoiding the set of patterns FP
#For example, try:
#AnP({[1,2,3]},10);
AnP:=proc(FP,N) local x,z,gu,pol,i,lu,pi,S,k,i1,n,gu1:
option remember:
if FP={} then
 print(`Set of patterns must be non-empty`):
fi:

S:={seq(nops(pi), pi in FP)}:

if nops(S)>1 then
 print(`the case where all patterns do not have the same length is  not implemented`):
 RETURN(FAIL):
fi:

k:=S[1]:

lu:=Um(x,z,FP):

gu:=[seq(i!,i=1..k-1)]:


pol:=Init(k-1,x,z):



for n from k to N do

pol:=ApplyUmbralMatrix(lu,pol,[seq(x[i1],i1=1..k-1),z]):
 pol:=expand(z^n*pol):
 gu1:=convert(subs({z=1,seq(x[i1]=1,i1=1..k-1)},pol),`+`):
 gu:=[op(gu),gu1]:
od:

gu:

end:


#AnPt(FP,N,t): The first N terms of the sequence
#weight-enumerating the  permutations of length n
#with weight product t[p]^#of occurrences of p for p in FP
#For example, try:
#AnPt({[1,2,3]},10,t);
AnPt:=proc(FP,N,t) local x,z,gu,pol,i,lu,pi,S,k,i1,n,gu1:
option remember:
if FP={} then
 print(`Set of patterns must be non-empty`):
fi:

S:={seq(nops(pi), pi in FP)}:

if nops(S)>1 then
 print(`All patterns must have the same length, not implemented`):
 RETURN(FAIL):
fi:

k:=S[1]:

lu:=Umt(x,z,FP,t):

gu:=[seq(i!,i=1..k-1)]:


pol:=Init(k-1,x,z):



for n from k to N do

pol:=ApplyUmbralMatrix(lu,pol,[seq(x[i1],i1=1..k-1),z]):
 pol:=expand(z^n*pol):
 gu1:=convert(subs({z=1,seq(x[i1]=1,i1=1..k-1)},pol),`+`):
 gu:=[op(gu),expand(gu1)]:
od:

gu:

end:






#UmS(x,z,FP): A simplified version of the global
#umbral evolution matrix operator, for enumerating
#permutations avoiding the patterns (of the same length) in FP
#in terms of an umbral matrix and catalytic variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation)
#For example, try:
#UmS(x,z,{[1,2,3]});
UmS:=proc(x,z,FP) local gu,i1,j1:
 gu:=
Um(x,z,FP):
[seq([seq(SimpU(gu[i1][j1]),j1=1..nops(gu[i1]))],i1=1..nops(gu))]:
end:



#UmNice(x,z,FP,f,g): A nice rendition of 
#an Umbral matrix operator with the input function
#being called f[i] and the output function called g[i]
#For example, try:
#UmNice(x,z,{[1,2,3]},f,g);
UmNice:=proc(x,z,FP,f,g) local gu,i1,j1,lu,mu,lu1:

gu:=UmS(x,z,FP):

for i1 from 1 to nops(gu) do
mu:=0:
 for j1 from 1 to nops(gu[i1]) do
 lu:=gu[i1][j1]:
 mu:=mu+add(lu1[1]*f[j1](op(lu1[2])),lu1 in lu):
 od:
 print(g[i1]=mu):
od:

end:





#ECGF(FP,K): The first K terms of the empirical 
#cluster exponential generating function. For example, try:
#ECGF({[1,3,2]},20);
ECGF:=proc(FP,K) local z,gu,f,i:
gu:=AnP(FP,K):
f:=1/(1+add(gu[i]*z^i/i!,i=1..K)):
f:=taylor(f,z=0,K+1):
if coeff(f,z,0)<>1 or coeff(f,z,1)<>-1 then
 RETURN(FAIL):
fi:

[0,seq(-coeff(f,z,i)*i!,i=2..K)]:
end:



####Stuff from Findrec




###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 or ope=1 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


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


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL and ope<>1 then
       if  SeqFromRec(ope,n,N,[op(1..ORDER,f)],nops(f))=f then
        RETURN(ope):
       fi:
      fi:
 od:
od:
FAIL:

end:





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

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

gu:=Ini:

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

gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:

####End Stuff from Findrec






#AnPlus(N,FP,N1,C,n,E): The first N1 terms of the sequence
#enumerating the number of permutations of length n
#avoiding the set of patterns FP using the first N terms
#and conjecturing a linear recurrence of complexity <=C
#satisfied by
#the coefficients of the cluster generating function
#For example, try:
#AnPlus(20,{[1,2,3]},40,6);
AnPlus:=proc(N,FP,N1,C,n,E) local ope,gu,f,z:
gu:=ECGF(FP,N):
ope:=Findrec(gu,n,E,C):

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

if SeqFromRec(ope,n,E,[op(1..degree(ope,E),gu)],nops(gu))<>gu then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:


gu:=SeqFromRec(ope,n,E,[op(1..degree(ope,E),gu)],N1):


f:=1-z+add(gu[i]*z^i/i!,i=1..N1):

f:=taylor(1/f,z=0,N1+1):
[seq(coeff(f,z,i)*i!,i=1..N1)],ope:


end:




####Brian Nakamura's program#########

Digits:=100:



#-----------------------------------------------------------------------
#  CAV: A Maple package to count permutations that (consecutively) avoid patterns.
#
#  Created: 
#  This version: 
#
#  Written by Brian Nakamura, Rutgers University
#
#-----------------------------------------------------------------------


#print(` CAV: A Maple package to count permutations that (consecutively) avoid patterns. `):
#print(` Written by: Brian Nakamura, Rutgers University `):
#print(` Version: Oct. 11, 2010 `):
print(``):
print(` For a list of the main procedures of CAV, type Help( ); `):
print(` For a list of the supporting procedures of CAV, type Help2( ); `):
print(` For help with a specific main procedure, type Help(procedure_name); `):
print(` For help with a specific supporting procedure, type Help2(procedure_name); `):
print(``):




#-----------------------------------------------------------------------
#
#   HELP PROCEDURES 
#
#-----------------------------------------------------------------------



Help:=proc()

if args=NULL then
  print(` The main procedures are: `):
  print(` CAV(B,N), CAVF(B,N), CAVT(B,N) `):
  print(` CAVt(B,N,t), CAVFt(B,N,t), CAVTt(B,N,t) `):
  print(` RhoApprox(p,N), AsymApprox(p,N) `):
  print(` MapWilfEq(n), WilfEqm(n,N,m), RhoApproxRank(n,N), AsymApproxRank(n,N) `):
  print(` For help with a specific main procedure, type Help(procedure_name); `):
  print(``):

elif nops([args])=1 and op(1,[args])=CAV then
  printf(`CAV(B,N): inputs a set of patterns B (where each pattern is entered as a list) and `):
  printf(`a positive integer N, and outputs a list L=[a1,a2, .., aN], where ai=# of length i `):
  printf(`permutations that consecutively avoid the patterns in B.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`Also, CAVT does the same thing (through a different approach) MUCH FASTER.\n\n`):
  printf(`For an optional verbose/detailed output, do CAV(B,N,true);\n`):
  printf(`For example, try CAV({[1,2,3],[3,2,1]},13); or CAV({[1,2,3,4]},13,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVt then
  printf(`CAVt(B,N,t): inputs a set of patterns B (where each pattern is entered as a list), `):
  printf(`a positive integer N, and a symbolic variable t.\n\n`):
  printf(`It outputs a list L=[b1,b2,.., bN], where bi is the sum of t^(# of consecutive `):
  printf(`occurrences in p of patterns in B) over all length i permutations p.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`Also, CAVTt does the same thing (through a different approach) MUCH FASTER.\n\n`):
  printf(`For an optional verbose/detailed output, do CAVt(B,N,t,true);\n`):
  printf(`For example, try CAVt({[1,2,3],[3,2,1]},13,t); or CAVt({[1,2,3,4]},13,t,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVF then
  printf(`CAVF(B,N): inputs a set of patterns B (where each pattern is entered as a list) and `):
  printf(`a positive integer N, and outputs a list L=[a1,a2, .., aN], where ai=# of length i `):
  printf(`permutations that consecutively avoid the patterns in B.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`Also, this is occasionally faster than CAV, but CAVT does the same thing `):
  printf(`(through a different approach) MUCH FASTER.\n\n`):
  printf(`For an optional verbose/detailed output, do CAVF(B,N,true);\n`):
  printf(`For example, try CAVF({[1,2,3],[3,2,1]},13); or CAVF({[1,2,3,4]},13,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVFt then
  printf(`CAVFt(B,N,t): inputs a set of patterns B (where each pattern is entered as a list), `):
  printf(`a positive integer N, and a symbolic variable t.\n\n`):
  printf(`It outputs a list L=[b1,b2,.., bN], where bi is the sum of t^(# of consecutive `):
  printf(`occurrences in p of patterns in B) over all length i permutations p.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`Also, this is occasionally faster than CAVt, but CAVTt does the same thing `):
  printf(`(through a different approach) MUCH FASTER.\n\n`):
  printf(`For an optional verbose/detailed output, do CAVFt(B,N,t,true);\n`):
  printf(`For example, try CAVFt({[1,2,3],[3,2,1]},13,t); or CAVFt({[1,2,3,4]},13,t,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVT then
  printf(`CAVT(B,N): inputs a set of patterns B (where each pattern is entered as a list) and `):
  printf(`a positive integer N, and outputs a list L=[a1,a2, .., aN], where ai=# of length i `):
  printf(`permutations that consecutively avoid the patterns in B.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`This procedure uses Cluster Tail functional equations and is generally MUCH FASTER `):
  printf(`than CAV and CAVF!\n\n`):
  printf(`For an optional verbose/detailed output, do CAVT(B,N,true);\n`):
  printf(`For example, try CAVT({[1,2,3],[3,2,1]},13); or CAVT({[1,2,3,4]},13,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVTt then
  printf(`CAVTt(B,N,t): inputs a set of patterns B (where each pattern is entered as a list), `):
  printf(`a positive integer N, and a symbolic variable t.\n\n`):
  printf(`It outputs a list L=[b1,b2,.., bN], where bi is the sum of t^(# of consecutive `):
  printf(`occurrences in p of patterns in B) over all length i permutations p.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`This procedure uses Cluster Tail functional equations and is generally MUCH FASTER `):
  printf(`than CAVt and CAVFt!\n\n`):
  printf(`For an optional verbose/detailed output, do CAVTt(B,N,t,true);\n`):
  printf(`For example, try CAVTt({[1,2,3],[3,2,1]},13,t); or CAVTt({[1,2,3,4]},13,t,true);\n\n`):

elif nops([args])=1 and op(1,[args])=RhoApprox then
  printf(`RhoApprox(p,N): inputs a pattern (as a list) and positive integer N>10 and `):
  printf(`computes approximations to rho, the growth constant for the number of permutations `):
  printf(`avoiding pattern p, by counting permutations avoiding p for length N-2, N-1, and N `):
  printf(`permutations.\n\n`):
  printf(`It outputs the approximate rho value followed by the error (difference between consecutive `):
  printf(`rho approximations).\n\n`):
  printf(`NOTE: If a(n)=# of length n permutations avoiding pattern p, then \n a(n) ~ gamma*(rho^n)*n!\n\n`):
  printf(`For an optional verbose/detailed output, do RhoApprox(p,N,true);\n`):
  printf(`For example, try RhoApprox([1,2,3],50); or RhoApprox([1,2,3],50,true);\n\n`):

elif nops([args])=1 and op(1,[args])=AsymApprox then
  printf(`AsymApprox(p,N): inputs a pattern (as a list) and positive integer N>10 and `):
  printf(`computes approximations to rho (the growth constant for the number of permutations `):
  printf(`avoiding pattern p) and gamma (rho's constant coefficient) by counting permutations `):
  printf(`avoiding p for length N-2, N-1, and N permutations.\n\n`):
  printf(`It outputs the approximate rho value, rho's approximate error (difference between `):
  printf(`consecutive rho approximations), the approximate gamma value, and gamma's approximate `):
  printf(`error (difference between consecutive gamma approximations).\n\n`):
  printf(`NOTE: If a(n)=# of length n permutations avoiding pattern p, then \n a(n) ~ gamma*(rho^n)*n!\n\n`):
  printf(`For an optional verbose/detailed output, do AsymApprox(p,N,true);\n`):
  printf(`For example, try AsymApprox([1,2,3],50); or AsymApprox([1,2,3],50,true);\n\n`):

elif nops([args])=1 and op(1,[args])=MapWilfEq then
  printf(`MapWilfEq(n): inputs a postive integer n>=3 and (rigorously!) finds Strong c-Wilf-Equivalences `):
  printf(`via equivalent OverlapMaps (after trivial equivalences of reversals and complementation `):
  printf(`are removed) for single pattern avoidance.\n\n`):
  printf(`It outputs a set of Strong c-Wilf-Equivalence `):
  printf(`classes, where each equivalence class is represented as a set of equivalent patterns.\n\n`):
  printf(`For an optional verbose/detailed output, do MapWilfEq(n,true);\n`):
  printf(`For example, try MapWilfEq(4); or MapWilfEq(4,true);\n\n`):

elif nops([args])=1 and op(1,[args])=WilfEqm then
  printf(`WilfEqm(n,N,m): inputs postive integers n>=3, N>=1, m>=1 `): 
  printf(`and searches for Strong c-Wilf-Equivalences `):
  printf(`via equivalent OverlapMaps (after trivial equivalences of reversals and complementation `):
  printf(`are removed) over sets of m patterns of length n.\n\n`):
  printf(`It outputs a list L=[L1, L2, .. ] and set V, where each Li=[ai, Ai] and `):
  printf(`Ai is the set of length n patterns that is avoided by exactly ai length N permutations.`):
  printf(`Set V is the collection of violators, i.e., the sets of patterns Ai that could not be`):
  printf(`rigorously proven to be Strong c-Wilf-Equivalent by the OverlapMaps approach.\n`):
  printf(`NOTE: if V={}, a rigorous classification was achieved!\n\n`):
  printf(`For an optional verbose/detailed output, do WilfEqm(n,N,m,true);\n`):
  printf(`For example, try WilfEqm(4,15,1); or WilfEqm(4,15,1,true);\n`):
  printf(`Or try WilfEqm(3,10,2); or WilfEqm(3,10,2,true);\n\n`):

elif nops([args])=1 and op(1,[args])=RhoApproxRank then
  printf(`RhoApproxRank(n,N): inputs positive integers n>=3 and N>10 and for each length n pattern p, `):
  printf(`computes approximations to rho, the growth constant for the number of permutations `):
  printf(`avoiding pattern p, by counting permutations avoiding p for length N-2, N-1, and N `):
  printf(`permutations. The rho values are then ranked in increasing order.\n\n`):
  printf(`It outputs a list L=[L1, .., Lk], where each Li is a list with Li=[approx. rho value, `):
  printf(`approx. rho error, set of single patterns that produce these values]. The "approx. rho error" `):
  printf(`is the difference between consecutive rho approximations.\n\n`):
  printf(`NOTE: If a(n)=# of length n permutations avoiding pattern p, then \n a(n) ~ gamma*(rho^n)*n!\n\n`):
  printf(`For an optional verbose/detailed output, do RhoApproxRank(n,N,true);\n`):
  printf(`For example, try RhoApproxRank(4,30); or RhoApproxRank(4,30,true);\n\n`):

elif nops([args])=1 and op(1,[args])=AsymApproxRank then
  printf(`AsymApproxRank(n,N): inputs positive integers n>=3 and N>10 and for each length n pattern p, `):
  printf(`computes approximations to rho (the growth constant for the number of permutations `):
  printf(`avoiding pattern p) and gamma (rho's constant coefficient) by counting permutations `):
  printf(`avoiding p for length N-2, N-1, and N permutations. The values are ranked in increasing `):
  printf(`order for rho.\n\n`):
  printf(`It outputs a list L=[L1, .., Lk], where each Li is a list with Li=[approx. rho value, `):
  printf(`approx. rho error, approx. gamma value, approx. gamma error, set of single patterns that `):
  printf(`produce these values]. The approx. rho (gamma) "error" `):
  printf(`is the difference between consecutive rho (gamma) approximations).\n\n`):
  printf(`NOTE: If a(n)=# of length n permutations avoiding pattern p, then \n a(n) ~ gamma*(rho^n)*n!\n\n`):
  printf(`For an optional verbose/detailed output, do AsymApproxRank(n,N,true);\n`):
  printf(`For example, try AsymApproxRank(4,30); or AsymApproxRank(4,30,true);\n\n`):

else
  ERROR(`Invalid Help parameters!`):
fi:

end:





Help2:=proc()

if args=NULL then
  print(` The primary supporting procedures are: `):
  print(` SCHEME(k,B,x,y,t), MakeTailFE(B,k,z,t), ScToTailSc(Sc,z) `):
  print(` CAVOID(Sc,N,k), CAVOIDtail(TSc,N,k) `):
  print(` OverlapMaps(p1,p2), CompactSets(S), RemoveTrivEq(S) `):
  print(` For help with a specific supporting procedure, type Help2(procedure_name); `):
  print(``):

elif nops([args])=1 and op(1,[args])=SCHEME then
  printf(`SCHEME(k,B,x,y,t): inputs a set of patterns B (where each pattern is entered as a list) `):
  printf(`and symbolic variables k,x,y,t, and outputs the basic scheme used to count permutations `):
  printf(`avoiding patterns in B. See the accompanying paper on more information regarding the scheme.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`If symbolic variable t is inputted, the scheme generated will be to count occurrences of `):
  printf(`patterns in B. If 0 is entered in place of t, the scheme will be for counting permutations that `):
  printf(`completely avoid B.\n`):
  printf(`This scheme can be used directly with CAVOID. To use CAVOIDtail, you must convert this scheme `):
  printf(`into a tail scheme via ScToTailSc first.\n\n`):
  printf(`For example, try SCHEME(k,{[1,2,3,4]},x,y,t); or SCHEME(k,{[1,2,3],[3,2,1]},x,y,0); \n\n`):

elif nops([args])=1 and op(1,[args])=MakeTailFE then
  printf(`MakeTailFE(B,k,z,t): inputs a set of patterns B (where each pattern is entered as a list) `):
  printf(`and symbolic variables k,z,t, and outputs the tail scheme (Cluster Tail GF functional `):
  printf(`equation) for counting permutations avoiding patterns in B. See the accompanying paper `):
  printf(`on more information regarding the tail scheme.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`If symbolic variable t is inputted, the tail scheme generated will be to count occurrences of `):
  printf(`patterns in B. If 0 is entered in place of t, the tail scheme will be for counting permutations `):
  printf(`that completely avoid B.\n`):
  printf(`This tail scheme can be used directly with CAVOIDtail. To use CAVOID, you must generate a basic `):
  printf(`scheme using SCHEME. However, the tail scheme approach is faster.\n\n`):
  printf(`For example, try MakeTailFE({[1,2,3]},k,z,t); or MakeTailFE({[1,2,3],[3,2,1]},k,z,0); \n\n`):

elif nops([args])=1 and op(1,[args])=ScToTailSc then
  printf(`ScToTailSc(Sc,z): inputs a basic scheme Sc (generated from SCHEME) and symbolic variable z, `):
  printf(`and outputs the tail scheme (Cluster Tail GF functional equation) for counting permutations `):
  printf(`avoiding patterns in B. See the accompanying paper on more information regarding the `):
  printf(`tail scheme.\n\n`):
  printf(`NOTE: The input of this procedure must be generated from SCHEME.\n`):
  printf(`The tail scheme generated can be used for CAVOIDtail, which is faster than CAVOID.\n\n`):
  printf(`For example, try ScToTailSc(Sc1,z); where Sc1 is a valid output out of SCHEME.\n\n`):

elif nops([args])=1 and op(1,[args])=CAVOID then
  printf(`CAVOID(Sc,N,k): inputs a basic scheme Sc (generated by SCHEME), a positive integer N, and `):
  printf(`symbolic variable k (that was used in the scheme Sc). It outputs the # of length N permutations `):
  printf(`avoiding a set of patterns (which is essentially encoded in the scheme). IF the scheme was `):
  printf(`created to count occurrences of patterns in B, then the output is b0+b1*t+b2*t^2+... , where `):
  printf(`bi=# of length N permutations with EXACTLY i occurrences of patterns encoded in Sc.\n\n`):
  printf(`NOTES: The input of this procedure must be generated from SCHEME.\n`):
  printf(`The procedure CAVOIDtail, which uses a tail scheme input, is generally much faster.\n\n`):
  printf(`For example, try CAVOID(Sc1,10,k); where Sc1 is a valid output out of SCHEME.\n\n`):

elif nops([args])=1 and op(1,[args])=CAVOIDtail then
  printf(`CAVOIDtail(TSc,N,k): inputs a tail scheme TSc (generated by either MakeTailFE directly or `):
  printf(`converted from a basic scheme by ScToTailSc), a positive integer N, and symbolic variable k `):
  printf(`(that was used in the scheme TSc). It outputs the # of length N permutations avoiding a set `):
  printf(`of patterns (which is essentially encoded in the scheme). IF the scheme was created to count `):
  printf(`occurrences of patterns in B, then the output is b0+b1*t+b2*t^2+... , where `):
  printf(`bi=# of length N permutations with EXACTLY i occurrences of patterns encoded in TSc.\n\n`):
  printf(`NOTES: The input of this procedure must be generated from MakeTailFE or converted from a basic `):
  printf(`scheme using ScToTailSc.\n`):
  printf(`This procedure is generally much faster than CAVOID.\n\n`):
  printf(`For example, try CAVOIDtail(TSc1,10,k); where TSc1 is a valid output of MakeTailFE.\n\n`):

elif nops([args])=1 and op(1,[args])=OverlapMaps then
  printf(`OverlapMaps(p1,p2): inputs two patterns p1 and p2, and outputs a set S={M1,M2, ..} of all `):
  printf(`possible overlaps where a tail of p1 overlaps with a head of p2. Each Mi={[b1,a1],[b2,a2], ..} `):
  printf(`is a mapping for a possible overlap, where p1's tail is a1,a2,.. and p2's head is b1,b2,... \n`):
  printf(`See the accompanying article for more information.\n\n`):
  printf(`For example, try OverlapMaps([1,3,2,4],[1,3,2,4]); or OverlapMaps([1,2,3,4],[2,4,1,3]);\n\n`):

elif nops([args])=1 and op(1,[args])=CompactSets then
  printf(`CompactSets(S): inputs a set S of sets of patterns (where each pattern is entered as a list), `):
  printf(`and outputs S modulo all trivial c-Wilf-Equivalent sets (reversals and complementations).\n\n`):
  printf(`For example, try CompactSets(S); where S is the set of all pairs of length 3 permutations.\n`):
  printf(`(i.e., S=choose(convert(permute(3),set),2); )\n\n`):

elif nops([args])=1 and op(1,[args])=RemoveTrivEq then
  printf(`RemoveTrivEq(S): inputs a set of single patterns S (where each pattern is entered as a list), `):
  printf(`and outputs S modulo all the trivial c-Wilf-Equivalences of reversals and complementations.\n\n`):
  printf(`For example, try RemoveTrivEq(S3); where S3=convert(permute(3),set);\n\n`):

else
  ERROR(`Invalid Help parameters!`):
fi:

end:





#-----------------------------------------------------------------------
#
#   MAIN PROCEDURES 
#
#-----------------------------------------------------------------------


# CAV(B,N): Computes the consecutive avoidance sequence through the general "scheme"
#   approach. NOTE: this is sometimes slower than CAVF and usually slower than CAVT
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the number of length k permutations
#   that avoid (consecutively) the set of patterns B
#
# Ex: CAV({[1,2,3],[3,2,1]}, 10);
#
CAV:=proc(B,N)
local S,i,verb,L:
global k,x,y:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=3 then
  if type(args[3],boolean) then
    verb:=args[3]:
  else
    ERROR(`Optional 3rd argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,0):

L:=[seq(CAVOID(S,i,k),i=1..N)]:


if verb then
  CAVprint(B,S,L):
else
  RETURN(L):
fi:

end:




# CAVt(B,N,t): Analogous to CAV except that the number of occurrences of bad patterns
#   are marked. NOTE: this is sometimes slower than CAVFt and usually slower than CAVTt
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#   t - variable
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the weight enumerated sum over length
#   k permutations, where the weight of permutation p is t^(# occurrences of bad patterns)
#
# Ex: CAVt({[1,2,3],[3,2,1]}, 10, t);
#
CAVt:=proc(B,N,t)
local S,i,verb,L:
global k,x,y:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=4 then
  if type(args[4],boolean) then
    verb:=args[4]:
  else
    ERROR(`Optional 4th argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,t):

L:=expand([seq(CAVOID(S,i,k),i=1..N)]):

if verb then
  CAVtprint(B,S,L,t):
else
  RETURN(L):
fi:

end:




# CAVF(B,N): Computes the consecutive avoidance sequence through the general "scheme"
#   approach, but by defining additional functions in the "scheme" to better utilize "option
#   remember". NOTE: this is sometimes faster than CAV and usually slower than CAVT
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the number of length k permutations
#   that avoid (consecutively) the set of patterns B
#
# Ex: CAVF({[1,2,3],[3,2,1]}, 10);
#
CAVF:=proc(B,N)
local S,i,verb,L:
global k,x,y:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=3 then
  if type(args[3],boolean) then
    verb:=args[3]:
  else
    ERROR(`Optional 3rd argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,0):

L:=[seq(CAVOIDF(S,i,k,x,y),i=1..N)]:

if verb then
  CAVFprint(B,S,L):
else
  RETURN(L):
fi:

end:





# CAVFt(B,N,t): Analogous to CAVF except that the number of occurrences of bad patterns
#   are marked. NOTE: this is sometimes faster than CAVt and usually slower than CAVTt
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#   t - variable
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the weight enumerated sum over length
#   k permutations, where the weight of permutation p is t^(# occurrences of bad patterns)
#
# Ex: CAVFt({[1,2,3],[3,2,1]}, 10, t);
#
CAVFt:=proc(B,N,t)
local S,i,verb,L:
global k,x,y:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=4 then
  if type(args[4],boolean) then
    verb:=args[4]:
  else
    ERROR(`Optional 4th argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,t):

L:=expand([seq(CAVOIDF(S,i,k,x,y),i=1..N)]):

if verb then
  CAVFtprint(B,S,L,t):
else
  RETURN(L):
fi:

end:




# CAVT(B,N): Computes the consecutive avoidance sequence through the general "scheme"
#   approach combined with the cluster tail generating function approach. 
#   NOTE: this is generally faster than CAV and CAVF
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the number of length k permutations
#   that avoid (consecutively) the set of patterns B
#
# Ex: CAVT({[1,2,3],[3,2,1]}, 10);
#
CAVT:=proc(B,N)
local S,i,FE,verb,L:
global k,x,y,z:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=3 then
  if type(args[3],boolean) then
    verb:=args[3]:
  else
    ERROR(`Optional 3rd argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,0):

FE:=ScToTailSc(S,z):

L:=[seq(CAVOIDtail(FE,i,k),i=1..N)]:

if verb then
  CAVTprint(B,S,FE,L):
else
  RETURN(L):
fi:

end:







# CAVTt(B,N,t): Analogous to CAVT except that the number of occurrences of bad patterns
#   are marked. NOTE: this is generally faster than CAV and CAVF
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#   t - variable
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the weight enumerated sum over length
#   k permutations, where the weight of permutation p is t^(# occurrences of bad patterns)
#
# Ex: CAVTt({[1,2,3],[3,2,1]}, 10, t);
#
CAVTt:=proc(B,N,t)
local S,i,FE,verb,L:
global k,x,y,z:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=4 then
  if type(args[4],boolean) then
    verb:=args[4]:
  else
    ERROR(`Optional 4th argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,t):

FE:=ScToTailSc(S,z):

L:=expand([seq(CAVOIDtail(FE,i,k),i=1..N)]):

if verb then
  CAVTtprint(B,S,FE,L,t):
else
  RETURN(L):
fi:

end:





# RhoApprox(p,N): inputs a SINGLE pattern p, and a pos. integer N
# and outputs the Elizalde constant rho (limit of t[n]/t[n-1]/n
# as n goes to infinity followed by the error.
RhoApprox:=proc(p,N) 
local S,FE,i,gu,mu,verb:
global k,x,y,z:

if N<=10 then
  ERROR(`N must be greater than 10`):
fi:

verb:=false:

if nops([args])>2 then
  if type(op(3,[args]),boolean) then
    verb:=op(3,[args]):
  else
    ERROR(`Optional 3rd input (verbose) must be a Boolean. Instead, received %1.`, op(3,[args])):
  fi:
fi:

S:=SCHEME(k,{p},x,y,0);
FE:=ScToTailSc(S,z):
gu:=[CAVOIDtail(FE,N-2,k),CAVOIDtail(FE,N-1,k),CAVOIDtail(FE,N,k)]:
mu:=[seq(evalf(gu[i+1]/gu[i]/(N-2+i)),i=1..2)]:

if verb then
  RhoApproxprint([mu[2],mu[2]-mu[1]],N,p):
else
  RETURN(mu[2],mu[2]-mu[1]):
fi:

end:






# AsymApprox(p,N): inputs a SINGLE pattern p, and a pos. integer N
# and outputs the Elizalde constant rho (limit of t[n]/t[n-1]/n
# as n goes to infinity, an approximation error for rho, the Elizalde
# coefficient gamma, and an approximation error for gamma
AsymApprox:=proc(p,N)
local S,FE,i,gu,mu,cu,verb:
global k,x,y,z:

if N<=10 then
  ERROR(`N must be greater than 10`):
fi:

verb:=false:

if nops([args])>2 then
  if type(op(3,[args]),boolean) then
    verb:=op(3,[args]):
  else
    ERROR(`Optional 3rd input (verbose) must be a Boolean. Instead, received %1.`, op(3,[args])):
  fi:
fi:

S:=SCHEME(k,{p},x,y,0);
FE:=ScToTailSc(S,z):
gu:=[CAVOIDtail(FE,N-2,k),CAVOIDtail(FE,N-1,k),CAVOIDtail(FE,N,k)]:
mu:=[seq(evalf(gu[i+1]/gu[i]/(N-2+i)),i=1..2)]:
cu:=[gu[2]/(mu[1]^(N-1))/(N-1)!, gu[3]/(mu[2]^N)/N!]:

if verb then
  AsymApproxprint([mu[2],mu[2]-mu[1],cu[2],cu[2]-cu[1]],N,p):
else
  RETURN(mu[2],mu[2]-mu[1],cu[2],cu[2]-cu[1]):
fi:

end:




MapWilfEq:=proc(n)
local Sn,Sn1,p,MapSet,M1,M2,L,eqset,i,verb:

verb:=false:

if nops([args])>1 then
  if type(op(2,[args]),boolean) then
    verb:=op(2,[args]):
  else
    ERROR(`Optional 2nd input (verbose) must be a Boolean. Instead, received %1.`, op(2,[args])):
  fi:
fi:


Sn:=convert(permute(n),set):
Sn1:=RemoveTrivEq(Sn):

MapSet:={}:

for p in Sn1 do
  MapSet:=MapSet union {[OverlapMaps(p,p),p]}:
od:

M1:={}:
L:=[]:
for i from 1 to nops(MapSet) do
  M2:=MapSet[i][1]:

  if M1=M2 then
    eqset:=L[nops(L)] union {MapSet[i][2]}:
    L:=[op(1..nops(L)-1,L), eqset]:
  else
    L:=[op(L), {MapSet[i][2]}]:
  fi:

  M1:=M2:

od:

# Check for verbose output
if verb then
  MapWEprint(n,convert(L,set)):
else
  RETURN(convert(L,set)):
fi:

end:






WilfEqm:=proc(n,N,m)
local Sn,eqSet,eqClass,p,Sc,FE,L,L2,a1,a2,V,i,Li,Lprev,Sm,verb:
global k,x,y,z:


verb:=false:

if nops([args])>3 then
  if type(op(4,[args]),boolean) then
    verb:=op(4,[args]):
  else
    ERROR(`Optional 4th input (verbose) must be a Boolean. Instead, received %1.`, op(4,[args])):
  fi:
fi:

L:={}:

if m=1 then 
  eqSet:=MapWilfEq(n):

  for eqClass in eqSet do
    p:=eqClass[1]:
    Sc:=SCHEME(k,{p},x,y,0):
    FE:=ScToTailSc(Sc,z):
    L:=L union {[CAVOIDtail(FE,N,k), eqClass]}:
  od:

  L:=convert(L,list):

  L2:={ [L[1][1], {L[1][2]}]}:

  for i from 2 to nops(L) do
    Li:=L[i]:
    Lprev:=L2[nops(L2)]:

    if Li[1]=Lprev[1] then
      L2:=L2 minus {Lprev}:
      L2:=L2 union {[Lprev[1], Lprev[2] union {Li[2]}]}:
    else
      L2:=L2 union {[Li[1], {Li[2]}]}:
    fi:
  od:

  V:={}:

  for Li in L2 do
    if nops(Li[2])>1 then
      V:=V union {Li[2]}:
    fi:
  od:

else
  Sm:=choose(convert(permute(n),set),m):
  Sm:=CompactSets(Sm):
  eqSet:={seq({p}, p in Sm)}:

  for eqClass in eqSet do
    p:=eqClass[1]:
    Sc:=SCHEME(k,p,x,y,0):
    FE:=ScToTailSc(Sc,z):
    L:=L union {[CAVOIDtail(FE,N,k), eqClass]}:
  od:

  L:=convert(L,list):

  L2:={ [L[1][1], L[1][2]]}:

  for i from 2 to nops(L) do
    Li:=L[i]:
    Lprev:=L2[nops(L2)]:

    if Li[1]=Lprev[1] then
      L2:=L2 minus {Lprev}:
      L2:=L2 union {[Lprev[1], Lprev[2] union Li[2]]}:
    else
      L2:=L2 union {[Li[1], Li[2]]}:
    fi:
  od:

  V:={}:

  for Li in L2 do
    if nops(Li[2])>1 then
      if SameScheme2(k,Li[2],x,y)=false then
        V:=V union {Li[2]}:
      fi:
    fi:
  od:

fi:


if verb then
  WilfEqmprint(L,V,n,N,m):
else
  RETURN(L,V):
fi:


end:





RhoApproxRank:=proc(n,N)
local eqSet,S,FE,i,gu,mu,p,eqClass,L,verb:
global k,x,y,z:


verb:=false:

if nops([args])>2 then
  if type(op(3,[args]),boolean) then
    verb:=op(3,[args]):
  else
    ERROR(`Optional 3rd input (verbose) must be a Boolean. Instead, received %1.`, op(3,[args])):
  fi:
fi:


eqSet:=MapWilfEq(n):

if N<=10 then
  ERROR(`N must be greater than 10`):
fi:

L:={}:

for eqClass in eqSet do
  p:=eqClass[1]:

  S:=SCHEME(k,{p},x,y,0);
  FE:=ScToTailSc(S,z):
  gu:=[CAVOIDtail(FE,N-2,k),CAVOIDtail(FE,N-1,k),CAVOIDtail(FE,N,k)]:
  mu:=[seq(evalf(gu[i+1]/gu[i]/(N-2+i)),i=1..2)]:
  L:=L union {[mu[2],mu[2]-mu[1],eqClass]}:

od:

if verb then
  RhoApproxRankprint(convert(L,list),N):
else
  RETURN(convert(L,list)):
fi:

end:








AsymApproxRank:=proc(n,N)
local eqSet,S,FE,i,gu,mu,p,eqClass,L,cu,verb:
global k,x,y,z:


verb:=false:

if nops([args])>2 then
  if type(op(3,[args]),boolean) then
    verb:=op(3,[args]):
  else
    ERROR(`Optional 3rd input (verbose) must be a Boolean. Instead, received %1.`, op(3,[args])):
  fi:
fi:

eqSet:=MapWilfEq(n):

if N<=10 then
  ERROR(`N must be greater than 10`):
fi:

L:={}:

for eqClass in eqSet do
  p:=eqClass[1]:

  S:=SCHEME(k,{p},x,y,0);
  FE:=ScToTailSc(S,z):
  gu:=[CAVOIDtail(FE,N-2,k),CAVOIDtail(FE,N-1,k),CAVOIDtail(FE,N,k)]:
  mu:=[seq(evalf(gu[i+1]/gu[i]/(N-2+i)),i=1..2)]:
  cu:=[gu[2]/(mu[1]^(N-1))/(N-1)!, gu[3]/(mu[2]^N)/N!]:
  L:=L union {[mu[2],mu[2]-mu[1],cu[2],cu[2]-cu[1],eqClass]}:

od:

if verb then
  AsymApproxRankprint(convert(L,list),N):
else
  RETURN(convert(L,list)):
fi:


end:






MakeTailFE:=proc(B,k,z,t)
local S,i,FE:
global x,y:

S:=SCHEME(k,B,x,y,t):

FE:=ScToTailSc(S,z):

end:





TailFE2:=proc(TSc,FL)
local S,receq,LHS,RHS:

S:={}:

for receq in TSc do
  LHS:=receq[2]:
  RHS:=receq[3]:
  S:=S union {[receq[1], FL[LHS[1]][LHS[2], LHS[3]] = TailFE2Exp(RHS, FL)]}:
od:

S:

end:




TailFE2Exp:=proc(RHS,FL)
local i:

if type(RHS,list) then
  if type(RHS[1],posint) then
    RETURN(FL[RHS[1]][RHS[2], RHS[3]]):
  elif RHS[1]=`+` then
    RETURN(add(TailFE2Exp(RHS[i],FL), i=2..nops(RHS))):
  elif RHS[1]=`*` then
    RETURN(mul(TailFE2Exp(RHS[i],FL), i=2..nops(RHS))):
  else
    ERROR(`Invalid Tail FE!`):
  fi:
else
  RETURN(RHS):
fi:  

end:






TailFEprint:=proc(TSc,FL)
local TCSc,TFE,n,i:

n:=NumFuncInScheme(TSc):

if n<>nops(FL) then
  ERROR(`The number of function names in list must equal the number of fuctions in Scheme.`):
fi:

TCSc:=ConvertScheme2(TSc,n):

printf(`We have the following Cluster Tail Generating Function(s):\n\n`):

for i from 1 to n do
  TFE:=TailFE2(TSc,FL):

  printf(`%a\n\n`, TFE[2][2]):
  printf(`with initial conditions: %a, if %a\n`, TFE[3][2], op(TFE[3][1])):
  printf(`and %a, if %a\n\n`, TFE[1][2], op(TFE[1][1])):
od:

printf(`These are encoded as the following scheme:\n\n%a`,TSc):
printf(`\n\n`):

end:





# Checks whether every pattern in S1 has the same scheme
SameScheme2:=proc(k,S1,x,y)
local Sc1,Sc2,i:

if nops(S1)=1 then
  RETURN(true):
fi:


Sc1:=SCHEME(k,S1[1],x,y,0):

for i from 2 to nops(S1) do
  Sc2:=SCHEME(k,S1[i],x,y,0):
  if evalb(Sc1<>Sc2) then
    RETURN(false):
  fi:
od:

RETURN(true):

end:





#-----------------------------------------------------------------------
#
#   CLUSTER TAIL SCHEME EXECUTION PROCEDURES 
#
#-----------------------------------------------------------------------


CAVOIDtail:=proc(Sc,N,k)
local j,CSc,n:
option remember:

if N<3 then
  RETURN(N!):
fi:

n:=NumFuncInScheme(Sc):

CSc:=ConvertScheme2(Sc,n):

N*CAVOIDtail(Sc,N-1,k)+
add(binomial(N,j)*ComputeTailCN(CSc,j,n,k)*CAVOIDtail(Sc,N-j,k),j=1..N):


end:






ComputeTailCN:=proc(CSc,N,n,k) local i,j,zvars,zcond:
option remember:

zvars:={}:
zcond:={}:

for i from 1 to n do
  zvars:=zvars union convert(CSc[i][1][2][3],set):
od:

zcond:={seq(zvars[i]=1, i=1..nops(zvars))}:

subs(zcond, simplify(expand(add(ComputeTailFiN(CSc,i,N,k),i=1..n)))):

end:


ComputeTailFiN:=proc(CSc,i,N,k)
local Ri,eq,RHS,total,j:
option remember:

Ri:=CSc[i]:

for eq in Ri do
  if evalb(subs(k=N,op(eq[1]))) then
    RHS:=eq[3]:

    if not type(RHS,list) then
      RETURN(subs(k=N,RHS)):
    else
      if member(RHS[1],{0,`+`,`*`}) then
        RETURN(EvalTailExpression(CSc,RHS,N,k)):
      else
        RETURN(ComputeTailFiN(CSc,RHS[1],subs(k=N,RHS[2]),k)):
      fi:
    fi:
  fi:

od:


end:



EvalTailExpression:=proc(CSc,RHS,N,k)
local j,i,zvars,zcond:
option remember:

if not type(RHS,list) then
  RETURN(subs(k=N,RHS)):
elif evalb(RHS[1]=`+`) then
  RETURN(add(EvalTailExpression(CSc,RHS[j],N,k),j=2..nops(RHS))):
elif evalb(RHS[1]=`*`) then
  RETURN(mul(EvalTailExpression(CSc,RHS[j],N,k),j=2..nops(RHS))):
else

  zcond:={}:

  i:=RHS[1]:

  zvars:=CSc[i][1][2][3]:

  zcond:={seq(zvars[i]=RHS[3][i], i=1..nops(zvars))}:

  RETURN(subs(zcond, simplify(expand(ComputeTailFiN(CSc,RHS[1],subs(k=N,RHS[2]),k))))):

fi:


end:




#-----------------------------------------------------------------------
#
#   (GENERAL) SCHEME EXECUTION PROCEDURES 
#
#-----------------------------------------------------------------------


CAVOID:=proc(Sc,N,k)
local j,CSc,n:
option remember:

if N<3 then
  RETURN(N!):
fi:

n:=NumFuncInScheme(Sc):

CSc:=ConvertScheme2(Sc,n):

N*CAVOID(Sc,N-1,k)+
add(binomial(N,j)*ComputeCN(CSc,j,n,k)*CAVOID(Sc,N-j,k),j=1..N):


end:





CAVOIDF:=proc(Sc,N,k,x,y)
local j,CESc,n1,n2,ESc:
option remember:

if N<3 then
  RETURN(N!):
fi:

n1:=NumFuncInScheme(Sc):

ESc:=ExpandScheme(Sc,k,x,y):

n2:=NumFuncInScheme(ESc):

CESc:=ConvertScheme2(ESc,n2):

CAVOIDF2(CESc,N,k,n1):

end:




CAVOIDF2:=proc(CESc,N,k,n1)
local j:
option remember:

if N<3 then
  RETURN(N!):
fi:

N*CAVOIDF2(CESc,N-1,k,n1)+
add(binomial(N,j)*ComputeCN(CESc,j,n1,k)*CAVOIDF2(CESc,N-j,k,n1),j=1..N):

end:



# computes the sum of weights of all length k clusters
ComputeCN:=proc(CSc,N,n,k) local i:
option remember:
add(ComputeFiN(CSc,i,N,k),i=1..n):
end:


# computes the sum of weights of all length k clusters that end in the i-th pattern
ComputeFiN:=proc(CSc,i,N,k) local xn,S,xvect:
option remember:
xn:=nops(CSc[i][2][2][3]):
S:={op(choose(N,xn))}:
add(ComputeFiNx(CSc,i,N,xvect,k), xvect in S):
end:




# computes the sum of weights of all length k clusters that end in the i-th pattern
# given by xvect
ComputeFiNx:=proc(CSc,i,N,xvect,k)
local Ri,eq,RHS,total,j:
option remember:

Ri:=CSc[i]:

for eq in Ri do
  if evalb(subs(k=N,op(eq[1]))) then
    RHS:=eq[3]:

    if not type(RHS,list) then
      RETURN(RHS):
    else
      if member(RHS[1],{0,`+`,`*`}) then
        RETURN(EvalExpression(CSc,RHS,N,xvect,k,eq[2][3])):
      else
        RETURN(ComputeFiNx(CSc,RHS[1],subs(k=N,RHS[2]),RHS[3],k)):
      fi:
    fi:
  fi:

od:

end:



# evalutaes a summation in a scheme
ComputeSummation:=proc(CSc,sum1,N,xvect,k,xlist)
local S,N2,Sfilt,sj:
option remember:

N2:=subs(k=N,sum1[2][2]):
S:={op(choose(N2,nops(sum1[2][3])))}:

Sfilt:={}:

for sj in S do
  if CheckCond2(xvect,sj,sum1[3],xlist,sum1[2][3]) then
    Sfilt:=Sfilt union {sj}:
  fi:
od:

add(ComputeFiNx(CSc,sum1[2][1],N2,sj,k), sj in Sfilt):

end:



# evaluates a general expression in a scheme
EvalExpression:=proc(CSc,RHS,N,xvect,k,xlist)
local j,cond:
option remember:

if not type(RHS,list) then
  RETURN(RHS):
elif RHS[1]=0 then
  RETURN(ComputeSummation(CSc,RHS,N,xvect,k,xlist)):
elif evalb(RHS[1]=`+`) then
  RETURN(add(EvalExpression(CSc,RHS[j],N,xvect,k,xlist),j=2..nops(RHS))):
elif evalb(RHS[1]=`*`) then
  RETURN(mul(EvalExpression(CSc,RHS[j],N,xvect,k,xlist),j=2..nops(RHS))):
else
  cond:={seq(xlist[j]=xvect[j],j=1..nops(xvect))}:
  RETURN(ComputeFiNx(CSc,RHS[1],subs(k=N,RHS[2]),subs(cond,RHS[3]),k)):

fi:

end:



# Converts a Scheme to a Converted Scheme so that it can quickly find the
# necessary recurrence for the i-th function (1<=i<=n)
ConvertScheme2:=proc(Sc,n)
local i,CSc:
option remember:

for i from 1 to n do
  CSc[i]:={}:
od:

for i from 1 to nops(Sc) do
  CSc[Sc[i][2][1]]:=CSc[Sc[i][2][1]] union {Sc[i]}:
od:

CSc:

end:



CheckCond2:=proc(xvect,ivect,cond,xlist,ilist)
local c1,S,j:
option remember:

S:={seq(xlist[j]=xvect[j],j=1..nops(xvect)), seq(ilist[j]=ivect[j],j=1..nops(ivect))}:

for c1 in cond do
  if not evalb(subs(S,c1)) then
    RETURN(false):
  fi:
od:

RETURN(true):

end:




#-----------------------------------------------------------------------
#
#   SCHEME GENERATION PROCEDURES 
#
#-----------------------------------------------------------------------


SCHEME:=proc(k,B,x,y,t)
local B1,Sc,i,j,r,xvect,EqFi,RecurFi:

CheckSetB(B):

# Bad pattern set B modulo redundancies
B1:=[op(RemoveRedundant(B))]:

Sc:={}:

for i from 1 to nops(B1) do
  r:=nops(B1[i]):
  xvect:=[seq(x[j],j=1..r)]:

# MODDED
  EqFi:={[{k<r}, [i,k,xvect], 0], [{k=r}, [i,k,xvect], t-1]}:

  RecurFi:=MakeRecurrenceFi(k,B1,x,y,i,t):

  # Need if statement to either add RecurFi to EqFi or not
  # Maybe this is not necessary?
  # If this fails, should all of EqFi be scrapped?


# MODDED
  if RecurFi<>[`*`, t-1, []] then
    EqFi:=EqFi union {[{k>r}, [i,k,xvect], RecurFi]}:
    Sc:=Sc union EqFi:
  fi:

od:

Sc:

end:




MakeRecurrenceFi:=proc(k,B1,x,y,i,t)
local RHS,RHSterms,j:

RHS:=[]:

for j from 1 to nops(B1) do
  RHSterms:=FindSums2(k,B1,x,y,j,i):
  RHS:=[op(RHS), op(RHSterms)]:
od:

if nops(RHS)>1 then
  RHS:=[`+`, op(RHS)]:
elif nops(RHS)=1 then
  RHS:=op(RHS):
fi:

# MODDED
[`*`, t-1, RHS]:

end:




FindSums2:=proc(k,B1,x,y,j,i)
local M,ir,jr,S,M1,m1,k2,b,chopped,t,cond:

M:=OverlapMaps(B1[j],B1[i]):

if M={} then
  RETURN({}):
fi:

# lengths of bad words
ir:=nops(B1[i]):
jr:=nops(B1[j]):

S:={}:

for M1 in M do
  k2:=k-ir+nops(M1):

  chopped:={seq(b,b=1..ir)} minus {seq(m1[1],m1 in M1)}:

  cond:={ seq(y[t]<y[t+1], t=1..jr-1) }:

  for m1 in M1 do
    # m1[1] -> m1[2]
    cond:=cond union {x[m1[1]]-nops(chopped intersect {seq(b,b=1..m1[1])})=y[m1[2]]}:
  od:

  S:=S union {[0, [j, k2, [seq(y[t], t=1..jr)]], cond]}:
od:

S:

end:





# OverlapMaps(p1,p2)
#
# Inputs: p1,p2 - permutation patterns
# Output: all "maps" where the tail/end of p1 can overlap with
#           the head/beginning of p2 in a larger permutation
#
# Comments:
#   The output is a set of sets {M1, M2, .., Mk} with each Mi
#   a set of lists. For example, if M1={[a1,b1], [a2,b2], ..},
#   then the a1-th term of p1 coincides with the b1-th term of p2,
#   the a2-th term of p1 coincides with the b2-th term of p2, etc.
#
# Example:
#   OverlapMaps([1,2,3], [1,2,3]) 
#     outputs {{[1, 3]}, {[1, 2], [2, 3]}}
# since it can have an overlap of length 1 (shown by {[1,3]})
#        1 2 3
#    1 2 3
# or an overlap of length 2 (shown by {[1, 2], [2, 3]})
#      1 2 3
#    1 2 3
#
OverlapMaps:=proc(p1,p2)
local i,Tails,Heads,OL,olap,S,S1,k1,k2,m:

k1:=nops(p1):
k2:=nops(p2):

Tails:={}:
Heads:={}:

for i from 1 to k1-1 do
  Tails:={op(Tails), ScalePerm([op(i+1..k1,p1)])}:
od:

for i from 1 to k2-1 do
  Heads:={op(Heads), ScalePerm([op(1..i,p2)])}:
od:

OL:=Heads intersect Tails:

S:={}:

for olap in OL do

  # length of the overlap (in terms of entries in p2)
  m:=nops(olap):

  S1:={}:

  for i from 1 to m do
    S1:={op(S1), [p2[i],p1[k1-m+i]]}:
  od:

  S:={op(S), S1}:

od:

S:

end:







# ScalePerm(p)
#
# Inputs: p - a subsequence of a permutation 
# Output: the canonical reduction of p
#
# Comments: 
#   If p contains k terms, the output will be the permutation
#   of {1,2,..,k} that follows the same respective ordering 
#   as in p
#
# Example:
#   ScalePerm([1,8,4,6]) returns [1,4,2,3]
ScalePerm:=proc(p)
local n,i,k,psorted,L:

n:=nops(p):

psorted:={op(p)}:
psorted:=[op(psorted)]:

L:=[]:

for i from 1 to n do
  k:=1:
  
  while k<=n do
    if psorted[k]=p[i] then
      L:=[op(L), k]:
      k:=n+1:
    else
      k:=k+1:
    fi:
  od:
  
od:

L:

end:



# RemoveRedundant(B)
#
# Inputs: B - set of bad patterns (each pattern given as a list)
# Output: "Reduced" set B' of bad patterns where the redundant 
#           patterns from B were removed 
#           (i.e., B' is a minimal subset of B where avoiding B'
#           is equivalent to avoiding B)
#
# Example:
#   RemoveRedundant({[1,2,3],[1,3,2],[1,2,3,4]}) 
#     returns {[1,2,3],[1,3,2]}
RemoveRedundant:=proc(B)
local Bsorted,B1,p1,p2,i,j,k:

Bsorted:=[op(B)]:

Bsorted:=[seq(Bsorted[nops(B)-i],i=0..nops(B)-1)]:

B1:=B:

for i from 1 to nops(B)-1 do
  p1:=Bsorted[i]:

  j:=i+1:
  while j<=nops(B) do
    p2:=Bsorted[j]:

    if nops(p1)>nops(p2) then
      if member(p2,{seq(ScalePerm([op(k..k+nops(p2)-1,p1)]),k=1..nops(p1)-nops(p2)+1)}) then
        B1:=B1 minus {p1}:
        j:=nops(B):
      fi:
    fi:

    j:=j+1:  
  od:
od:

B1:

end:




#-----------------------------------------------------------------------
#
#   SCHEME MODIFICATION PROCEDURES 
#
#-----------------------------------------------------------------------


# inputs a RHS and counts the number of summations
# (sigmas!) in RHS
CountSums:=proc(RHS)
local i:
option remember:

if not type(RHS,list) then
  RETURN(0):
elif type(RHS[1],posint) then
  RETURN(0):
elif RHS[1]=0 then
  RETURN(1+CountSums(RHS[2])):
elif member(RHS[1],{`+`,`*`}) then
  RETURN(add(CountSums(RHS[i]),i=2..nops(RHS))):
else
  RETURN(0):
fi:

end:



# finds the number of functions defined in a scheme
NumFuncInScheme:=proc(Sc)
local n,eq:

n:=1:

for eq in Sc do
  n:=max(n,eq[2][1]):
od:

n:

end:




ExpandScheme:=proc(Sc,k,x,y)
#local eq,n,ESc,b,i,RHS,possibleEqs,dtree,dt,tlist,dtterm,dttermvars,L,fnum:
local eq,n,ESc,b,i,RHS,dtree,dt,tlist,dtterm,dttermvars,L,fnum,taillens,Stmp,xvars,yvars,newRHS:

n:=NumFuncInScheme(Sc):

ESc:=Sc:

# next function name
b:=n+1:

Stmp:={}:
for eq in Sc do
  Stmp:=Stmp union {[eq[2][1], nops(eq[2][3])]}:
od:

taillens:=[seq(Stmp[i][2],i=1..n)]:

for eq in Sc do
  RHS:=eq[3]:

  if CountSums(RHS)>1 then
    ESc:=ESc minus {eq}:

    dtree:=MakeDescentTree(RHS):

    tlist:=[]:

    for dt in dtree do
      dtterm:=ExtractDescent(RHS,dt):
      # dtterm:=ShiftTok2(dtterm,k):
      dttermvars:=ExtractNewVars1(dtterm[3],dtterm[2][3]):
      dttermvars:=[k, [op(dttermvars)]]:
      tlist:=[op(tlist), [dttermvars,dtterm,dt]]:
    od:

    for i from 1 to nops(tlist) do
      fnum:=IsRedundant(ESc,tlist[i][1],tlist[i][2]):
      if fnum>0 then
        RHS:=SubsDescent(RHS,tlist[i][3],[fnum,tlist[i][2][2][2],tlist[i][1][2]]):
      else
        taillens:=[op(taillens), nops(tlist[i][1][2])]:
        RHS:=SubsDescent(RHS,tlist[i][3],[b,tlist[i][2][2][2],tlist[i][1][2]]):
        #ESc:=ESc union {[{true}, [b,op(tlist[i][1])], tlist[i][2]]}:
        xvars:={seq(x[j],j=1..taillens[eq[2][1]])}:
        yvars:=convert(tlist[i][2][2][3],set):
        newRHS:=ShiftTok2(tlist[i][2],k):
        ESc:=ESc union {[{true},[b,tlist[i][1][1],
                         [op(RemoveConstants(tlist[i][1][2],xvars))]],
                         [newRHS[1], newRHS[2], RemoveConstantsEqSet(newRHS[3], {op(xvars), op(yvars)})]]}:
        b:=b+1:
      fi:
    od: 

    ESc:=ESc union {[eq[1], eq[2], RHS]}:

    # Check if a new function name would be redundant
    # Make necessary substitutions
  fi:

od:

ESc:

end:



# Find descent trees
# For each one, see if a substitution would be redundant
# 




# if redundant, return the function "name" where it was previuosly
# defined; if NOT redundant, return -1.
IsRedundant:=proc(Sc,Fin,FRHS)
local n,i:

n:=NumFuncInScheme(Sc):

for i from 1 to n do
  if member([{true}, [i, op(Fin)], FRHS], Sc) then
    RETURN(i):
  fi:
od:

RETURN(-1):

end:





# MAY BE BUGGY!!!


# Ex: S:={x[1]=i[2], x[2]=i[3], i[1]<i[2], i[2]<i[3]}:
#     ExtractNewVars(S,[i[1],i[2],i[3]])
#   will output {x[1], x[2]}
ExtractNewVars1:=proc(S,cvars)
local c,newvars,newvars2,c1:

newvars:={}:

for c in S do
  if type(c,`=`) then
    #newvars:=newvars union RemoveConstants(convert(c,set)):
    newvars:=newvars union convert(c,set):
  fi:
od:

#newvars minus RemoveConstants(convert(cvars,set)):
newvars minus convert(cvars,set):


end:





RemoveConstants:=proc(S,varset)
local S1,s:

S1:={}:

for s in S do
  if member(s,varset) then
    S1:=S1 union {s}:
  elif convert(s,set) intersect varset<>{} then
    S1:=S1 union (convert(s,set) intersect varset):
  fi:
od:

S1:

end:



RemoveConstantsEqSet:=proc(S,varset)
local S1,c,cset,eq,i:

S1:={}:

for c in S do
  if type(c,`=`) then
    cset:=convert(c,list):
    eq:=[]:

    for i from 1 to 2 do
      if member(cset[i],varset) then
        eq:=[op(eq), cset[i]]:
      else
        eq:=[op(eq), op(convert(cset[i],set) intersect varset)]:
      fi:
    od:

    S1:=S1 union {eq[1]=eq[2]}:
  else
    S1:=S1 union {c}:
  fi:

od:

S1:

end:



ShiftTok2:=proc(sum1,k)
[sum1[1], [sum1[2][1], k, sum1[2][3]], sum1[3]]:
end:






# creates a set of descent trees in RHS to the summations in RHS
MakeDescentTree:=proc(RHS)
local nextlevel,i,S,dt:

if not type(RHS,list) then
  RETURN(-1):
elif RHS[1]=0 then
  RETURN({[]}):
elif member(RHS[1],{`+`,`*`}) then
  S:={}:
  for i from 2 to nops(RHS) do
    nextlevel:=MakeDescentTree(RHS[i]):
    if nextlevel<>-1 then
      S:=S union {seq([i,op(dt)], dt in nextlevel)}:
    fi:
  od:
  RETURN(S):
fi:

end:




# extracts the entry in RHS given by the descent tree dt
ExtractDescent:=proc(RHS,dt)

if dt=[] then
  RETURN(RHS):
else
  RETURN(ExtractDescent(RHS[dt[1]], [op(2..nops(dt),dt)])):
fi:

end:



# substitutes the entry in RHS given by the descent tree dt
# with A
SubsDescent:=proc(RHS,dt,A)
local i,L:

if dt=[] then
  RETURN(A):
else
  L:=[seq(RHS[i],i=1..dt[1]-1)]:
  L:=[op(L), SubsDescent(RHS[dt[1]],[op(2..nops(dt),dt)],A)]:
  L:=[op(L), seq(RHS[i],i=dt[1]+1..nops(RHS))]:
  RETURN(L):
fi:

end:




#-----------------------------------------------------------------------
#
#   CLUSTER TAIL GF GENERATION PROCEDURES 
#
#-----------------------------------------------------------------------


ScToTailSc:=proc(Sc,z)
local i,j,l,n,CSc,vars1,vars2,xlist,yset,curRecur,zvars,newfvars,newRecur,RHS:

n:=NumFuncInScheme(Sc):
CSc:=ConvertScheme2(Sc,n):

# Check that all patterns of same size
if n>1 then

  j:=CSc[1][1][2][3]:

  for i from 2 to n do
    if j<>CSc[i][1][2][3] then
      ERROR(`This procedure can only handle sets with patterns of the same length!`):
    fi:
  od:

fi:

newRecur:={}:

# for future generalization with multi-pattern sets
for i from 1 to n do
  curRecur:=CSc[i]:
  xlist:=curRecur[1][2][3]:
  zvars:=[seq(z[j], j=1..nops(xlist))]:
  newfvars:=[curRecur[1][2][1], curRecur[1][2][2], zvars]:

  for j from 1 to nops(curRecur) do

    RHS:=curRecur[j][3]:
    
    if whattype(RHS)<>`list` then
      newRecur:=newRecur union
        {[curRecur[j][1], newfvars, RHS*mul(zvars[l]^l,l=1..nops(xlist))]}:
    else
      newRecur:=newRecur union {[curRecur[j][1], newfvars, ConvertToZ(RHS,xlist,z)]}:
    fi:

  od:

  # WARNING: hard coded- use Descent Trees instead for greater flexibility!

od:

newRecur:


end:



# converts a RHS of a Scheme to the cluster tail form (in z)
ConvertToZ:=proc(RHS,xlist,z)
local i:


if RHS[1]=`+` or RHS[1]=`*` then
  RETURN([RHS[1], seq(ConvertToZ(RHS[i],xlist,z), i=2..nops(RHS))]):
elif RHS[1]=0 then
  # HIT A SUMMATION!
  RETURN(ConvertSumToZ(RHS,xlist,z)):
else
  RETURN(RHS):
fi:


end:




# converts a summation in a scheme to the cluster tail form (via geometric series)
ConvertSumToZ:=proc(RHS,xlist,z)
local i,j,xset,ylist,zvars,cond,vars1,vars2,expr,Slist,simplist,numerlist,curexp,keepinsum,condeqset,newvars,newRHS,yvars,zprod,coeffadj,xi,xiexp,junk:


ylist:=RHS[2][3]:
zvars:=[seq(z[i],i=1..nops(ylist))]:
cond:=RHS[3]:

xset:=convert(xlist,set):

vars1:=ExtractNewVars2(cond,ylist):

condeqset:={}:
yvars:={}:

for i from 1 to nops(cond) do
  if type(cond[i],`=`) then
    condeqset:=condeqset union {convert(cond[i],list)}:
    yvars:=yvars union (convert(cond[i],set) intersect convert(ylist,set)):
  fi:
od:

vars2:={}:



for expr in vars1 do
  if type(expr,`+`) then
    vars2:=vars2 union (convert(expr,set) intersect xset):
  else
    vars2:=vars2 union {expr}:
  fi:
od:


# create Slist
Slist:=MakeSlist(xset,vars2,z):

simplist:=convert(expand(SimplifySums(Slist)),list):


newRHS:=[]:

for i from 1 to nops(simplist) do

  numerlist:=convert(numer(simplist[i]),list):

  keepinsum:=[]:

  for j from 1 to nops(numerlist) do
    if type(numerlist[j],`^`) then
      curexp:=convert(numerlist[j],list):
      if convert(curexp,set) intersect vars2<>{} then
        keepinsum:=[op(keepinsum), curexp]:
      fi:
    fi:

  od:

  newvars:=[]:
  coeffadj:=1:



  for j from 1 to nops(ylist) do
    if member(ylist[j],yvars) then
      zprod:=1:

      # find matching xi for the ylist[j]
      
      for junk in condeqset do
        if junk[2]=ylist[j] then
          xiexp:=junk[1]:
          if type(xiexp,`+`) then
            xi:=op(convert(xiexp,set) intersect xset):
          else
            xi:=xiexp:
          fi:
        elif junk[1]=ylist[j] then
          xiexp:=junk[2]:
          if type(xiexp,`+`) then
            xi:=op(convert(xiexp,set) intersect xset):
          else
            xi:=xiexp:
          fi:
        fi:
      od:



      for curexp in keepinsum do
#        if member(curexp[2],vars2) then
        if curexp[2]=xi then
          zprod:=zprod*curexp[1]:
          coeffadj:=coeffadj*curexp[1]^xiexp:
        fi:
      od:
      newvars:=[op(newvars),zprod]:
    else
      newvars:=[op(newvars),1]:
    fi:

  od:

  newvars:=[RHS[2][1], RHS[2][2], newvars]:

  newRHS:=[op(newRHS), [`*`, expand(simplist[i]/coeffadj), newvars]]:

od:

[`+`, op(newRHS)]:

end:





MakeSlist:=proc(xset,xskip,z)
local usedx,Slist,xi,i,j,upB,upperxskip:

usedx:=xset minus xskip:

#Slist:=[mul(z[op(xi)]^xi, xi in usedx)]:

Slist:=[]:

upperxskip:=xskip:

for j from 1 to nops(usedx) do
  xi:=usedx[j]:
  i:=op(xi):

  upperxskip:=upperxskip minus {op(1..i,xset)}:

  if upperxskip={} then
    upB:=k-(nops(xset)-i):
  else
    upB:=upperxskip[1]:
    upB:=upB-(op(upB)-i):
  fi:

  if i=1 then
    Slist:=[ [xi,1,upB], op(Slist)]:
  else
    Slist:=[ [xi,xset[i-1]+1,upB], op(Slist)]:
  fi:

od:


#[mul(z[op(xi)]^xi, xi in usedx), op(Slist)]:

[mul(z[op(xi)]^xi, xi in xset), op(Slist)]:


end:




##########################
# Extract x[i]'s utilized in SCHEME's RHS (says set R)
# - Use ExtractNewVars1 for "x[i] expressions"
# - Extract actual x[i]'s from expressions

# Create "Slist" for function SimplifySums

# For each summand in "simplified output", factor out remaining z[i]^x[i]'s
# 
# Which yi's have `=` in cond?
# - if not -> 1 entry in recursion
# - if so, find all zj^xj with corresponding xj's
# - divide out by zj^(xj-k) for appropriate k
# 
##########################






# SimplifySums(Slist): inputs a list Slist representing a nested summation (of k sums)
# in the following format:
# 
# Slist = [summand, [sum1_var, sum1_LB, sum1_UB], .., [sumk_var, sumk_LB, sumk_UB] ]
#
# It outputs the summand simplified via Maple's "sum" command.
#
SimplifySums:=proc(Slist) local L1:

if nops(Slist)=1 then
  RETURN(Slist[1]):
else
  L1:=Slist[2]:
  RETURN(SimplifySums([sum(Slist[1], L1[1]=L1[2]..L1[3]), op(3..nops(Slist),Slist)])):
fi:

end:






# Ex: S:={x[1]=i[2], x[2]=i[3], i[1]<i[2], i[2]<i[3]}:
#     ExtractNewVars2(S,[i[1],i[2],i[3]])
#   will output {x[1], x[2]}
ExtractNewVars2:=proc(S,cvars)
local c,newvars,newvars2,c1:

newvars:={}:

for c in S do
  if type(c,`=`) then
    newvars:=newvars union convert(c,set):
  fi:
od:

newvars minus convert(cvars,set):

end:





#-----------------------------------------------------------------------
#  Support/Miscellaneous Procedures
#-----------------------------------------------------------------------

CompactSets:=proc(S)
local Sout,S1,Stmp:

Sout:={}:

for S1 in S do
  Stmp:=TrivEqSets(S1):
  Sout:=Sout union {Stmp[1]}:
od:

Sout:

end:



TrivEqSets:=proc(S)
local Sf,Sc,Sfc,n,p,pf:

Sf:={}:
Sc:={}:
Sfc:={}:

for p in S do
  n:=nops(p):
  Sf:=Sf union {[seq(p[n-i],i=0..n-1)]}:
  Sc:=Sc union {[seq(n-p[i]+1,i=1..n)]}:
  pf:=[seq(p[n-i],i=0..n-1)]:
  Sfc:=Sfc union {[seq(n-pf[i]+1,i=1..n)]}:
od:

{S,Sf,Sc,Sfc}:

end:



# Removes trivial equivalences from set S of single patterns
RemoveTrivEq:=proc(S)
local p,pc,pf,pfc,i,S1,n,Sp:

S1:=S:

for p in S1 do
  n:=nops(p):
  pf:=[seq(p[n-i],i=0..n-1)]:
  pc:=[seq(n-p[i]+1,i=1..n)]:
  pfc:=[seq(n-pf[i]+1,i=1..n)]:
  Sp:={p,pf,pc,pfc}:
  S1:=S1 minus Sp union {Sp[1]}:
od:

S1:

end:



# Checks that B is a set of bad patterns in the proper format
CheckSetB:=proc(B)
local p,n:

if not type(B,set) then
  ERROR(`Input B must be a SET of patterns.`):
fi:

if B={} then
  ERROR(`Input B must be a non-empty set.`):
fi:

for p in B do
  if not type(p,list) then
    ERROR(`All elements of B must be lists. B contained non-list %1.`, p):
  fi:

  n:=nops(p):

  if convert(p,set)<>{seq(i,i=1..n)} then
    ERROR(`Invalid pattern %1.`, p):
  fi:

  if n<3 then
    ERROR(`Patterns must be of length at least 3. Set B contained %1.`,p):
  fi:

od:

end:



#-----------------------------------------------------------------------
#  Print Procedures
#-----------------------------------------------------------------------



# Verbose output for CAV(B,N)
CAVprint:=proc(B,S,L)
printf(`Consecutive Avoidance in Permutations (using CAV) with the Pattern Set %a\n\n\n`, B):
printf(`The following scheme was utilized to count permutations that consecutively avoid %a:\n\n%a\n\n`, B, S):
printf(`The number of permutations of length i (1<=i<=%d) consecutively avoiding %a `, nops(L),B):
printf(`is given by the sequence:\n\n%a\n\n\n`, L):
end:



# Verbose output for CAVt(B,N,t)
CAVtprint:=proc(B,S,L,t)
printf(`Consecutive Avoidance in Permutations (using CAVt) with the Pattern Set %a\n\n\n`, B):
printf(`NOTE: This version keeps track of errors in each permutation as well.\n\n`):
printf(`The following scheme was utilized to count permutations that consecutively avoid %a:\n\n%a\n\n`, B, S):
printf(`The number of permutations of length i (1<=i<=%d) with j instances of patterns in %a `, nops(L),B):
printf(`is given by the coefficient of %a^j in the i-th term in the sequence:\n\n%a\n\n\n`, t, L):
end:



# Verbose output for CAVF(B,N)
CAVFprint:=proc(B,S,L)
printf(`Consecutive Avoidance in Permutations (using CAVF) with the Pattern Set %a\n\n\n`, B):
printf(`The following scheme was utilized to count permutations that consecutively avoid %a:\n\n%a\n\n`, B, S):
printf(`The number of permutations of length i (1<=i<=%d) consecutively avoiding %a `, nops(L),B):
printf(`is given by the sequence:\n\n%a\n\n\n`, L):
end:



# Verbose output for CAVFt(B,N,t)
CAVFtprint:=proc(B,S,L,t)
printf(`Consecutive Avoidance in Permutations (using CAVFt) with the Pattern Set %a\n\n\n`, B):
printf(`NOTE: This version keeps track of errors in each permutation as well.\n\n`):
printf(`The following scheme was utilized to count permutations that consecutively avoid %a:\n\n%a\n\n`, B, S):
printf(`The number of permutations of length i (1<=i<=%d) with j instances of patterns in %a `, nops(L),B):
printf(`is given by the coefficient of %a^j in the i-th term in the sequence:\n\n%a\n\n\n`, t, L):
end:



# Verbose output for CAVT(B,N)
CAVTprint:=proc(B,S,FE,L) global F:

printf(`Consecutive Avoidance in Permutations (using CAVT) with the Pattern Set %a\n\n\n`, B):
printf(`CAVT created the following initial scheme for consecutively avoiding %a:\n\n%a\n\n\n`, B, S):
printf(`Using this, the Cluster Tail GF functional equation was created. `):
printf(`It is encoded as:\n\n%a\n\n\n`,FE):
printf(`A slightly more human friendly form can be written as:`):
printf(`\n\n%a\n\n\n`, TailFE2(FE,[seq(F[i],i=1..NumFuncInScheme(S))])):
printf(`This Cluster Tail GF functional equation was used to count permutations avoiding %a.\n`, B):
printf(`The number of permutations of length i (1<=i<=%d) consecutively avoiding %a `, nops(L),B):
printf(`is given by the sequence:\n\n%a\n\n\n`, L):

end:



# Verbose output for CAVTt(B,N,t)
CAVTtprint:=proc(B,S,FE,L,t) global F:

printf(`Consecutive Avoidance in Permutations (using CAVTt) with the Pattern Set %a\n\n\n`, B):
printf(`NOTE: This version keeps track of errors in each permutation as well.\n\n`):
printf(`CAVTt created the following initial scheme for consecutively avoiding %a:\n\n%a\n\n\n`, B, S):
printf(`Using this, the Cluster Tail GF functional equation was created. `):
printf(`It is encoded as:\n\n%a\n\n\n`,FE):
printf(`A slightly more human friendly form can be written as:`):
printf(`\n\n%a\n\n\n`, TailFE2(FE,[seq(F[i],i=1..NumFuncInScheme(S))])):
printf(`This Cluster Tail GF functional equation was used to count permutations `):
printf(`with various numbers of occurrences/errors from %a.\n`, B):
printf(`The number of permutations of length i (1<=i<=%d) with j instances of patterns in %a `, nops(L),B):
printf(`is given by the coefficient of %a^j in the i-th term in the sequence:\n\n%a\n\n\n`, t, L):

end:



# Verbose output for MapWE(n)
MapWEprint:=proc(n,S) local S1:

printf(`(Non-trivial) Strong Wilf-Equivalences for Length %d Patterns \n\n`, n):

printf(`NOTE: We will disregard the trivial equivalences of reversals and complementations.\n\n`):
printf(`Based off of the OverlapMaps equivalence, we have the following Strong Wilf-Equivalences: \n\n`):

for S1 in S do
  if nops(S1)=1 then
    printf(`Did not find a Strong c-Wilf-Equivalence for %a.\n\n`, op(S1)):
  elif nops(S1)>1 then
    printf(`The following single patterns are Strongly c-Wilf-Equivalent: %a \n`, S1):
    printf(`because they have the same OverlapMap: %a \n\n`, OverlapMaps(S1[1],S1[1])):
  fi:
od:

end:



# Verbose output for RhoApprox(p,N)
RhoApproxprint:=proc(L,N,p)

printf(`Approximate Rho Value for the Single Pattern %a \n\n`, p):

printf(`NOTES: For a single pattern, we assume A(n) ~ %a*(%a^n)*n!\n`, gamma, rho):
printf(`The number of decimal places is currently set to %d.\n\n\n`, Digits):
printf(`The single pattern %a has an approximate rho value: %a\n`, p, L[1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[2]):

end:



# Verbose output for AsymApprox(p,N)
AsymApproxprint:=proc(L,N,p)

printf(`Approximate Asymptotics for the Single Pattern %a \n\n`, p):

printf(`NOTES: For a single pattern, we assume A(n) ~ %a*(%a^n)*n!\n`, gamma, rho):
printf(`The number of decimal places is currently set to %d.\n\n\n`, Digits):
printf(`The single pattern %a has an approximate rho value: %a\n`, p, L[1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n`, N-1, N, L[2]):
printf(`It also has an approximate gamma value: %a\n`, L[3]):
printf(`The error between the n=%d and n=%d approximated gamma terms is %a.\n\n\n`, N-1, N, L[4]):

end:



# Verbose output for RhoApproxRank(n,N)
RhoApproxRankprint:=proc(L,N)
local i,m,k:

m:=nops(L):
k:=nops(L[1][3][1]):

printf(`Ranking of Approximate Rho Values for Length %d Patterns \n\n`, k):

printf(`NOTES: For a single pattern, we assume A(n) ~ %a*(%a^n)*n!\n`, gamma, rho):
printf(`The number of decimal places is currently set to %d.\n\n\n`, Digits):

printf(`The set of single patterns %a produces the largest rho value: %a\n`, L[m][3], L[m][1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[m][2]):

if m>2 then
  for i from 2 to m-1 do
    printf(`The set of single patterns %a produces the next smallest rho value: %a\n`, L[m+1-i][3], L[m+1-i][1]):
    printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[m+1-i][2]):
  od:
fi:

printf(`The set of single patterns %a produces the smallest rho value: %a\n`, L[1][3], L[1][1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[1][2]):


end:



# Verbose output for AsymApproxRank(n,N)
AsymApproxRankprint:=proc(L,N)
local i,m,k:

m:=nops(L):
k:=nops(L[1][5][1]):

printf(`Approximate Asymptotics (Ranked by Rho) for Length %d Patterns \n\n`, k):

printf(`NOTES: For a single pattern, we assume A(n) ~ %a*(%a^n)*n!\n`, gamma, rho):
printf(`The number of decimal places is currently set to %d.\n\n\n`, Digits):

printf(`The set of single patterns %a produces the largest rho value: %a\n`, L[m][5], L[m][1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n`, N-1, N, L[m][2]):
printf(`The coefficient gamma is approximately: %a\n`, L[m][3]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[m][4]):

if m>2 then
  for i from 2 to m-1 do
    printf(`The set of single patterns %a produces the largest rho value: %a\n`, L[m+1-i][5], L[m+1-i][1]):
    printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n`, N-1, N, L[m+1-i][2]):
    printf(`The coefficient gamma is approximately: %a\n`, L[m+1-i][3]):
    printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[m+1-i][4]):
  od:
fi:

printf(`The set of single patterns %a produces the smallest rho value: %a\n`, L[1][5], L[1][1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n`, N-1, N, L[1][2]):
printf(`The coefficient gamma is approximately: %a\n`, L[1][3]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[1][4]):

end:





WilfEqmprint:=proc(L,V,n,N,m)
local Llen,i,Li:

if m>1 then
  printf(`(Non-trivial) Strong c-Wilf-Equivalences for Sets of %d Patterns of Length %d\n\n`, m, n):
else
  printf(`(Non-trivial) Strong c-Wilf-Equivalences for Single Patterns of Length %d\n\n`, n):
fi:

printf(`NOTE: We will disregard the trivial equivalences of reversals and complementations.\n\n\n`):


if V ={} then
  printf(`We were able to attain a TOTAL CLASSIFICATION of Wilf-Equivalences!\n`):
  printf(`The Strong c-Wilf-Equivalences are below `):
  printf(`(and there were no Wilf-Equivalences that were not strong):\n\n\n`):
fi:

Llen:=nops(L):

for i from 1 to Llen do
  Li:=L[Llen+1-i]:
  
  if m=1 then
    if nops(Li[2])>1 then
      printf(`The following single patterns are Strong c-Wilf-Equivalent by OverlapMaps: %q\n`, op(Li[2])):
      printf(`They have the same OverlapMaps: %a\n`, OverlapMaps(Li[2][1],Li[2][1])):
      printf(`For each single pattern, A(%d)=%a\n\n\n`, N, Li[1]):
    else
      printf(`No Strong c-Wilf-Equivalences could be proven for the single pattern: %a\n`, Li[2][1]):
      printf(`It has OverlapMaps: %a\n`, OverlapMaps(Li[2][1],Li[2][1])):
      printf(`For this single pattern, A(%d)=%a\n\n\n`,N,Li[1]):
    fi:
  elif m>1 then
    if nops(Li[2])>1 then
      printf(`The following %d-sets of patterns are Strong c-Wilf-Equivalent by OverlapMaps: %q\n`, m, op(Li[2])):
      printf(`They have the same OverlapMaps: %a\n`, OverlapMaps(Li[2][1],Li[2][1])):
      printf(`For each set of patterns, A(%d)=%a\n\n\n`, N, Li[1]):
    else
      printf(`No Strong c-Wilf-Equivalences could be proven for the set of patterns: %a\n`, Li[2][1]):
      printf(`It has OverlapMaps: %a\n`, OverlapMaps(Li[2][1],Li[2][1])):
      printf(`For this set of patterns, A(%d)=%a\n\n\n`,N,Li[1]):
    fi:


  fi:

od:


if V<>{} then
  if m=1 then
    printf(`The following single patterns MIGHT be Wilf-Equivalent, `):
    printf(`but this could not be proven or disproven.\n`):
    printf(`Trying a larger N value may help in disproving some of the below ambiguities.\n\n`):

    for i from 1 to nops(V) do
      printf(`The following (proven) single pattern Strong c-Wilf-Equivalence classes `):
      printf(`MIGHT also be Wilf-Equivalent: %q\n\n`, op(V[i])):
    od:

  else
    printf(`The following sets of patterns MIGHT be Wilf-Equivalent, `):
    printf(`but this could not be proven or disproven.\n`):
    printf(`Trying a larger N value may help in disproving some of the below ambiguities.\n\n`):

    for i from 1 to nops(V) do
      printf(`MIGHT be Wilf-Equivalent: %q\n\n`, op(V[i])):
    od:
  fi:

  printf(`\n`):

fi:


end:





#-----------------------------------------------------------------------
#  Testing Code
#-----------------------------------------------------------------------


Test132:=proc(n)
local F,C,k,x1,x2,x3,N,j:
option remember:

if n<3 then
  RETURN(n!):
fi:

F:=(k,x1,x2,x3)->piecewise(k=3,-1,k<3,0,-(sum(sum(F(k-2,y1,x1,y3), y1=1..x1-1), y3=1+x1..k-2))):

C:=N->sum(sum(sum(F(N,x1,x2,x3), x3=x2+1..N), x2=x1+1..N-1), x1=1..N-2):

n*Test123(n-1)+add(binomial(n,j)*C(j)*Test123(n-j),j=3..n):

end:




TailGFNaive:=proc(Sc,i,k,N,z)
local j,n,CSc,xn,S,N1,TGF:

n:=NumFuncInScheme(Sc):
CSc:=ConvertScheme2(Sc,n):

xn:=nops(CSc[i][2][2][3]):

S:={op(choose(N,xn))}:
add(ComputeFiNx(CSc,i,N,xvect,k)*mul(z[j]^xvect[j],j=1..xn), xvect in S):

end:




TestRHS132:=proc(Sc,i,k,N,z)
local n,CSc,TGF,x1,x2,x3,y1,y3:

n:=NumFuncInScheme(Sc):
CSc:=ConvertScheme2(Sc,n):

TGF:=0:

for x1 from 1 to N-2 do
  for x2 from x1+1 to N-1 do
    for x3 from x2+1 to N do
      for y1 from 1 to x1-1 do
        for y3 from x1+1 to N-2 do
          TGF:=TGF+ComputeFiNx(CSc,i,N-2,[y1,x1,y3],k)*z[1]^x1*z[2]^x2*z[3]^x3:
        od:
      od:
    od:
  od:
od:

-TGF:

end:





WilfEq:=proc(n,N)
local Sn,Sn1,L,p,Sc,FE,S,Li,i,Lprev,violators,verb:
global k,x,y,z:

verb:=false:

if nops([args])=3 then
  if type(args[3],boolean) then
    verb:=args[3]:
  else
    ERROR(`Optional 3rd argument (verbose) must be of type boolean.`):
  fi:
fi:


Sn:=convert(permute(n),set):
Sn1:=RemoveTrivEq(Sn):

L:={}:

for p in Sn1 do
  Sc:=SCHEME(k,{p},x,y,0):
  FE:=ScToTailSc(Sc,z):
  L:=L union {[CAVOIDtail(FE,N,k), p]}:
od:

L:=convert(L,list):

S:={ [L[1][1], {L[1][2]}]}:

for i from 2 to nops(L) do
  Li:=L[i]:
  Lprev:=S[nops(S)]:
  
  if Li[1]=Lprev[1] then
    S:=S minus {Lprev}:
    S:=S union {[Lprev[1], Lprev[2] union {Li[2]}]}:
  else
    S:=S union {[Li[1], {Li[2]}]}:
  fi:
od:


violators:={}:

for Li in S do
  if nops(Li[2])>1 then
    if SameScheme(k,Li[2],x,y)=false then
      violators:=violators union {Li[2]}:
    fi:
  fi:
od:

if verb then
  print(` Under construction. Please use VERBOSE=false. `):
else
  if violators={} then
    RETURN(S,true):
  else
    RETURN(S,false,violators):
  fi:
fi:

end:





SameScheme:=proc(k,S1,x,y)
local Sc1,Sc2,i:

if nops(S1)=1 then
  RETURN(true):
fi:


Sc1:=SCHEME(k,{S1[1]},x,y,0):

for i from 2 to nops(S1) do
  Sc2:=SCHEME(k,{S1[i]},x,y,0):
  if evalb(Sc1<>Sc2) then
    RETURN(false):
  fi:
od:

RETURN(true):

end:





####End Brian Nakamura's program#########


#ECGFbrian(FP,K): The first K terms of the empirical 
#cluster exponential generating function. For example, try:
#ECGFbrian({[1,3,2]},20);
ECGFbrian:=proc(FP,K) local z,gu,f,i:
gu:=CAVT(FP,K):
f:=1/(1+add(gu[i]*z^i/i!,i=1..K)):
f:=taylor(f,z=0,K+1):
if coeff(f,z,0)<>1 or coeff(f,z,1)<>-1 then
 RETURN(FAIL):
fi:

[0,seq(-coeff(f,z,i)*i!,i=2..K)]:
end:

#AnPlusBrian(N,FP,N1,C,n,E): The first N1 terms of the sequence
#enumerating the number of permutations of length n
#avoiding the set of patterns FP using the first N terms
#and conjecturing a linear recurrence of complexity <=C
#satisfied by
#the coefficients of the cluster generating function
#For example, try:
#AnPlusBrian(20,{[1,2,3]},40,6,n,E);
AnPlusBrian:=proc(N,FP,N1,C,n,E) local ope,gu,f,z,i:
gu:=ECGFbrian(FP,N):
ope:=Findrec(gu,n,E,C):

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

if SeqFromRec(ope,n,E,[op(1..degree(ope,E),gu)],nops(gu))<>gu then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:


gu:=SeqFromRec(ope,n,E,[op(1..degree(ope,E),gu)],N1):


f:=1-z+add(gu[i]*z^i/i!,i=1..N1):

f:=taylor(1/f,z=0,N1+1):
[seq(coeff(f,z,i)*i!,i=1..N1)],ope:


end:




#AnPam(a,m,N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the number of
#occurrences of the  consecutive pattern 1...a tau (a+1) in S_{m+1}
#for any tau.
#For example, try:
#AnPam(1,2,10,t);
AnPam:=proc(a,m,N,t) local gu,x,i,j,k:
gu:=1-x-add((t-1)^k*x^(k*m+1)/(k*m+1)!*mul(binomial(j*m-a,m-a),j=1..k),
k=1..trunc((N-1)/m)+2):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:


#DKfkp(k,p,m,a,b): implementing the recursion (5) of 
#Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance. 
#It is the number of k-clusters for any consecutive pattern
#that starts with a+1 and ends at b+1 such that the cluster
#starts at p+1
#For example, try:
#DKfkp(2,6,1,3,2);
DKfkp:=proc(k,p,m,a,b) local q:
option remember:
if a=b then
 RETURN(0):
fi:
if a>b then
RETURN(DKfkp(k,p,m,b,a)):
fi:

if k=1 then
 if p=a then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:

add(DKfkp(k-1,q-b,m,a,b)*
binomial(p,a)*binomial(k*m-q,m-b)*binomial(q-p-1,b-a-1),q=p+1..m*k):
end:




#DKfk(k,m,a,b): The number of k-clusters of length m*k+1
#for a pattern tau of length m+1 that starts with
#a+1 and ends with b+1 in
#Anton Khoroshkin's paper "Anick-type resolutions and consecutive pattern
#avoidance. 
#For example, try:
#DKfk(2,3,1,2);
DKfk:=proc(k,m,a,b) local p:
add(DKfkp(k,p,m,a,b),p=0..k*m):
end:



#AnPDK(a,b,m,N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the weight enumerator for
#occurrences of the  consecutive pattern tau of length m+1
#with exactly one overlap and
#tau[1]=a+1 and tau[m+1]=b+1 following the formula
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance."
#AnPDK(1,2,3,10,t);
AnPDK:=proc(a,b,m,N,t) local gu,x,i,k:
gu:=1-x-add((t-1)^k*x^(k*m+1)/(k*m+1)!*
DKfk(k,m,a,b),k=1..trunc((N-1)/m)+2):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:


#DK1324cnl(n,l): Formula (6) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK1324cnl(10,3);
DK1324cnl:=proc(n,l) local k:
option remember:
if n=1 then
 if l=0 then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:

if n=2 or n=3 then
 RETURN(0):
fi:

add(binomial(2*k,k)/(k+1)*DK1324cnl(n-2*k-1,l-k),k=1..trunc((n-2)/2)):
end:


#DK1324cn(n,t): Summing the DK1324cnl(n,l) of Formula (6) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK1324cn(10,t);
DK1324cn:=proc(n,t) local l:
expand(add(DK1324cnl(n,l)*(t-1)^l,l=1..n)):
end:


#AnP1324(N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the weight enumerator for
#occurrences of the  consecutive pattern 1324
#taken from Theorem 4 in
#Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For examle, try:
#AnP1324(10,t);
AnP1324:=proc(N,t) local gu,x,i,n:
gu:=1-x-add(x^n/n!*DK1324cn(n,t),n=1..N):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:





#DK1423cnl(n,l): Formula (8) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK1423cnl(10,3);
DK1423cnl:=proc(n,l) local k:
option remember:
if n=1 then
 if l=0 then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:

if n=2 or n=3 then
 RETURN(0):
fi:

add(binomial(n-k-2,k)*DK1423cnl(n-2*k-1,l-k),k=1..trunc((n-2)/2)):
end:


#DK1423cn(n,t): Summing the DK1324cnl(n,l) of Formula (8) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK1423cn(10,t);
DK1423cn:=proc(n,t) local l:
expand(add(DK1423cnl(n,l)*(t-1)^l,l=1..n)):
end:


#AnP1423(N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the weight enumerator for
#occurrences of the  consecutive pattern 1324
#taken from Theorem 5 in
#Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For examle, try:
#AnP1423(10,t);
AnP1423:=proc(N,t) local gu,x,i,n:
gu:=1-x-add(x^n/n!*DK1423cn(n,t),n=1..N):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:


#AnPincr(a,N): The first N terms of the sequence
#enumerating the permutations of length n
#avoiding the consecutive pattern [1,2,3,.., a]
#using the formula of Sergi Elizalde and Marc Noy. For exaple, try:
#AnPincr(3,20);
AnPincr:=proc(a,N) local gu,x,i,k:
gu:=add(x^(k*a)/(k*a)!,k=0..trunc(N/a)+1)-
add(x^(k*a+1)/(k*a+1)!,k=0..trunc(N/a)+1):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:




#DK2413cnlpq(n,l,p,q): Formula (14) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK2413cnlpq(10,3,1,2);
DK2413cnlpq:=proc(n,l,p,q) local r,s:
option remember:
if n=2 or n=3 then
 RETURN(0):
fi:
if n=4 then

 if l=1 and p=2 and q=4 then
    RETURN(1):
  else
    RETURN(0):
 fi:
fi:

add(add(DK2413cnlpq(n-2,l-1,r,s-1),r=1..p-1),s=p+1..q-1)+
(p-1)*add(add(DK2413cnlpq(n-3,l-1,r-1,s-1),r=p+1..s-1),s=p+2..q-1)+
(p-1)*add(add(DK2413cnlpq(n-3,l-1,r-1,s-2),r=p+1..q-1),s=q+1..n):
end:

#DK2413cnl(n,l): used in Theorem 7 of
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK2413cnl(10,3);
DK2413cnl:=proc(n,l) local p,q:
add(add(DK2413cnlpq(n,l,p,q),p=2..q-2),q=3..n-1): 
end:


#DK2413cn(n,t): Summing the DK2413cnl(n,l) of Formula (14) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK2413cn(10,t);
DK2413cn:=proc(n,t) local l:
expand(add(DK2413cnl(n,l)*(t-1)^l,l=1..n)):
end:


#AnP2413(N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the weight enumerator for
#occurrences of the  consecutive pattern 2413
#taken from Theorem 7 in
#Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For examle, try:
#AnP2413(10,t);
AnP2413:=proc(N,t) local gu,x,i,n:
print(`Wrong!`):
gu:=1-x-add(x^n/n!*DK2413cn(n,t),n=1..N):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:






#khaverim(p): all the friends of the pattern p.
#For example, try: khaverim([1,2,3,4]);
khaverim:=proc(p) local i,n:
n:=nops(p):
{p,[seq(p[n+1-i],i=1..nops(p))],
[seq(n+1-p[i],i=1..n)],
[seq(n+1-p[n+1-i],i=1..nops(p))]}:
end:

#Khaverim(S): all the friends of the set of patterns S.
#of the same length
#For example, try: Khaverim({[1,2,3,4]});
Khaverim:=proc(S) local i,n,p:
n:={seq(nops(p), p in S)}:
if nops(n)<>1 then
print(`All patterns in `, S, `should be the same length`):
RETURN(FAIL):
fi:
n:=n[1]:
{S,
{seq([seq(p[n+1-i],i=1..nops(p))], p in S)},
{seq([seq(n+1-p[i],i=1..n)],p in S)},
{seq([seq(n+1-p[n+1-i],i=1..nops(p))],p in S)}
}:
end:







#Theorem(FP): The theorem that enables fast enumeration
#of permutations avoiding the set of patterns FP
#of the same length.
#For example, try:
#Theorem({[1,2,3]});
Theorem:=proc(FP) local x,z,pol,i,lu,pi,S,k,i1,g,mu:
option remember:
if FP={} then
 print(`Set of patterns must be non-empty`):
fi:

S:={seq(nops(pi), pi in FP)}:

if nops(S)>1 then
 print(`the case where all patterns do not have the same length is  not implemented`):
 RETURN(FAIL):
fi:

k:=S[1]:

lu:=Um(x,z,FP):

mu:=ApplyUmbralMatrix(lu,[seq(g[i](seq(x[i1],i1=1..k-1),z),i=1..(k-1)!)],
[seq(x[i1],i1=1..k-1),z]):

pol:=Init(k-1,x,z):

print(`The (ordinary) generating function`):
print(`whose coeff. of z^n tells you the number of permutations of`):
print(`length n avoiding the set of patterns`, FP):
print(` equals `):
print(Sum(g[i](1$(k-1),z),i=1..(k-1)!)):
print(`where the vector of formal power series of length`, (k-1)!):
print(`in the variables`, seq(x[i],i=1..(k-1)!),z):
print(`is the unique vector that`):
print(`satisfies the following system of functional equations`):

for i from 1 to (k-1)! do
print(g[i](seq(x[i1],i1=1..k-1),z)=pol[i]+mu[i]):
od:

end:









#####stuff with general t[q]########
#UmNicet(x,z,k,t,f,g): A nice rendition of 
#an Umbral matrix operator with the input function
#being called f[i] and the output function called g[i]
#for all general patterns of length k
#For example, try:
#UmNicet(x,z,3,t,f,g);
UmNicet:=proc(x,z,k,t,f,g) local gu,i1,j1,lu,mu,lu1:

gu:=UmSt(x,z,k,t):

for i1 from 1 to nops(gu) do
mu:=0:
 for j1 from 1 to nops(gu[i1]) do
 lu:=gu[i1][j1]:
 mu:=mu+add(lu1[1]*f[j1](op(lu1[2])),lu1 in lu):
 od:
 print(g[i1]=mu):
od:

end:


#UmS(x,z,FP): A simplified version of the global
#umbral evolution matrix operator, for enumerating
#permutations avoiding the patterns (of the same length) in FP
#in terms of an umbral matrix and catalytic variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation)
#For example, try:
#UmS(x,z,{[1,2,3]});
UmS:=proc(x,z,FP) local gu,i1,j1:
 gu:=
Um(x,z,FP):
[seq([seq(SimpU(gu[i1][j1]),j1=1..nops(gu[i1]))],i1=1..nops(gu))]:
end:



#UmNice(x,z,FP,f,g): A nice rendition of 
#an Umbral matrix operator with the input function
#being called f[i] and the output function called g[i]
#For example, try:
#UmNice(x,z,{[1,2,3]},f,g);
UmNice:=proc(x,z,FP,f,g) local gu,i1,j1,lu,mu,lu1:

gu:=UmS(x,z,FP):

for i1 from 1 to nops(gu) do
mu:=0:
 for j1 from 1 to nops(gu[i1]) do
 lu:=gu[i1][j1]:
 mu:=mu+add(lu1[1]*f[j1](op(lu1[2])),lu1 in lu):
 od:
 print(g[i1]=mu):
od:

end:





#ECGF(FP,K): The first K terms of the empirical 
#cluster exponential generating function. For example, try:
#ECGF({[1,3,2]},20);
ECGF:=proc(FP,K) local z,gu,f,i:
gu:=AnP(FP,K):
f:=1/(1+add(gu[i]*z^i/i!,i=1..K)):
f:=taylor(f,z=0,K+1):
if coeff(f,z,0)<>1 or coeff(f,z,1)<>-1 then
 RETURN(FAIL):
fi:

[0,seq(-coeff(f,z,i)*i!,i=2..K)]:
end:



####Stuff from Findrec




###Findrec 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrecVerbose:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
if (1+DEGREE)*(1+ORDER)+3+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
print(`There is some slack, there are `, nops(ope)):
print(ope):
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 
#findrec(f,DEGREE,ORDER,n,N): guesses a recurrence operator annihilating
#the sequence f of degree DEGREE and order ORDER
#For example, try: findrec([seq(i,i=1..10)],0,2,n,N);
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,i,j,n0,kv,var1,eq1,mu,a:
option remember:


if not findrecEx(f,DEGREE,ORDER,ithprime(20)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(40)) then
 RETURN(FAIL):
fi:

if not findrecEx(f,DEGREE,ORDER,ithprime(80)) then
 RETURN(FAIL):
fi:


if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:

ope:=subs(var1,ope):

if ope=0 or ope=1 then
  RETURN(FAIL):
fi:

ope:={seq(coeff(expand(ope),kv[i],1),i=1..nops(kv))} minus {0}:

if nops(ope)>1 then
RETURN(Yafe(ope[1],N)[2]):
elif nops(ope)=1 then
RETURN(Yafe(ope[1],N)[2]):
else
 RETURN(FAIL):
fi:

end:
 


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


#Findrec(f,n,N,MaxC): Given a list f tries to find a linear recurrence equation with
#poly coffs.
#of maximum DEGREE+ORDER<=MaxC
#e.g. try Findrec([1,1,2,3,5,8,13,21,34,55,89],n,N,2);
Findrec:=proc(f,n,N,MaxC)
local DEGREE, ORDER,ope,L:

for L from 0 to MaxC do
for ORDER from 0 to L do
 DEGREE:=L-ORDER:
if (2+DEGREE)*(1+ORDER)+4>=nops(f) then
 print(`Insufficient data for degree`, DEGREE, `and order `,ORDER):
 RETURN(FAIL):
fi:
 ope:=findrec([op(1..(2+DEGREE)*(1+ORDER)+4,f)],DEGREE,ORDER,n,N):
     if ope<>FAIL and ope<>1 then
       if  SeqFromRec(ope,n,N,[op(1..ORDER,f)],nops(f))=f then
        RETURN(ope):
       fi:
      fi:
 od:
od:
FAIL:

end:





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

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

gu:=Ini:

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

gu:

end:
 
#end Findrec


with(linalg):

#findrecEx(f,DEGREE,ORDER,m1): Explores whether thre
#is a good chance that there is a recurrence of degree DEGREE
#and order ORDER, using the prime m1
#For example, try: findrecEx([seq(i,i=1..10)],0,2,n,N,1003);
findrecEx:=proc(f,DEGREE,ORDER,m1)
local ope,var,eq,i,j,n0,eq1,a,A1,
D1,E1,Eq,Var,f1,n,N:
option remember:
f1:=f mod m1:
if (1+DEGREE)*(1+ORDER)+5+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f1) mod m1:
  od:
 
   eq:= eq union {eq1}:
od:


Eq:= convert(eq,list):
Var:= convert(var,list):


D1:=nops(Var):
E1:=nops(Eq):
if E1<D1 then
  RETURN(true):
fi:

A1:=matrix(D1,D1):

for i from 1 to D1-1 do
for j from 1 to D1 do
  A1[i,j]:=coeff(Eq[i],Var[j]):
od:
od:

for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1],Var[j]):
od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:

if E1-D1>=1 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+1],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

if E1-D1>=2 then
 for j from 1 to nops(Var) do
  A1[D1,j]:=coeff(Eq[D1+2],Var[j]):
 od:

if det(A1) mod m1 <>0  then
 RETURN(false):
fi:
fi:

true:

end:

####End Stuff from Findrec






#AnPlus(N,FP,N1,C,n,E): The first N1 terms of the sequence
#enumerating the number of permutations of length n
#avoiding the set of patterns FP using the first N terms
#and conjecturing a linear recurrence of complexity <=C
#satisfied by
#the coefficients of the cluster generating function
#For example, try:
#AnPlus(20,{[1,2,3]},40,6);
AnPlus:=proc(N,FP,N1,C,n,E) local ope,gu,f,z:
gu:=ECGF(FP,N):
ope:=Findrec(gu,n,E,C):

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

if SeqFromRec(ope,n,E,[op(1..degree(ope,E),gu)],nops(gu))<>gu then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:


gu:=SeqFromRec(ope,n,E,[op(1..degree(ope,E),gu)],N1):


f:=1-z+add(gu[i]*z^i/i!,i=1..N1):

f:=taylor(1/f,z=0,N1+1):
[seq(coeff(f,z,i)*i!,i=1..N1)],ope:


end:




####Brian Nakamura's program#########

Digits:=100:



#-----------------------------------------------------------------------
#  CAV: A Maple package to count permutations that (consecutively) avoid patterns.
#
#  Created: 
#  This version: 
#
#  Written by Brian Nakamura, Rutgers University
#
#-----------------------------------------------------------------------


#print(` CAV: A Maple package to count permutations that (consecutively) avoid patterns. `):
#print(` Written by: Brian Nakamura, Rutgers University `):
#print(` Version: Oct. 11, 2010 `):
print(``):

#-----------------------------------------------------------------------
#
#   HELP PROCEDURES 
#
#-----------------------------------------------------------------------



Help:=proc()

if args=NULL then
  print(` The main procedures are: `):
  print(` CAV(B,N), CAVF(B,N), CAVT(B,N) `):
  print(` CAVt(B,N,t), CAVFt(B,N,t), CAVTt(B,N,t) `):
  print(` RhoApprox(p,N), AsymApprox(p,N) `):
  print(` MapWilfEq(n), WilfEqm(n,N,m), RhoApproxRank(n,N), AsymApproxRank(n,N) `):
  print(` For help with a specific main procedure, type Help(procedure_name); `):
  print(``):

elif nops([args])=1 and op(1,[args])=CAV then
  printf(`CAV(B,N): inputs a set of patterns B (where each pattern is entered as a list) and `):
  printf(`a positive integer N, and outputs a list L=[a1,a2, .., aN], where ai=# of length i `):
  printf(`permutations that consecutively avoid the patterns in B.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`Also, CAVT does the same thing (through a different approach) MUCH FASTER.\n\n`):
  printf(`For an optional verbose/detailed output, do CAV(B,N,true);\n`):
  printf(`For example, try CAV({[1,2,3],[3,2,1]},13); or CAV({[1,2,3,4]},13,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVt then
  printf(`CAVt(B,N,t): inputs a set of patterns B (where each pattern is entered as a list), `):
  printf(`a positive integer N, and a symbolic variable t.\n\n`):
  printf(`It outputs a list L=[b1,b2,.., bN], where bi is the sum of t^(# of consecutive `):
  printf(`occurrences in p of patterns in B) over all length i permutations p.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`Also, CAVTt does the same thing (through a different approach) MUCH FASTER.\n\n`):
  printf(`For an optional verbose/detailed output, do CAVt(B,N,t,true);\n`):
  printf(`For example, try CAVt({[1,2,3],[3,2,1]},13,t); or CAVt({[1,2,3,4]},13,t,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVF then
  printf(`CAVF(B,N): inputs a set of patterns B (where each pattern is entered as a list) and `):
  printf(`a positive integer N, and outputs a list L=[a1,a2, .., aN], where ai=# of length i `):
  printf(`permutations that consecutively avoid the patterns in B.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`Also, this is occasionally faster than CAV, but CAVT does the same thing `):
  printf(`(through a different approach) MUCH FASTER.\n\n`):
  printf(`For an optional verbose/detailed output, do CAVF(B,N,true);\n`):
  printf(`For example, try CAVF({[1,2,3],[3,2,1]},13); or CAVF({[1,2,3,4]},13,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVFt then
  printf(`CAVFt(B,N,t): inputs a set of patterns B (where each pattern is entered as a list), `):
  printf(`a positive integer N, and a symbolic variable t.\n\n`):
  printf(`It outputs a list L=[b1,b2,.., bN], where bi is the sum of t^(# of consecutive `):
  printf(`occurrences in p of patterns in B) over all length i permutations p.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`Also, this is occasionally faster than CAVt, but CAVTt does the same thing `):
  printf(`(through a different approach) MUCH FASTER.\n\n`):
  printf(`For an optional verbose/detailed output, do CAVFt(B,N,t,true);\n`):
  printf(`For example, try CAVFt({[1,2,3],[3,2,1]},13,t); or CAVFt({[1,2,3,4]},13,t,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVT then
  printf(`CAVT(B,N): inputs a set of patterns B (where each pattern is entered as a list) and `):
  printf(`a positive integer N, and outputs a list L=[a1,a2, .., aN], where ai=# of length i `):
  printf(`permutations that consecutively avoid the patterns in B.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`This procedure uses Cluster Tail functional equations and is generally MUCH FASTER `):
  printf(`than CAV and CAVF!\n\n`):
  printf(`For an optional verbose/detailed output, do CAVT(B,N,true);\n`):
  printf(`For example, try CAVT({[1,2,3],[3,2,1]},13); or CAVT({[1,2,3,4]},13,true);\n\n`):

elif nops([args])=1 and op(1,[args])=CAVTt then
  printf(`CAVTt(B,N,t): inputs a set of patterns B (where each pattern is entered as a list), `):
  printf(`a positive integer N, and a symbolic variable t.\n\n`):
  printf(`It outputs a list L=[b1,b2,.., bN], where bi is the sum of t^(# of consecutive `):
  printf(`occurrences in p of patterns in B) over all length i permutations p.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`This procedure uses Cluster Tail functional equations and is generally MUCH FASTER `):
  printf(`than CAVt and CAVFt!\n\n`):
  printf(`For an optional verbose/detailed output, do CAVTt(B,N,t,true);\n`):
  printf(`For example, try CAVTt({[1,2,3],[3,2,1]},13,t); or CAVTt({[1,2,3,4]},13,t,true);\n\n`):

elif nops([args])=1 and op(1,[args])=RhoApprox then
  printf(`RhoApprox(p,N): inputs a pattern (as a list) and positive integer N>10 and `):
  printf(`computes approximations to rho, the growth constant for the number of permutations `):
  printf(`avoiding pattern p, by counting permutations avoiding p for length N-2, N-1, and N `):
  printf(`permutations.\n\n`):
  printf(`It outputs the approximate rho value followed by the error (difference between consecutive `):
  printf(`rho approximations).\n\n`):
  printf(`NOTE: If a(n)=# of length n permutations avoiding pattern p, then \n a(n) ~ gamma*(rho^n)*n!\n\n`):
  printf(`For an optional verbose/detailed output, do RhoApprox(p,N,true);\n`):
  printf(`For example, try RhoApprox([1,2,3],50); or RhoApprox([1,2,3],50,true);\n\n`):

elif nops([args])=1 and op(1,[args])=AsymApprox then
  printf(`AsymApprox(p,N): inputs a pattern (as a list) and positive integer N>10 and `):
  printf(`computes approximations to rho (the growth constant for the number of permutations `):
  printf(`avoiding pattern p) and gamma (rho's constant coefficient) by counting permutations `):
  printf(`avoiding p for length N-2, N-1, and N permutations.\n\n`):
  printf(`It outputs the approximate rho value, rho's approximate error (difference between `):
  printf(`consecutive rho approximations), the approximate gamma value, and gamma's approximate `):
  printf(`error (difference between consecutive gamma approximations).\n\n`):
  printf(`NOTE: If a(n)=# of length n permutations avoiding pattern p, then \n a(n) ~ gamma*(rho^n)*n!\n\n`):
  printf(`For an optional verbose/detailed output, do AsymApprox(p,N,true);\n`):
  printf(`For example, try AsymApprox([1,2,3],50); or AsymApprox([1,2,3],50,true);\n\n`):

elif nops([args])=1 and op(1,[args])=MapWilfEq then
  printf(`MapWilfEq(n): inputs a postive integer n>=3 and (rigorously!) finds Strong c-Wilf-Equivalences `):
  printf(`via equivalent OverlapMaps (after trivial equivalences of reversals and complementation `):
  printf(`are removed) for single pattern avoidance.\n\n`):
  printf(`It outputs a set of Strong c-Wilf-Equivalence `):
  printf(`classes, where each equivalence class is represented as a set of equivalent patterns.\n\n`):
  printf(`For an optional verbose/detailed output, do MapWilfEq(n,true);\n`):
  printf(`For example, try MapWilfEq(4); or MapWilfEq(4,true);\n\n`):

elif nops([args])=1 and op(1,[args])=WilfEqm then
  printf(`WilfEqm(n,N,m): inputs postive integers n>=3, N>=1, m>=1 `): 
  printf(`and searches for Strong c-Wilf-Equivalences `):
  printf(`via equivalent OverlapMaps (after trivial equivalences of reversals and complementation `):
  printf(`are removed) over sets of m patterns of length n.\n\n`):
  printf(`It outputs a list L=[L1, L2, .. ] and set V, where each Li=[ai, Ai] and `):
  printf(`Ai is the set of length n patterns that is avoided by exactly ai length N permutations.`):
  printf(`Set V is the collection of violators, i.e., the sets of patterns Ai that could not be`):
  printf(`rigorously proven to be Strong c-Wilf-Equivalent by the OverlapMaps approach.\n`):
  printf(`NOTE: if V={}, a rigorous classification was achieved!\n\n`):
  printf(`For an optional verbose/detailed output, do WilfEqm(n,N,m,true);\n`):
  printf(`For example, try WilfEqm(4,15,1); or WilfEqm(4,15,1,true);\n`):
  printf(`Or try WilfEqm(3,10,2); or WilfEqm(3,10,2,true);\n\n`):

elif nops([args])=1 and op(1,[args])=RhoApproxRank then
  printf(`RhoApproxRank(n,N): inputs positive integers n>=3 and N>10 and for each length n pattern p, `):
  printf(`computes approximations to rho, the growth constant for the number of permutations `):
  printf(`avoiding pattern p, by counting permutations avoiding p for length N-2, N-1, and N `):
  printf(`permutations. The rho values are then ranked in increasing order.\n\n`):
  printf(`It outputs a list L=[L1, .., Lk], where each Li is a list with Li=[approx. rho value, `):
  printf(`approx. rho error, set of single patterns that produce these values]. The "approx. rho error" `):
  printf(`is the difference between consecutive rho approximations.\n\n`):
  printf(`NOTE: If a(n)=# of length n permutations avoiding pattern p, then \n a(n) ~ gamma*(rho^n)*n!\n\n`):
  printf(`For an optional verbose/detailed output, do RhoApproxRank(n,N,true);\n`):
  printf(`For example, try RhoApproxRank(4,30); or RhoApproxRank(4,30,true);\n\n`):

elif nops([args])=1 and op(1,[args])=AsymApproxRank then
  printf(`AsymApproxRank(n,N): inputs positive integers n>=3 and N>10 and for each length n pattern p, `):
  printf(`computes approximations to rho (the growth constant for the number of permutations `):
  printf(`avoiding pattern p) and gamma (rho's constant coefficient) by counting permutations `):
  printf(`avoiding p for length N-2, N-1, and N permutations. The values are ranked in increasing `):
  printf(`order for rho.\n\n`):
  printf(`It outputs a list L=[L1, .., Lk], where each Li is a list with Li=[approx. rho value, `):
  printf(`approx. rho error, approx. gamma value, approx. gamma error, set of single patterns that `):
  printf(`produce these values]. The approx. rho (gamma) "error" `):
  printf(`is the difference between consecutive rho (gamma) approximations).\n\n`):
  printf(`NOTE: If a(n)=# of length n permutations avoiding pattern p, then \n a(n) ~ gamma*(rho^n)*n!\n\n`):
  printf(`For an optional verbose/detailed output, do AsymApproxRank(n,N,true);\n`):
  printf(`For example, try AsymApproxRank(4,30); or AsymApproxRank(4,30,true);\n\n`):

else
  ERROR(`Invalid Help parameters!`):
fi:

end:





Help2:=proc()

if args=NULL then
  print(` The primary supporting procedures are: `):
  print(` SCHEME(k,B,x,y,t), MakeTailFE(B,k,z,t), ScToTailSc(Sc,z) `):
  print(` CAVOID(Sc,N,k), CAVOIDtail(TSc,N,k) `):
  print(` OverlapMaps(p1,p2), CompactSets(S), RemoveTrivEq(S) `):
  print(` For help with a specific supporting procedure, type Help2(procedure_name); `):
  print(``):

elif nops([args])=1 and op(1,[args])=SCHEME then
  printf(`SCHEME(k,B,x,y,t): inputs a set of patterns B (where each pattern is entered as a list) `):
  printf(`and symbolic variables k,x,y,t, and outputs the basic scheme used to count permutations `):
  printf(`avoiding patterns in B. See the accompanying paper on more information regarding the scheme.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`If symbolic variable t is inputted, the scheme generated will be to count occurrences of `):
  printf(`patterns in B. If 0 is entered in place of t, the scheme will be for counting permutations that `):
  printf(`completely avoid B.\n`):
  printf(`This scheme can be used directly with CAVOID. To use CAVOIDtail, you must convert this scheme `):
  printf(`into a tail scheme via ScToTailSc first.\n\n`):
  printf(`For example, try SCHEME(k,{[1,2,3,4]},x,y,t); or SCHEME(k,{[1,2,3],[3,2,1]},x,y,0); \n\n`):

elif nops([args])=1 and op(1,[args])=MakeTailFE then
  printf(`MakeTailFE(B,k,z,t): inputs a set of patterns B (where each pattern is entered as a list) `):
  printf(`and symbolic variables k,z,t, and outputs the tail scheme (Cluster Tail GF functional `):
  printf(`equation) for counting permutations avoiding patterns in B. See the accompanying paper `):
  printf(`on more information regarding the tail scheme.\n\n`):
  printf(`NOTES: Patterns in B should be length at least 3.\n`):
  printf(`If symbolic variable t is inputted, the tail scheme generated will be to count occurrences of `):
  printf(`patterns in B. If 0 is entered in place of t, the tail scheme will be for counting permutations `):
  printf(`that completely avoid B.\n`):
  printf(`This tail scheme can be used directly with CAVOIDtail. To use CAVOID, you must generate a basic `):
  printf(`scheme using SCHEME. However, the tail scheme approach is faster.\n\n`):
  printf(`For example, try MakeTailFE({[1,2,3]},k,z,t); or MakeTailFE({[1,2,3],[3,2,1]},k,z,0); \n\n`):

elif nops([args])=1 and op(1,[args])=ScToTailSc then
  printf(`ScToTailSc(Sc,z): inputs a basic scheme Sc (generated from SCHEME) and symbolic variable z, `):
  printf(`and outputs the tail scheme (Cluster Tail GF functional equation) for counting permutations `):
  printf(`avoiding patterns in B. See the accompanying paper on more information regarding the `):
  printf(`tail scheme.\n\n`):
  printf(`NOTE: The input of this procedure must be generated from SCHEME.\n`):
  printf(`The tail scheme generated can be used for CAVOIDtail, which is faster than CAVOID.\n\n`):
  printf(`For example, try ScToTailSc(Sc1,z); where Sc1 is a valid output out of SCHEME.\n\n`):

elif nops([args])=1 and op(1,[args])=CAVOID then
  printf(`CAVOID(Sc,N,k): inputs a basic scheme Sc (generated by SCHEME), a positive integer N, and `):
  printf(`symbolic variable k (that was used in the scheme Sc). It outputs the # of length N permutations `):
  printf(`avoiding a set of patterns (which is essentially encoded in the scheme). IF the scheme was `):
  printf(`created to count occurrences of patterns in B, then the output is b0+b1*t+b2*t^2+... , where `):
  printf(`bi=# of length N permutations with EXACTLY i occurrences of patterns encoded in Sc.\n\n`):
  printf(`NOTES: The input of this procedure must be generated from SCHEME.\n`):
  printf(`The procedure CAVOIDtail, which uses a tail scheme input, is generally much faster.\n\n`):
  printf(`For example, try CAVOID(Sc1,10,k); where Sc1 is a valid output out of SCHEME.\n\n`):

elif nops([args])=1 and op(1,[args])=CAVOIDtail then
  printf(`CAVOIDtail(TSc,N,k): inputs a tail scheme TSc (generated by either MakeTailFE directly or `):
  printf(`converted from a basic scheme by ScToTailSc), a positive integer N, and symbolic variable k `):
  printf(`(that was used in the scheme TSc). It outputs the # of length N permutations avoiding a set `):
  printf(`of patterns (which is essentially encoded in the scheme). IF the scheme was created to count `):
  printf(`occurrences of patterns in B, then the output is b0+b1*t+b2*t^2+... , where `):
  printf(`bi=# of length N permutations with EXACTLY i occurrences of patterns encoded in TSc.\n\n`):
  printf(`NOTES: The input of this procedure must be generated from MakeTailFE or converted from a basic `):
  printf(`scheme using ScToTailSc.\n`):
  printf(`This procedure is generally much faster than CAVOID.\n\n`):
  printf(`For example, try CAVOIDtail(TSc1,10,k); where TSc1 is a valid output of MakeTailFE.\n\n`):

elif nops([args])=1 and op(1,[args])=OverlapMaps then
  printf(`OverlapMaps(p1,p2): inputs two patterns p1 and p2, and outputs a set S={M1,M2, ..} of all `):
  printf(`possible overlaps where a tail of p1 overlaps with a head of p2. Each Mi={[b1,a1],[b2,a2], ..} `):
  printf(`is a mapping for a possible overlap, where p1's tail is a1,a2,.. and p2's head is b1,b2,... \n`):
  printf(`See the accompanying article for more information.\n\n`):
  printf(`For example, try OverlapMaps([1,3,2,4],[1,3,2,4]); or OverlapMaps([1,2,3,4],[2,4,1,3]);\n\n`):

elif nops([args])=1 and op(1,[args])=CompactSets then
  printf(`CompactSets(S): inputs a set S of sets of patterns (where each pattern is entered as a list), `):
  printf(`and outputs S modulo all trivial c-Wilf-Equivalent sets (reversals and complementations).\n\n`):
  printf(`For example, try CompactSets(S); where S is the set of all pairs of length 3 permutations.\n`):
  printf(`(i.e., S=choose(convert(permute(3),set),2); )\n\n`):

elif nops([args])=1 and op(1,[args])=RemoveTrivEq then
  printf(`RemoveTrivEq(S): inputs a set of single patterns S (where each pattern is entered as a list), `):
  printf(`and outputs S modulo all the trivial c-Wilf-Equivalences of reversals and complementations.\n\n`):
  printf(`For example, try RemoveTrivEq(S3); where S3=convert(permute(3),set);\n\n`):

else
  ERROR(`Invalid Help parameters!`):
fi:

end:





#-----------------------------------------------------------------------
#
#   MAIN PROCEDURES 
#
#-----------------------------------------------------------------------


# CAV(B,N): Computes the consecutive avoidance sequence through the general "scheme"
#   approach. NOTE: this is sometimes slower than CAVF and usually slower than CAVT
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the number of length k permutations
#   that avoid (consecutively) the set of patterns B
#
# Ex: CAV({[1,2,3],[3,2,1]}, 10);
#
CAV:=proc(B,N)
local S,i,verb,L:
global k,x,y:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=3 then
  if type(args[3],boolean) then
    verb:=args[3]:
  else
    ERROR(`Optional 3rd argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,0):

L:=[seq(CAVOID(S,i,k),i=1..N)]:


if verb then
  CAVprint(B,S,L):
else
  RETURN(L):
fi:

end:




# CAVt(B,N,t): Analogous to CAV except that the number of occurrences of bad patterns
#   are marked. NOTE: this is sometimes slower than CAVFt and usually slower than CAVTt
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#   t - variable
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the weight enumerated sum over length
#   k permutations, where the weight of permutation p is t^(# occurrences of bad patterns)
#
# Ex: CAVt({[1,2,3],[3,2,1]}, 10, t);
#
CAVt:=proc(B,N,t)
local S,i,verb,L:
global k,x,y:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=4 then
  if type(args[4],boolean) then
    verb:=args[4]:
  else
    ERROR(`Optional 4th argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,t):

L:=expand([seq(CAVOID(S,i,k),i=1..N)]):

if verb then
  CAVtprint(B,S,L,t):
else
  RETURN(L):
fi:

end:




# CAVF(B,N): Computes the consecutive avoidance sequence through the general "scheme"
#   approach, but by defining additional functions in the "scheme" to better utilize "option
#   remember". NOTE: this is sometimes faster than CAV and usually slower than CAVT
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the number of length k permutations
#   that avoid (consecutively) the set of patterns B
#
# Ex: CAVF({[1,2,3],[3,2,1]}, 10);
#
CAVF:=proc(B,N)
local S,i,verb,L:
global k,x,y:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=3 then
  if type(args[3],boolean) then
    verb:=args[3]:
  else
    ERROR(`Optional 3rd argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,0):

L:=[seq(CAVOIDF(S,i,k,x,y),i=1..N)]:

if verb then
  CAVFprint(B,S,L):
else
  RETURN(L):
fi:

end:





# CAVFt(B,N,t): Analogous to CAVF except that the number of occurrences of bad patterns
#   are marked. NOTE: this is sometimes faster than CAVt and usually slower than CAVTt
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#   t - variable
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the weight enumerated sum over length
#   k permutations, where the weight of permutation p is t^(# occurrences of bad patterns)
#
# Ex: CAVFt({[1,2,3],[3,2,1]}, 10, t);
#
CAVFt:=proc(B,N,t)
local S,i,verb,L:
global k,x,y:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=4 then
  if type(args[4],boolean) then
    verb:=args[4]:
  else
    ERROR(`Optional 4th argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,t):

L:=expand([seq(CAVOIDF(S,i,k,x,y),i=1..N)]):

if verb then
  CAVFtprint(B,S,L,t):
else
  RETURN(L):
fi:

end:




# CAVT(B,N): Computes the consecutive avoidance sequence through the general "scheme"
#   approach combined with the cluster tail generating function approach. 
#   NOTE: this is generally faster than CAV and CAVF
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the number of length k permutations
#   that avoid (consecutively) the set of patterns B
#
# Ex: CAVT({[1,2,3],[3,2,1]}, 10);
#
CAVT:=proc(B,N)
local S,i,FE,verb,L:
global k,x,y,z:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=3 then
  if type(args[3],boolean) then
    verb:=args[3]:
  else
    ERROR(`Optional 3rd argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,0):

FE:=ScToTailSc(S,z):

L:=[seq(CAVOIDtail(FE,i,k),i=1..N)]:

if verb then
  CAVTprint(B,S,FE,L):
else
  RETURN(L):
fi:

end:







# CAVTt(B,N,t): Analogous to CAVT except that the number of occurrences of bad patterns
#   are marked. NOTE: this is generally faster than CAV and CAVF
#
# Inputs:
#   B - set of patterns to avoid (where patterns are entered as lists)
#   N - positive integer
#   t - variable
#
# Output:
#   Returns a list [a1,a2, .. , aN] where ak is the weight enumerated sum over length
#   k permutations, where the weight of permutation p is t^(# occurrences of bad patterns)
#
# Ex: CAVTt({[1,2,3],[3,2,1]}, 10, t);
#
CAVTt:=proc(B,N,t)
local S,i,FE,verb,L:
global k,x,y,z:

if not type(N,posint) then
  ERROR(`Input N must be positive integer, but received %1`, N):
fi:

verb:=false:

if nops([args])>=4 then
  if type(args[4],boolean) then
    verb:=args[4]:
  else
    ERROR(`Optional 4th argument (verbose) must be of type boolean.`):
  fi:
fi:

S:=SCHEME(k,B,x,y,t):

FE:=ScToTailSc(S,z):

L:=expand([seq(CAVOIDtail(FE,i,k),i=1..N)]):

if verb then
  CAVTtprint(B,S,FE,L,t):
else
  RETURN(L):
fi:

end:





# RhoApprox(p,N): inputs a SINGLE pattern p, and a pos. integer N
# and outputs the Elizalde constant rho (limit of t[n]/t[n-1]/n
# as n goes to infinity followed by the error.
RhoApprox:=proc(p,N) 
local S,FE,i,gu,mu,verb:
global k,x,y,z:

if N<=10 then
  ERROR(`N must be greater than 10`):
fi:

verb:=false:

if nops([args])>2 then
  if type(op(3,[args]),boolean) then
    verb:=op(3,[args]):
  else
    ERROR(`Optional 3rd input (verbose) must be a Boolean. Instead, received %1.`, op(3,[args])):
  fi:
fi:

S:=SCHEME(k,{p},x,y,0);
FE:=ScToTailSc(S,z):
gu:=[CAVOIDtail(FE,N-2,k),CAVOIDtail(FE,N-1,k),CAVOIDtail(FE,N,k)]:
mu:=[seq(evalf(gu[i+1]/gu[i]/(N-2+i)),i=1..2)]:

if verb then
  RhoApproxprint([mu[2],mu[2]-mu[1]],N,p):
else
  RETURN(mu[2],mu[2]-mu[1]):
fi:

end:






# AsymApprox(p,N): inputs a SINGLE pattern p, and a pos. integer N
# and outputs the Elizalde constant rho (limit of t[n]/t[n-1]/n
# as n goes to infinity, an approximation error for rho, the Elizalde
# coefficient gamma, and an approximation error for gamma
AsymApprox:=proc(p,N)
local S,FE,i,gu,mu,cu,verb:
global k,x,y,z:

if N<=10 then
  ERROR(`N must be greater than 10`):
fi:

verb:=false:

if nops([args])>2 then
  if type(op(3,[args]),boolean) then
    verb:=op(3,[args]):
  else
    ERROR(`Optional 3rd input (verbose) must be a Boolean. Instead, received %1.`, op(3,[args])):
  fi:
fi:

S:=SCHEME(k,{p},x,y,0);
FE:=ScToTailSc(S,z):
gu:=[CAVOIDtail(FE,N-2,k),CAVOIDtail(FE,N-1,k),CAVOIDtail(FE,N,k)]:
mu:=[seq(evalf(gu[i+1]/gu[i]/(N-2+i)),i=1..2)]:
cu:=[gu[2]/(mu[1]^(N-1))/(N-1)!, gu[3]/(mu[2]^N)/N!]:

if verb then
  AsymApproxprint([mu[2],mu[2]-mu[1],cu[2],cu[2]-cu[1]],N,p):
else
  RETURN(mu[2],mu[2]-mu[1],cu[2],cu[2]-cu[1]):
fi:

end:




MapWilfEq:=proc(n)
local Sn,Sn1,p,MapSet,M1,M2,L,eqset,i,verb:

verb:=false:

if nops([args])>1 then
  if type(op(2,[args]),boolean) then
    verb:=op(2,[args]):
  else
    ERROR(`Optional 2nd input (verbose) must be a Boolean. Instead, received %1.`, op(2,[args])):
  fi:
fi:


Sn:=convert(permute(n),set):
Sn1:=RemoveTrivEq(Sn):

MapSet:={}:

for p in Sn1 do
  MapSet:=MapSet union {[OverlapMaps(p,p),p]}:
od:

M1:={}:
L:=[]:
for i from 1 to nops(MapSet) do
  M2:=MapSet[i][1]:

  if M1=M2 then
    eqset:=L[nops(L)] union {MapSet[i][2]}:
    L:=[op(1..nops(L)-1,L), eqset]:
  else
    L:=[op(L), {MapSet[i][2]}]:
  fi:

  M1:=M2:

od:

# Check for verbose output
if verb then
  MapWEprint(n,convert(L,set)):
else
  RETURN(convert(L,set)):
fi:

end:






WilfEqm:=proc(n,N,m)
local Sn,eqSet,eqClass,p,Sc,FE,L,L2,a1,a2,V,i,Li,Lprev,Sm,verb:
global k,x,y,z:


verb:=false:

if nops([args])>3 then
  if type(op(4,[args]),boolean) then
    verb:=op(4,[args]):
  else
    ERROR(`Optional 4th input (verbose) must be a Boolean. Instead, received %1.`, op(4,[args])):
  fi:
fi:

L:={}:

if m=1 then 
  eqSet:=MapWilfEq(n):

  for eqClass in eqSet do
    p:=eqClass[1]:
    Sc:=SCHEME(k,{p},x,y,0):
    FE:=ScToTailSc(Sc,z):
    L:=L union {[CAVOIDtail(FE,N,k), eqClass]}:
  od:

  L:=convert(L,list):

  L2:={ [L[1][1], {L[1][2]}]}:

  for i from 2 to nops(L) do
    Li:=L[i]:
    Lprev:=L2[nops(L2)]:

    if Li[1]=Lprev[1] then
      L2:=L2 minus {Lprev}:
      L2:=L2 union {[Lprev[1], Lprev[2] union {Li[2]}]}:
    else
      L2:=L2 union {[Li[1], {Li[2]}]}:
    fi:
  od:

  V:={}:

  for Li in L2 do
    if nops(Li[2])>1 then
      V:=V union {Li[2]}:
    fi:
  od:

else
  Sm:=choose(convert(permute(n),set),m):
  Sm:=CompactSets(Sm):
  eqSet:={seq({p}, p in Sm)}:

  for eqClass in eqSet do
    p:=eqClass[1]:
    Sc:=SCHEME(k,p,x,y,0):
    FE:=ScToTailSc(Sc,z):
    L:=L union {[CAVOIDtail(FE,N,k), eqClass]}:
  od:

  L:=convert(L,list):

  L2:={ [L[1][1], L[1][2]]}:

  for i from 2 to nops(L) do
    Li:=L[i]:
    Lprev:=L2[nops(L2)]:

    if Li[1]=Lprev[1] then
      L2:=L2 minus {Lprev}:
      L2:=L2 union {[Lprev[1], Lprev[2] union Li[2]]}:
    else
      L2:=L2 union {[Li[1], Li[2]]}:
    fi:
  od:

  V:={}:

  for Li in L2 do
    if nops(Li[2])>1 then
      if SameScheme2(k,Li[2],x,y)=false then
        V:=V union {Li[2]}:
      fi:
    fi:
  od:

fi:


if verb then
  WilfEqmprint(L,V,n,N,m):
else
  RETURN(L,V):
fi:


end:





RhoApproxRank:=proc(n,N)
local eqSet,S,FE,i,gu,mu,p,eqClass,L,verb:
global k,x,y,z:


verb:=false:

if nops([args])>2 then
  if type(op(3,[args]),boolean) then
    verb:=op(3,[args]):
  else
    ERROR(`Optional 3rd input (verbose) must be a Boolean. Instead, received %1.`, op(3,[args])):
  fi:
fi:


eqSet:=MapWilfEq(n):

if N<=10 then
  ERROR(`N must be greater than 10`):
fi:

L:={}:

for eqClass in eqSet do
  p:=eqClass[1]:

  S:=SCHEME(k,{p},x,y,0);
  FE:=ScToTailSc(S,z):
  gu:=[CAVOIDtail(FE,N-2,k),CAVOIDtail(FE,N-1,k),CAVOIDtail(FE,N,k)]:
  mu:=[seq(evalf(gu[i+1]/gu[i]/(N-2+i)),i=1..2)]:
  L:=L union {[mu[2],mu[2]-mu[1],eqClass]}:

od:

if verb then
  RhoApproxRankprint(convert(L,list),N):
else
  RETURN(convert(L,list)):
fi:

end:








AsymApproxRank:=proc(n,N)
local eqSet,S,FE,i,gu,mu,p,eqClass,L,cu,verb:
global k,x,y,z:


verb:=false:

if nops([args])>2 then
  if type(op(3,[args]),boolean) then
    verb:=op(3,[args]):
  else
    ERROR(`Optional 3rd input (verbose) must be a Boolean. Instead, received %1.`, op(3,[args])):
  fi:
fi:

eqSet:=MapWilfEq(n):

if N<=10 then
  ERROR(`N must be greater than 10`):
fi:

L:={}:

for eqClass in eqSet do
  p:=eqClass[1]:

  S:=SCHEME(k,{p},x,y,0);
  FE:=ScToTailSc(S,z):
  gu:=[CAVOIDtail(FE,N-2,k),CAVOIDtail(FE,N-1,k),CAVOIDtail(FE,N,k)]:
  mu:=[seq(evalf(gu[i+1]/gu[i]/(N-2+i)),i=1..2)]:
  cu:=[gu[2]/(mu[1]^(N-1))/(N-1)!, gu[3]/(mu[2]^N)/N!]:
  L:=L union {[mu[2],mu[2]-mu[1],cu[2],cu[2]-cu[1],eqClass]}:

od:

if verb then
  AsymApproxRankprint(convert(L,list),N):
else
  RETURN(convert(L,list)):
fi:


end:






MakeTailFE:=proc(B,k,z,t)
local S,i,FE:
global x,y:

S:=SCHEME(k,B,x,y,t):

FE:=ScToTailSc(S,z):

end:





TailFE2:=proc(TSc,FL)
local S,receq,LHS,RHS:

S:={}:

for receq in TSc do
  LHS:=receq[2]:
  RHS:=receq[3]:
  S:=S union {[receq[1], FL[LHS[1]][LHS[2], LHS[3]] = TailFE2Exp(RHS, FL)]}:
od:

S:

end:




TailFE2Exp:=proc(RHS,FL)
local i:

if type(RHS,list) then
  if type(RHS[1],posint) then
    RETURN(FL[RHS[1]][RHS[2], RHS[3]]):
  elif RHS[1]=`+` then
    RETURN(add(TailFE2Exp(RHS[i],FL), i=2..nops(RHS))):
  elif RHS[1]=`*` then
    RETURN(mul(TailFE2Exp(RHS[i],FL), i=2..nops(RHS))):
  else
    ERROR(`Invalid Tail FE!`):
  fi:
else
  RETURN(RHS):
fi:  

end:






TailFEprint:=proc(TSc,FL)
local TCSc,TFE,n,i:

n:=NumFuncInScheme(TSc):

if n<>nops(FL) then
  ERROR(`The number of function names in list must equal the number of fuctions in Scheme.`):
fi:

TCSc:=ConvertScheme2(TSc,n):

printf(`We have the following Cluster Tail Generating Function(s):\n\n`):

for i from 1 to n do
  TFE:=TailFE2(TSc,FL):

  printf(`%a\n\n`, TFE[2][2]):
  printf(`with initial conditions: %a, if %a\n`, TFE[3][2], op(TFE[3][1])):
  printf(`and %a, if %a\n\n`, TFE[1][2], op(TFE[1][1])):
od:

printf(`These are encoded as the following scheme:\n\n%a`,TSc):
printf(`\n\n`):

end:





# Checks whether every pattern in S1 has the same scheme
SameScheme2:=proc(k,S1,x,y)
local Sc1,Sc2,i:

if nops(S1)=1 then
  RETURN(true):
fi:


Sc1:=SCHEME(k,S1[1],x,y,0):

for i from 2 to nops(S1) do
  Sc2:=SCHEME(k,S1[i],x,y,0):
  if evalb(Sc1<>Sc2) then
    RETURN(false):
  fi:
od:

RETURN(true):

end:





#-----------------------------------------------------------------------
#
#   CLUSTER TAIL SCHEME EXECUTION PROCEDURES 
#
#-----------------------------------------------------------------------


CAVOIDtail:=proc(Sc,N,k)
local j,CSc,n:
option remember:

if N<3 then
  RETURN(N!):
fi:

n:=NumFuncInScheme(Sc):

CSc:=ConvertScheme2(Sc,n):

N*CAVOIDtail(Sc,N-1,k)+
add(binomial(N,j)*ComputeTailCN(CSc,j,n,k)*CAVOIDtail(Sc,N-j,k),j=1..N):


end:






ComputeTailCN:=proc(CSc,N,n,k) local i,j,zvars,zcond:
option remember:

zvars:={}:
zcond:={}:

for i from 1 to n do
  zvars:=zvars union convert(CSc[i][1][2][3],set):
od:

zcond:={seq(zvars[i]=1, i=1..nops(zvars))}:

subs(zcond, simplify(expand(add(ComputeTailFiN(CSc,i,N,k),i=1..n)))):

end:


ComputeTailFiN:=proc(CSc,i,N,k)
local Ri,eq,RHS,total,j:
option remember:

Ri:=CSc[i]:

for eq in Ri do
  if evalb(subs(k=N,op(eq[1]))) then
    RHS:=eq[3]:

    if not type(RHS,list) then
      RETURN(subs(k=N,RHS)):
    else
      if member(RHS[1],{0,`+`,`*`}) then
        RETURN(EvalTailExpression(CSc,RHS,N,k)):
      else
        RETURN(ComputeTailFiN(CSc,RHS[1],subs(k=N,RHS[2]),k)):
      fi:
    fi:
  fi:

od:


end:



EvalTailExpression:=proc(CSc,RHS,N,k)
local j,i,zvars,zcond:
option remember:

if not type(RHS,list) then
  RETURN(subs(k=N,RHS)):
elif evalb(RHS[1]=`+`) then
  RETURN(add(EvalTailExpression(CSc,RHS[j],N,k),j=2..nops(RHS))):
elif evalb(RHS[1]=`*`) then
  RETURN(mul(EvalTailExpression(CSc,RHS[j],N,k),j=2..nops(RHS))):
else

  zcond:={}:

  i:=RHS[1]:

  zvars:=CSc[i][1][2][3]:

  zcond:={seq(zvars[i]=RHS[3][i], i=1..nops(zvars))}:

  RETURN(subs(zcond, simplify(expand(ComputeTailFiN(CSc,RHS[1],subs(k=N,RHS[2]),k))))):

fi:


end:




#-----------------------------------------------------------------------
#
#   (GENERAL) SCHEME EXECUTION PROCEDURES 
#
#-----------------------------------------------------------------------


CAVOID:=proc(Sc,N,k)
local j,CSc,n:
option remember:

if N<3 then
  RETURN(N!):
fi:

n:=NumFuncInScheme(Sc):

CSc:=ConvertScheme2(Sc,n):

N*CAVOID(Sc,N-1,k)+
add(binomial(N,j)*ComputeCN(CSc,j,n,k)*CAVOID(Sc,N-j,k),j=1..N):


end:





CAVOIDF:=proc(Sc,N,k,x,y)
local j,CESc,n1,n2,ESc:
option remember:

if N<3 then
  RETURN(N!):
fi:

n1:=NumFuncInScheme(Sc):

ESc:=ExpandScheme(Sc,k,x,y):

n2:=NumFuncInScheme(ESc):

CESc:=ConvertScheme2(ESc,n2):

CAVOIDF2(CESc,N,k,n1):

end:




CAVOIDF2:=proc(CESc,N,k,n1)
local j:
option remember:

if N<3 then
  RETURN(N!):
fi:

N*CAVOIDF2(CESc,N-1,k,n1)+
add(binomial(N,j)*ComputeCN(CESc,j,n1,k)*CAVOIDF2(CESc,N-j,k,n1),j=1..N):

end:



# computes the sum of weights of all length k clusters
ComputeCN:=proc(CSc,N,n,k) local i:
option remember:
add(ComputeFiN(CSc,i,N,k),i=1..n):
end:


# computes the sum of weights of all length k clusters that end in the i-th pattern
ComputeFiN:=proc(CSc,i,N,k) local xn,S,xvect:
option remember:
xn:=nops(CSc[i][2][2][3]):
S:={op(choose(N,xn))}:
add(ComputeFiNx(CSc,i,N,xvect,k), xvect in S):
end:




# computes the sum of weights of all length k clusters that end in the i-th pattern
# given by xvect
ComputeFiNx:=proc(CSc,i,N,xvect,k)
local Ri,eq,RHS,total,j:
option remember:

Ri:=CSc[i]:

for eq in Ri do
  if evalb(subs(k=N,op(eq[1]))) then
    RHS:=eq[3]:

    if not type(RHS,list) then
      RETURN(RHS):
    else
      if member(RHS[1],{0,`+`,`*`}) then
        RETURN(EvalExpression(CSc,RHS,N,xvect,k,eq[2][3])):
      else
        RETURN(ComputeFiNx(CSc,RHS[1],subs(k=N,RHS[2]),RHS[3],k)):
      fi:
    fi:
  fi:

od:

end:



# evalutaes a summation in a scheme
ComputeSummation:=proc(CSc,sum1,N,xvect,k,xlist)
local S,N2,Sfilt,sj:
option remember:

N2:=subs(k=N,sum1[2][2]):
S:={op(choose(N2,nops(sum1[2][3])))}:

Sfilt:={}:

for sj in S do
  if CheckCond2(xvect,sj,sum1[3],xlist,sum1[2][3]) then
    Sfilt:=Sfilt union {sj}:
  fi:
od:

add(ComputeFiNx(CSc,sum1[2][1],N2,sj,k), sj in Sfilt):

end:



# evaluates a general expression in a scheme
EvalExpression:=proc(CSc,RHS,N,xvect,k,xlist)
local j,cond:
option remember:

if not type(RHS,list) then
  RETURN(RHS):
elif RHS[1]=0 then
  RETURN(ComputeSummation(CSc,RHS,N,xvect,k,xlist)):
elif evalb(RHS[1]=`+`) then
  RETURN(add(EvalExpression(CSc,RHS[j],N,xvect,k,xlist),j=2..nops(RHS))):
elif evalb(RHS[1]=`*`) then
  RETURN(mul(EvalExpression(CSc,RHS[j],N,xvect,k,xlist),j=2..nops(RHS))):
else
  cond:={seq(xlist[j]=xvect[j],j=1..nops(xvect))}:
  RETURN(ComputeFiNx(CSc,RHS[1],subs(k=N,RHS[2]),subs(cond,RHS[3]),k)):

fi:

end:



# Converts a Scheme to a Converted Scheme so that it can quickly find the
# necessary recurrence for the i-th function (1<=i<=n)
ConvertScheme2:=proc(Sc,n)
local i,CSc:
option remember:

for i from 1 to n do
  CSc[i]:={}:
od:

for i from 1 to nops(Sc) do
  CSc[Sc[i][2][1]]:=CSc[Sc[i][2][1]] union {Sc[i]}:
od:

CSc:

end:



CheckCond2:=proc(xvect,ivect,cond,xlist,ilist)
local c1,S,j:
option remember:

S:={seq(xlist[j]=xvect[j],j=1..nops(xvect)), seq(ilist[j]=ivect[j],j=1..nops(ivect))}:

for c1 in cond do
  if not evalb(subs(S,c1)) then
    RETURN(false):
  fi:
od:

RETURN(true):

end:




#-----------------------------------------------------------------------
#
#   SCHEME GENERATION PROCEDURES 
#
#-----------------------------------------------------------------------


SCHEME:=proc(k,B,x,y,t)
local B1,Sc,i,j,r,xvect,EqFi,RecurFi:

CheckSetB(B):

# Bad pattern set B modulo redundancies
B1:=[op(RemoveRedundant(B))]:

Sc:={}:

for i from 1 to nops(B1) do
  r:=nops(B1[i]):
  xvect:=[seq(x[j],j=1..r)]:

# MODDED
  EqFi:={[{k<r}, [i,k,xvect], 0], [{k=r}, [i,k,xvect], t-1]}:

  RecurFi:=MakeRecurrenceFi(k,B1,x,y,i,t):

  # Need if statement to either add RecurFi to EqFi or not
  # Maybe this is not necessary?
  # If this fails, should all of EqFi be scrapped?


# MODDED
  if RecurFi<>[`*`, t-1, []] then
    EqFi:=EqFi union {[{k>r}, [i,k,xvect], RecurFi]}:
    Sc:=Sc union EqFi:
  fi:

od:

Sc:

end:




MakeRecurrenceFi:=proc(k,B1,x,y,i,t)
local RHS,RHSterms,j:

RHS:=[]:

for j from 1 to nops(B1) do
  RHSterms:=FindSums2(k,B1,x,y,j,i):
  RHS:=[op(RHS), op(RHSterms)]:
od:

if nops(RHS)>1 then
  RHS:=[`+`, op(RHS)]:
elif nops(RHS)=1 then
  RHS:=op(RHS):
fi:

# MODDED
[`*`, t-1, RHS]:

end:




FindSums2:=proc(k,B1,x,y,j,i)
local M,ir,jr,S,M1,m1,k2,b,chopped,t,cond:

M:=OverlapMaps(B1[j],B1[i]):

if M={} then
  RETURN({}):
fi:

# lengths of bad words
ir:=nops(B1[i]):
jr:=nops(B1[j]):

S:={}:

for M1 in M do
  k2:=k-ir+nops(M1):

  chopped:={seq(b,b=1..ir)} minus {seq(m1[1],m1 in M1)}:

  cond:={ seq(y[t]<y[t+1], t=1..jr-1) }:

  for m1 in M1 do
    # m1[1] -> m1[2]
    cond:=cond union {x[m1[1]]-nops(chopped intersect {seq(b,b=1..m1[1])})=y[m1[2]]}:
  od:

  S:=S union {[0, [j, k2, [seq(y[t], t=1..jr)]], cond]}:
od:

S:

end:





# OverlapMaps(p1,p2)
#
# Inputs: p1,p2 - permutation patterns
# Output: all "maps" where the tail/end of p1 can overlap with
#           the head/beginning of p2 in a larger permutation
#
# Comments:
#   The output is a set of sets {M1, M2, .., Mk} with each Mi
#   a set of lists. For example, if M1={[a1,b1], [a2,b2], ..},
#   then the a1-th term of p1 coincides with the b1-th term of p2,
#   the a2-th term of p1 coincides with the b2-th term of p2, etc.
#
# Example:
#   OverlapMaps([1,2,3], [1,2,3]) 
#     outputs {{[1, 3]}, {[1, 2], [2, 3]}}
# since it can have an overlap of length 1 (shown by {[1,3]})
#        1 2 3
#    1 2 3
# or an overlap of length 2 (shown by {[1, 2], [2, 3]})
#      1 2 3
#    1 2 3
#
OverlapMaps:=proc(p1,p2)
local i,Tails,Heads,OL,olap,S,S1,k1,k2,m:

k1:=nops(p1):
k2:=nops(p2):

Tails:={}:
Heads:={}:

for i from 1 to k1-1 do
  Tails:={op(Tails), ScalePerm([op(i+1..k1,p1)])}:
od:

for i from 1 to k2-1 do
  Heads:={op(Heads), ScalePerm([op(1..i,p2)])}:
od:

OL:=Heads intersect Tails:

S:={}:

for olap in OL do

  # length of the overlap (in terms of entries in p2)
  m:=nops(olap):

  S1:={}:

  for i from 1 to m do
    S1:={op(S1), [p2[i],p1[k1-m+i]]}:
  od:

  S:={op(S), S1}:

od:

S:

end:







# ScalePerm(p)
#
# Inputs: p - a subsequence of a permutation 
# Output: the canonical reduction of p
#
# Comments: 
#   If p contains k terms, the output will be the permutation
#   of {1,2,..,k} that follows the same respective ordering 
#   as in p
#
# Example:
#   ScalePerm([1,8,4,6]) returns [1,4,2,3]
ScalePerm:=proc(p)
local n,i,k,psorted,L:

n:=nops(p):

psorted:={op(p)}:
psorted:=[op(psorted)]:

L:=[]:

for i from 1 to n do
  k:=1:
  
  while k<=n do
    if psorted[k]=p[i] then
      L:=[op(L), k]:
      k:=n+1:
    else
      k:=k+1:
    fi:
  od:
  
od:

L:

end:



# RemoveRedundant(B)
#
# Inputs: B - set of bad patterns (each pattern given as a list)
# Output: "Reduced" set B' of bad patterns where the redundant 
#           patterns from B were removed 
#           (i.e., B' is a minimal subset of B where avoiding B'
#           is equivalent to avoiding B)
#
# Example:
#   RemoveRedundant({[1,2,3],[1,3,2],[1,2,3,4]}) 
#     returns {[1,2,3],[1,3,2]}
RemoveRedundant:=proc(B)
local Bsorted,B1,p1,p2,i,j,k:

Bsorted:=[op(B)]:

Bsorted:=[seq(Bsorted[nops(B)-i],i=0..nops(B)-1)]:

B1:=B:

for i from 1 to nops(B)-1 do
  p1:=Bsorted[i]:

  j:=i+1:
  while j<=nops(B) do
    p2:=Bsorted[j]:

    if nops(p1)>nops(p2) then
      if member(p2,{seq(ScalePerm([op(k..k+nops(p2)-1,p1)]),k=1..nops(p1)-nops(p2)+1)}) then
        B1:=B1 minus {p1}:
        j:=nops(B):
      fi:
    fi:

    j:=j+1:  
  od:
od:

B1:

end:




#-----------------------------------------------------------------------
#
#   SCHEME MODIFICATION PROCEDURES 
#
#-----------------------------------------------------------------------


# inputs a RHS and counts the number of summations
# (sigmas!) in RHS
CountSums:=proc(RHS)
local i:
option remember:

if not type(RHS,list) then
  RETURN(0):
elif type(RHS[1],posint) then
  RETURN(0):
elif RHS[1]=0 then
  RETURN(1+CountSums(RHS[2])):
elif member(RHS[1],{`+`,`*`}) then
  RETURN(add(CountSums(RHS[i]),i=2..nops(RHS))):
else
  RETURN(0):
fi:

end:



# finds the number of functions defined in a scheme
NumFuncInScheme:=proc(Sc)
local n,eq:

n:=1:

for eq in Sc do
  n:=max(n,eq[2][1]):
od:

n:

end:




ExpandScheme:=proc(Sc,k,x,y)
#local eq,n,ESc,b,i,RHS,possibleEqs,dtree,dt,tlist,dtterm,dttermvars,L,fnum:
local eq,n,ESc,b,i,RHS,dtree,dt,tlist,dtterm,dttermvars,L,fnum,taillens,Stmp,xvars,yvars,newRHS:

n:=NumFuncInScheme(Sc):

ESc:=Sc:

# next function name
b:=n+1:

Stmp:={}:
for eq in Sc do
  Stmp:=Stmp union {[eq[2][1], nops(eq[2][3])]}:
od:

taillens:=[seq(Stmp[i][2],i=1..n)]:

for eq in Sc do
  RHS:=eq[3]:

  if CountSums(RHS)>1 then
    ESc:=ESc minus {eq}:

    dtree:=MakeDescentTree(RHS):

    tlist:=[]:

    for dt in dtree do
      dtterm:=ExtractDescent(RHS,dt):
      # dtterm:=ShiftTok2(dtterm,k):
      dttermvars:=ExtractNewVars1(dtterm[3],dtterm[2][3]):
      dttermvars:=[k, [op(dttermvars)]]:
      tlist:=[op(tlist), [dttermvars,dtterm,dt]]:
    od:

    for i from 1 to nops(tlist) do
      fnum:=IsRedundant(ESc,tlist[i][1],tlist[i][2]):
      if fnum>0 then
        RHS:=SubsDescent(RHS,tlist[i][3],[fnum,tlist[i][2][2][2],tlist[i][1][2]]):
      else
        taillens:=[op(taillens), nops(tlist[i][1][2])]:
        RHS:=SubsDescent(RHS,tlist[i][3],[b,tlist[i][2][2][2],tlist[i][1][2]]):
        #ESc:=ESc union {[{true}, [b,op(tlist[i][1])], tlist[i][2]]}:
        xvars:={seq(x[j],j=1..taillens[eq[2][1]])}:
        yvars:=convert(tlist[i][2][2][3],set):
        newRHS:=ShiftTok2(tlist[i][2],k):
        ESc:=ESc union {[{true},[b,tlist[i][1][1],
                         [op(RemoveConstants(tlist[i][1][2],xvars))]],
                         [newRHS[1], newRHS[2], RemoveConstantsEqSet(newRHS[3], {op(xvars), op(yvars)})]]}:
        b:=b+1:
      fi:
    od: 

    ESc:=ESc union {[eq[1], eq[2], RHS]}:

    # Check if a new function name would be redundant
    # Make necessary substitutions
  fi:

od:

ESc:

end:



# Find descent trees
# For each one, see if a substitution would be redundant
# 




# if redundant, return the function "name" where it was previuosly
# defined; if NOT redundant, return -1.
IsRedundant:=proc(Sc,Fin,FRHS)
local n,i:

n:=NumFuncInScheme(Sc):

for i from 1 to n do
  if member([{true}, [i, op(Fin)], FRHS], Sc) then
    RETURN(i):
  fi:
od:

RETURN(-1):

end:





# MAY BE BUGGY!!!


# Ex: S:={x[1]=i[2], x[2]=i[3], i[1]<i[2], i[2]<i[3]}:
#     ExtractNewVars(S,[i[1],i[2],i[3]])
#   will output {x[1], x[2]}
ExtractNewVars1:=proc(S,cvars)
local c,newvars,newvars2,c1:

newvars:={}:

for c in S do
  if type(c,`=`) then
    #newvars:=newvars union RemoveConstants(convert(c,set)):
    newvars:=newvars union convert(c,set):
  fi:
od:

#newvars minus RemoveConstants(convert(cvars,set)):
newvars minus convert(cvars,set):


end:





RemoveConstants:=proc(S,varset)
local S1,s:

S1:={}:

for s in S do
  if member(s,varset) then
    S1:=S1 union {s}:
  elif convert(s,set) intersect varset<>{} then
    S1:=S1 union (convert(s,set) intersect varset):
  fi:
od:

S1:

end:



RemoveConstantsEqSet:=proc(S,varset)
local S1,c,cset,eq,i:

S1:={}:

for c in S do
  if type(c,`=`) then
    cset:=convert(c,list):
    eq:=[]:

    for i from 1 to 2 do
      if member(cset[i],varset) then
        eq:=[op(eq), cset[i]]:
      else
        eq:=[op(eq), op(convert(cset[i],set) intersect varset)]:
      fi:
    od:

    S1:=S1 union {eq[1]=eq[2]}:
  else
    S1:=S1 union {c}:
  fi:

od:

S1:

end:



ShiftTok2:=proc(sum1,k)
[sum1[1], [sum1[2][1], k, sum1[2][3]], sum1[3]]:
end:






# creates a set of descent trees in RHS to the summations in RHS
MakeDescentTree:=proc(RHS)
local nextlevel,i,S,dt:

if not type(RHS,list) then
  RETURN(-1):
elif RHS[1]=0 then
  RETURN({[]}):
elif member(RHS[1],{`+`,`*`}) then
  S:={}:
  for i from 2 to nops(RHS) do
    nextlevel:=MakeDescentTree(RHS[i]):
    if nextlevel<>-1 then
      S:=S union {seq([i,op(dt)], dt in nextlevel)}:
    fi:
  od:
  RETURN(S):
fi:

end:




# extracts the entry in RHS given by the descent tree dt
ExtractDescent:=proc(RHS,dt)

if dt=[] then
  RETURN(RHS):
else
  RETURN(ExtractDescent(RHS[dt[1]], [op(2..nops(dt),dt)])):
fi:

end:



# substitutes the entry in RHS given by the descent tree dt
# with A
SubsDescent:=proc(RHS,dt,A)
local i,L:

if dt=[] then
  RETURN(A):
else
  L:=[seq(RHS[i],i=1..dt[1]-1)]:
  L:=[op(L), SubsDescent(RHS[dt[1]],[op(2..nops(dt),dt)],A)]:
  L:=[op(L), seq(RHS[i],i=dt[1]+1..nops(RHS))]:
  RETURN(L):
fi:

end:




#-----------------------------------------------------------------------
#
#   CLUSTER TAIL GF GENERATION PROCEDURES 
#
#-----------------------------------------------------------------------


ScToTailSc:=proc(Sc,z)
local i,j,l,n,CSc,vars1,vars2,xlist,yset,curRecur,zvars,newfvars,newRecur,RHS:

n:=NumFuncInScheme(Sc):
CSc:=ConvertScheme2(Sc,n):

# Check that all patterns of same size
if n>1 then

  j:=CSc[1][1][2][3]:

  for i from 2 to n do
    if j<>CSc[i][1][2][3] then
      ERROR(`This procedure can only handle sets with patterns of the same length!`):
    fi:
  od:

fi:

newRecur:={}:

# for future generalization with multi-pattern sets
for i from 1 to n do
  curRecur:=CSc[i]:
  xlist:=curRecur[1][2][3]:
  zvars:=[seq(z[j], j=1..nops(xlist))]:
  newfvars:=[curRecur[1][2][1], curRecur[1][2][2], zvars]:

  for j from 1 to nops(curRecur) do

    RHS:=curRecur[j][3]:
    
    if whattype(RHS)<>`list` then
      newRecur:=newRecur union
        {[curRecur[j][1], newfvars, RHS*mul(zvars[l]^l,l=1..nops(xlist))]}:
    else
      newRecur:=newRecur union {[curRecur[j][1], newfvars, ConvertToZ(RHS,xlist,z)]}:
    fi:

  od:

  # WARNING: hard coded- use Descent Trees instead for greater flexibility!

od:

newRecur:


end:



# converts a RHS of a Scheme to the cluster tail form (in z)
ConvertToZ:=proc(RHS,xlist,z)
local i:


if RHS[1]=`+` or RHS[1]=`*` then
  RETURN([RHS[1], seq(ConvertToZ(RHS[i],xlist,z), i=2..nops(RHS))]):
elif RHS[1]=0 then
  # HIT A SUMMATION!
  RETURN(ConvertSumToZ(RHS,xlist,z)):
else
  RETURN(RHS):
fi:


end:




# converts a summation in a scheme to the cluster tail form (via geometric series)
ConvertSumToZ:=proc(RHS,xlist,z)
local i,j,xset,ylist,zvars,cond,vars1,vars2,expr,Slist,simplist,numerlist,curexp,keepinsum,condeqset,newvars,newRHS,yvars,zprod,coeffadj,xi,xiexp,junk:


ylist:=RHS[2][3]:
zvars:=[seq(z[i],i=1..nops(ylist))]:
cond:=RHS[3]:

xset:=convert(xlist,set):

vars1:=ExtractNewVars2(cond,ylist):

condeqset:={}:
yvars:={}:

for i from 1 to nops(cond) do
  if type(cond[i],`=`) then
    condeqset:=condeqset union {convert(cond[i],list)}:
    yvars:=yvars union (convert(cond[i],set) intersect convert(ylist,set)):
  fi:
od:

vars2:={}:



for expr in vars1 do
  if type(expr,`+`) then
    vars2:=vars2 union (convert(expr,set) intersect xset):
  else
    vars2:=vars2 union {expr}:
  fi:
od:


# create Slist
Slist:=MakeSlist(xset,vars2,z):

simplist:=convert(expand(SimplifySums(Slist)),list):


newRHS:=[]:

for i from 1 to nops(simplist) do

  numerlist:=convert(numer(simplist[i]),list):

  keepinsum:=[]:

  for j from 1 to nops(numerlist) do
    if type(numerlist[j],`^`) then
      curexp:=convert(numerlist[j],list):
      if convert(curexp,set) intersect vars2<>{} then
        keepinsum:=[op(keepinsum), curexp]:
      fi:
    fi:

  od:

  newvars:=[]:
  coeffadj:=1:



  for j from 1 to nops(ylist) do
    if member(ylist[j],yvars) then
      zprod:=1:

      # find matching xi for the ylist[j]
      
      for junk in condeqset do
        if junk[2]=ylist[j] then
          xiexp:=junk[1]:
          if type(xiexp,`+`) then
            xi:=op(convert(xiexp,set) intersect xset):
          else
            xi:=xiexp:
          fi:
        elif junk[1]=ylist[j] then
          xiexp:=junk[2]:
          if type(xiexp,`+`) then
            xi:=op(convert(xiexp,set) intersect xset):
          else
            xi:=xiexp:
          fi:
        fi:
      od:



      for curexp in keepinsum do
#        if member(curexp[2],vars2) then
        if curexp[2]=xi then
          zprod:=zprod*curexp[1]:
          coeffadj:=coeffadj*curexp[1]^xiexp:
        fi:
      od:
      newvars:=[op(newvars),zprod]:
    else
      newvars:=[op(newvars),1]:
    fi:

  od:

  newvars:=[RHS[2][1], RHS[2][2], newvars]:

  newRHS:=[op(newRHS), [`*`, expand(simplist[i]/coeffadj), newvars]]:

od:

[`+`, op(newRHS)]:

end:





MakeSlist:=proc(xset,xskip,z)
local usedx,Slist,xi,i,j,upB,upperxskip:

usedx:=xset minus xskip:

#Slist:=[mul(z[op(xi)]^xi, xi in usedx)]:

Slist:=[]:

upperxskip:=xskip:

for j from 1 to nops(usedx) do
  xi:=usedx[j]:
  i:=op(xi):

  upperxskip:=upperxskip minus {op(1..i,xset)}:

  if upperxskip={} then
    upB:=k-(nops(xset)-i):
  else
    upB:=upperxskip[1]:
    upB:=upB-(op(upB)-i):
  fi:

  if i=1 then
    Slist:=[ [xi,1,upB], op(Slist)]:
  else
    Slist:=[ [xi,xset[i-1]+1,upB], op(Slist)]:
  fi:

od:


#[mul(z[op(xi)]^xi, xi in usedx), op(Slist)]:

[mul(z[op(xi)]^xi, xi in xset), op(Slist)]:


end:




##########################
# Extract x[i]'s utilized in SCHEME's RHS (says set R)
# - Use ExtractNewVars1 for "x[i] expressions"
# - Extract actual x[i]'s from expressions

# Create "Slist" for function SimplifySums

# For each summand in "simplified output", factor out remaining z[i]^x[i]'s
# 
# Which yi's have `=` in cond?
# - if not -> 1 entry in recursion
# - if so, find all zj^xj with corresponding xj's
# - divide out by zj^(xj-k) for appropriate k
# 
##########################






# SimplifySums(Slist): inputs a list Slist representing a nested summation (of k sums)
# in the following format:
# 
# Slist = [summand, [sum1_var, sum1_LB, sum1_UB], .., [sumk_var, sumk_LB, sumk_UB] ]
#
# It outputs the summand simplified via Maple's "sum" command.
#
SimplifySums:=proc(Slist) local L1:

if nops(Slist)=1 then
  RETURN(Slist[1]):
else
  L1:=Slist[2]:
  RETURN(SimplifySums([sum(Slist[1], L1[1]=L1[2]..L1[3]), op(3..nops(Slist),Slist)])):
fi:

end:






# Ex: S:={x[1]=i[2], x[2]=i[3], i[1]<i[2], i[2]<i[3]}:
#     ExtractNewVars2(S,[i[1],i[2],i[3]])
#   will output {x[1], x[2]}
ExtractNewVars2:=proc(S,cvars)
local c,newvars,newvars2,c1:

newvars:={}:

for c in S do
  if type(c,`=`) then
    newvars:=newvars union convert(c,set):
  fi:
od:

newvars minus convert(cvars,set):

end:





#-----------------------------------------------------------------------
#  Support/Miscellaneous Procedures
#-----------------------------------------------------------------------

CompactSets:=proc(S)
local Sout,S1,Stmp:

Sout:={}:

for S1 in S do
  Stmp:=TrivEqSets(S1):
  Sout:=Sout union {Stmp[1]}:
od:

Sout:

end:



TrivEqSets:=proc(S)
local Sf,Sc,Sfc,n,p,pf:

Sf:={}:
Sc:={}:
Sfc:={}:

for p in S do
  n:=nops(p):
  Sf:=Sf union {[seq(p[n-i],i=0..n-1)]}:
  Sc:=Sc union {[seq(n-p[i]+1,i=1..n)]}:
  pf:=[seq(p[n-i],i=0..n-1)]:
  Sfc:=Sfc union {[seq(n-pf[i]+1,i=1..n)]}:
od:

{S,Sf,Sc,Sfc}:

end:



# Removes trivial equivalences from set S of single patterns
RemoveTrivEq:=proc(S)
local p,pc,pf,pfc,i,S1,n,Sp:

S1:=S:

for p in S1 do
  n:=nops(p):
  pf:=[seq(p[n-i],i=0..n-1)]:
  pc:=[seq(n-p[i]+1,i=1..n)]:
  pfc:=[seq(n-pf[i]+1,i=1..n)]:
  Sp:={p,pf,pc,pfc}:
  S1:=S1 minus Sp union {Sp[1]}:
od:

S1:

end:



# Checks that B is a set of bad patterns in the proper format
CheckSetB:=proc(B)
local p,n:

if not type(B,set) then
  ERROR(`Input B must be a SET of patterns.`):
fi:

if B={} then
  ERROR(`Input B must be a non-empty set.`):
fi:

for p in B do
  if not type(p,list) then
    ERROR(`All elements of B must be lists. B contained non-list %1.`, p):
  fi:

  n:=nops(p):

  if convert(p,set)<>{seq(i,i=1..n)} then
    ERROR(`Invalid pattern %1.`, p):
  fi:

  if n<3 then
    ERROR(`Patterns must be of length at least 3. Set B contained %1.`,p):
  fi:

od:

end:



#-----------------------------------------------------------------------
#  Print Procedures
#-----------------------------------------------------------------------



# Verbose output for CAV(B,N)
CAVprint:=proc(B,S,L)
printf(`Consecutive Avoidance in Permutations (using CAV) with the Pattern Set %a\n\n\n`, B):
printf(`The following scheme was utilized to count permutations that consecutively avoid %a:\n\n%a\n\n`, B, S):
printf(`The number of permutations of length i (1<=i<=%d) consecutively avoiding %a `, nops(L),B):
printf(`is given by the sequence:\n\n%a\n\n\n`, L):
end:



# Verbose output for CAVt(B,N,t)
CAVtprint:=proc(B,S,L,t)
printf(`Consecutive Avoidance in Permutations (using CAVt) with the Pattern Set %a\n\n\n`, B):
printf(`NOTE: This version keeps track of errors in each permutation as well.\n\n`):
printf(`The following scheme was utilized to count permutations that consecutively avoid %a:\n\n%a\n\n`, B, S):
printf(`The number of permutations of length i (1<=i<=%d) with j instances of patterns in %a `, nops(L),B):
printf(`is given by the coefficient of %a^j in the i-th term in the sequence:\n\n%a\n\n\n`, t, L):
end:



# Verbose output for CAVF(B,N)
CAVFprint:=proc(B,S,L)
printf(`Consecutive Avoidance in Permutations (using CAVF) with the Pattern Set %a\n\n\n`, B):
printf(`The following scheme was utilized to count permutations that consecutively avoid %a:\n\n%a\n\n`, B, S):
printf(`The number of permutations of length i (1<=i<=%d) consecutively avoiding %a `, nops(L),B):
printf(`is given by the sequence:\n\n%a\n\n\n`, L):
end:



# Verbose output for CAVFt(B,N,t)
CAVFtprint:=proc(B,S,L,t)
printf(`Consecutive Avoidance in Permutations (using CAVFt) with the Pattern Set %a\n\n\n`, B):
printf(`NOTE: This version keeps track of errors in each permutation as well.\n\n`):
printf(`The following scheme was utilized to count permutations that consecutively avoid %a:\n\n%a\n\n`, B, S):
printf(`The number of permutations of length i (1<=i<=%d) with j instances of patterns in %a `, nops(L),B):
printf(`is given by the coefficient of %a^j in the i-th term in the sequence:\n\n%a\n\n\n`, t, L):
end:



# Verbose output for CAVT(B,N)
CAVTprint:=proc(B,S,FE,L) global F:

printf(`Consecutive Avoidance in Permutations (using CAVT) with the Pattern Set %a\n\n\n`, B):
printf(`CAVT created the following initial scheme for consecutively avoiding %a:\n\n%a\n\n\n`, B, S):
printf(`Using this, the Cluster Tail GF functional equation was created. `):
printf(`It is encoded as:\n\n%a\n\n\n`,FE):
printf(`A slightly more human friendly form can be written as:`):
printf(`\n\n%a\n\n\n`, TailFE2(FE,[seq(F[i],i=1..NumFuncInScheme(S))])):
printf(`This Cluster Tail GF functional equation was used to count permutations avoiding %a.\n`, B):
printf(`The number of permutations of length i (1<=i<=%d) consecutively avoiding %a `, nops(L),B):
printf(`is given by the sequence:\n\n%a\n\n\n`, L):

end:



# Verbose output for CAVTt(B,N,t)
CAVTtprint:=proc(B,S,FE,L,t) global F:

printf(`Consecutive Avoidance in Permutations (using CAVTt) with the Pattern Set %a\n\n\n`, B):
printf(`NOTE: This version keeps track of errors in each permutation as well.\n\n`):
printf(`CAVTt created the following initial scheme for consecutively avoiding %a:\n\n%a\n\n\n`, B, S):
printf(`Using this, the Cluster Tail GF functional equation was created. `):
printf(`It is encoded as:\n\n%a\n\n\n`,FE):
printf(`A slightly more human friendly form can be written as:`):
printf(`\n\n%a\n\n\n`, TailFE2(FE,[seq(F[i],i=1..NumFuncInScheme(S))])):
printf(`This Cluster Tail GF functional equation was used to count permutations `):
printf(`with various numbers of occurrences/errors from %a.\n`, B):
printf(`The number of permutations of length i (1<=i<=%d) with j instances of patterns in %a `, nops(L),B):
printf(`is given by the coefficient of %a^j in the i-th term in the sequence:\n\n%a\n\n\n`, t, L):

end:



# Verbose output for MapWE(n)
MapWEprint:=proc(n,S) local S1:

printf(`(Non-trivial) Strong Wilf-Equivalences for Length %d Patterns \n\n`, n):

printf(`NOTE: We will disregard the trivial equivalences of reversals and complementations.\n\n`):
printf(`Based off of the OverlapMaps equivalence, we have the following Strong Wilf-Equivalences: \n\n`):

for S1 in S do
  if nops(S1)=1 then
    printf(`Did not find a Strong c-Wilf-Equivalence for %a.\n\n`, op(S1)):
  elif nops(S1)>1 then
    printf(`The following single patterns are Strongly c-Wilf-Equivalent: %a \n`, S1):
    printf(`because they have the same OverlapMap: %a \n\n`, OverlapMaps(S1[1],S1[1])):
  fi:
od:

end:



# Verbose output for RhoApprox(p,N)
RhoApproxprint:=proc(L,N,p)

printf(`Approximate Rho Value for the Single Pattern %a \n\n`, p):

printf(`NOTES: For a single pattern, we assume A(n) ~ %a*(%a^n)*n!\n`, gamma, rho):
printf(`The number of decimal places is currently set to %d.\n\n\n`, Digits):
printf(`The single pattern %a has an approximate rho value: %a\n`, p, L[1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[2]):

end:



# Verbose output for AsymApprox(p,N)
AsymApproxprint:=proc(L,N,p)

printf(`Approximate Asymptotics for the Single Pattern %a \n\n`, p):

printf(`NOTES: For a single pattern, we assume A(n) ~ %a*(%a^n)*n!\n`, gamma, rho):
printf(`The number of decimal places is currently set to %d.\n\n\n`, Digits):
printf(`The single pattern %a has an approximate rho value: %a\n`, p, L[1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n`, N-1, N, L[2]):
printf(`It also has an approximate gamma value: %a\n`, L[3]):
printf(`The error between the n=%d and n=%d approximated gamma terms is %a.\n\n\n`, N-1, N, L[4]):

end:



# Verbose output for RhoApproxRank(n,N)
RhoApproxRankprint:=proc(L,N)
local i,m,k:

m:=nops(L):
k:=nops(L[1][3][1]):

printf(`Ranking of Approximate Rho Values for Length %d Patterns \n\n`, k):

printf(`NOTES: For a single pattern, we assume A(n) ~ %a*(%a^n)*n!\n`, gamma, rho):
printf(`The number of decimal places is currently set to %d.\n\n\n`, Digits):

printf(`The set of single patterns %a produces the largest rho value: %a\n`, L[m][3], L[m][1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[m][2]):

if m>2 then
  for i from 2 to m-1 do
    printf(`The set of single patterns %a produces the next smallest rho value: %a\n`, L[m+1-i][3], L[m+1-i][1]):
    printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[m+1-i][2]):
  od:
fi:

printf(`The set of single patterns %a produces the smallest rho value: %a\n`, L[1][3], L[1][1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[1][2]):


end:



# Verbose output for AsymApproxRank(n,N)
AsymApproxRankprint:=proc(L,N)
local i,m,k:

m:=nops(L):
k:=nops(L[1][5][1]):

printf(`Approximate Asymptotics (Ranked by Rho) for Length %d Patterns \n\n`, k):

printf(`NOTES: For a single pattern, we assume A(n) ~ %a*(%a^n)*n!\n`, gamma, rho):
printf(`The number of decimal places is currently set to %d.\n\n\n`, Digits):

printf(`The set of single patterns %a produces the largest rho value: %a\n`, L[m][5], L[m][1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n`, N-1, N, L[m][2]):
printf(`The coefficient gamma is approximately: %a\n`, L[m][3]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[m][4]):

if m>2 then
  for i from 2 to m-1 do
    printf(`The set of single patterns %a produces the largest rho value: %a\n`, L[m+1-i][5], L[m+1-i][1]):
    printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n`, N-1, N, L[m+1-i][2]):
    printf(`The coefficient gamma is approximately: %a\n`, L[m+1-i][3]):
    printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[m+1-i][4]):
  od:
fi:

printf(`The set of single patterns %a produces the smallest rho value: %a\n`, L[1][5], L[1][1]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n`, N-1, N, L[1][2]):
printf(`The coefficient gamma is approximately: %a\n`, L[1][3]):
printf(`The error between the n=%d and n=%d approximated rho terms is %a.\n\n\n`, N-1, N, L[1][4]):

end:





WilfEqmprint:=proc(L,V,n,N,m)
local Llen,i,Li:

if m>1 then
  printf(`(Non-trivial) Strong c-Wilf-Equivalences for Sets of %d Patterns of Length %d\n\n`, m, n):
else
  printf(`(Non-trivial) Strong c-Wilf-Equivalences for Single Patterns of Length %d\n\n`, n):
fi:

printf(`NOTE: We will disregard the trivial equivalences of reversals and complementations.\n\n\n`):


if V ={} then
  printf(`We were able to attain a TOTAL CLASSIFICATION of Wilf-Equivalences!\n`):
  printf(`The Strong c-Wilf-Equivalences are below `):
  printf(`(and there were no Wilf-Equivalences that were not strong):\n\n\n`):
fi:

Llen:=nops(L):

for i from 1 to Llen do
  Li:=L[Llen+1-i]:
  
  if m=1 then
    if nops(Li[2])>1 then
      printf(`The following single patterns are Strong c-Wilf-Equivalent by OverlapMaps: %q\n`, op(Li[2])):
      printf(`They have the same OverlapMaps: %a\n`, OverlapMaps(Li[2][1],Li[2][1])):
      printf(`For each single pattern, A(%d)=%a\n\n\n`, N, Li[1]):
    else
      printf(`No Strong c-Wilf-Equivalences could be proven for the single pattern: %a\n`, Li[2][1]):
      printf(`It has OverlapMaps: %a\n`, OverlapMaps(Li[2][1],Li[2][1])):
      printf(`For this single pattern, A(%d)=%a\n\n\n`,N,Li[1]):
    fi:
  elif m>1 then
    if nops(Li[2])>1 then
      printf(`The following %d-sets of patterns are Strong c-Wilf-Equivalent by OverlapMaps: %q\n`, m, op(Li[2])):
      printf(`They have the same OverlapMaps: %a\n`, OverlapMaps(Li[2][1],Li[2][1])):
      printf(`For each set of patterns, A(%d)=%a\n\n\n`, N, Li[1]):
    else
      printf(`No Strong c-Wilf-Equivalences could be proven for the set of patterns: %a\n`, Li[2][1]):
      printf(`It has OverlapMaps: %a\n`, OverlapMaps(Li[2][1],Li[2][1])):
      printf(`For this set of patterns, A(%d)=%a\n\n\n`,N,Li[1]):
    fi:


  fi:

od:


if V<>{} then
  if m=1 then
    printf(`The following single patterns MIGHT be Wilf-Equivalent, `):
    printf(`but this could not be proven or disproven.\n`):
    printf(`Trying a larger N value may help in disproving some of the below ambiguities.\n\n`):

    for i from 1 to nops(V) do
      printf(`The following (proven) single pattern Strong c-Wilf-Equivalence classes `):
      printf(`MIGHT also be Wilf-Equivalent: %q\n\n`, op(V[i])):
    od:

  else
    printf(`The following sets of patterns MIGHT be Wilf-Equivalent, `):
    printf(`but this could not be proven or disproven.\n`):
    printf(`Trying a larger N value may help in disproving some of the below ambiguities.\n\n`):

    for i from 1 to nops(V) do
      printf(`MIGHT be Wilf-Equivalent: %q\n\n`, op(V[i])):
    od:
  fi:

  printf(`\n`):

fi:


end:





#-----------------------------------------------------------------------
#  Testing Code
#-----------------------------------------------------------------------


Test132:=proc(n)
local F,C,k,x1,x2,x3,N,j:
option remember:

if n<3 then
  RETURN(n!):
fi:

F:=(k,x1,x2,x3)->piecewise(k=3,-1,k<3,0,-(sum(sum(F(k-2,y1,x1,y3), y1=1..x1-1), y3=1+x1..k-2))):

C:=N->sum(sum(sum(F(N,x1,x2,x3), x3=x2+1..N), x2=x1+1..N-1), x1=1..N-2):

n*Test123(n-1)+add(binomial(n,j)*C(j)*Test123(n-j),j=3..n):

end:




TailGFNaive:=proc(Sc,i,k,N,z)
local j,n,CSc,xn,S,N1,TGF:

n:=NumFuncInScheme(Sc):
CSc:=ConvertScheme2(Sc,n):

xn:=nops(CSc[i][2][2][3]):

S:={op(choose(N,xn))}:
add(ComputeFiNx(CSc,i,N,xvect,k)*mul(z[j]^xvect[j],j=1..xn), xvect in S):

end:




TestRHS132:=proc(Sc,i,k,N,z)
local n,CSc,TGF,x1,x2,x3,y1,y3:

n:=NumFuncInScheme(Sc):
CSc:=ConvertScheme2(Sc,n):

TGF:=0:

for x1 from 1 to N-2 do
  for x2 from x1+1 to N-1 do
    for x3 from x2+1 to N do
      for y1 from 1 to x1-1 do
        for y3 from x1+1 to N-2 do
          TGF:=TGF+ComputeFiNx(CSc,i,N-2,[y1,x1,y3],k)*z[1]^x1*z[2]^x2*z[3]^x3:
        od:
      od:
    od:
  od:
od:

-TGF:

end:





WilfEq:=proc(n,N)
local Sn,Sn1,L,p,Sc,FE,S,Li,i,Lprev,violators,verb:
global k,x,y,z:

verb:=false:

if nops([args])=3 then
  if type(args[3],boolean) then
    verb:=args[3]:
  else
    ERROR(`Optional 3rd argument (verbose) must be of type boolean.`):
  fi:
fi:


Sn:=convert(permute(n),set):
Sn1:=RemoveTrivEq(Sn):

L:={}:

for p in Sn1 do
  Sc:=SCHEME(k,{p},x,y,0):
  FE:=ScToTailSc(Sc,z):
  L:=L union {[CAVOIDtail(FE,N,k), p]}:
od:

L:=convert(L,list):

S:={ [L[1][1], {L[1][2]}]}:

for i from 2 to nops(L) do
  Li:=L[i]:
  Lprev:=S[nops(S)]:
  
  if Li[1]=Lprev[1] then
    S:=S minus {Lprev}:
    S:=S union {[Lprev[1], Lprev[2] union {Li[2]}]}:
  else
    S:=S union {[Li[1], {Li[2]}]}:
  fi:
od:


violators:={}:

for Li in S do
  if nops(Li[2])>1 then
    if SameScheme(k,Li[2],x,y)=false then
      violators:=violators union {Li[2]}:
    fi:
  fi:
od:

if verb then
  print(` Under construction. Please use VERBOSE=false. `):
else
  if violators={} then
    RETURN(S,true):
  else
    RETURN(S,false,violators):
  fi:
fi:

end:





SameScheme:=proc(k,S1,x,y)
local Sc1,Sc2,i:

if nops(S1)=1 then
  RETURN(true):
fi:


Sc1:=SCHEME(k,{S1[1]},x,y,0):

for i from 2 to nops(S1) do
  Sc2:=SCHEME(k,{S1[i]},x,y,0):
  if evalb(Sc1<>Sc2) then
    RETURN(false):
  fi:
od:

RETURN(true):

end:





####End Brian Nakamura's program#########


#ECGFbrian(FP,K): The first K terms of the empirical 
#cluster exponential generating function. For example, try:
#ECGFbrian({[1,3,2]},20);
ECGFbrian:=proc(FP,K) local z,gu,f,i:
gu:=CAVT(FP,K):
f:=1/(1+add(gu[i]*z^i/i!,i=1..K)):
f:=taylor(f,z=0,K+1):
if coeff(f,z,0)<>1 or coeff(f,z,1)<>-1 then
 RETURN(FAIL):
fi:

[0,seq(-coeff(f,z,i)*i!,i=2..K)]:
end:

#AnPlusBrian(N,FP,N1,C,n,E): The first N1 terms of the sequence
#enumerating the number of permutations of length n
#avoiding the set of patterns FP using the first N terms
#and conjecturing a linear recurrence of complexity <=C
#satisfied by
#the coefficients of the cluster generating function
#For example, try:
#AnPlusBrian(20,{[1,2,3]},40,6,n,E);
AnPlusBrian:=proc(N,FP,N1,C,n,E) local ope,gu,f,z,i:
gu:=ECGFbrian(FP,N):
ope:=Findrec(gu,n,E,C):

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

if SeqFromRec(ope,n,E,[op(1..degree(ope,E),gu)],nops(gu))<>gu then
 print(`Something is wrong`):
 RETURN(FAIL):
fi:


gu:=SeqFromRec(ope,n,E,[op(1..degree(ope,E),gu)],N1):


f:=1-z+add(gu[i]*z^i/i!,i=1..N1):

f:=taylor(1/f,z=0,N1+1):
[seq(coeff(f,z,i)*i!,i=1..N1)],ope:


end:




#AnPam(a,m,N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the number of
#occurrences of the  consecutive pattern 1...a tau (a+1) in S_{m+1}
#for any tau.
#For example, try:
#AnPam(1,2,10,t);
AnPam:=proc(a,m,N,t) local gu,x,i,j,k:
gu:=1-x-add((t-1)^k*x^(k*m+1)/(k*m+1)!*mul(binomial(j*m-a,m-a),j=1..k),
k=1..trunc((N-1)/m)+2):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:


#DKfkp(k,p,m,a,b): implementing the recursion (5) of 
#Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance. 
#It is the number of k-clusters for any consecutive pattern
#that starts with a+1 and ends at b+1 such that the cluster
#starts at p+1
#For example, try:
#DKfkp(2,6,1,3,2);
DKfkp:=proc(k,p,m,a,b) local q:
option remember:
if a=b then
 RETURN(0):
fi:
if a>b then
RETURN(DKfkp(k,p,m,b,a)):
fi:

if k=1 then
 if p=a then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:

add(DKfkp(k-1,q-b,m,a,b)*
binomial(p,a)*binomial(k*m-q,m-b)*binomial(q-p-1,b-a-1),q=p+1..m*k):
end:




#DKfk(k,m,a,b): The number of k-clusters of length m*k+1
#for a pattern tau of length m+1 that starts with
#a+1 and ends with b+1 in
#Anton Khoroshkin's paper "Anick-type resolutions and consecutive pattern
#avoidance. 
#For example, try:
#DKfk(2,3,1,2);
DKfk:=proc(k,m,a,b) local p:
add(DKfkp(k,p,m,a,b),p=0..k*m):
end:



#AnPDK(a,b,m,N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the weight enumerator for
#occurrences of the  consecutive pattern tau of length m+1
#with exactly one overlap and
#tau[1]=a+1 and tau[m+1]=b+1 following the formula
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance."
#AnPDK(1,2,3,10,t);
AnPDK:=proc(a,b,m,N,t) local gu,x,i,k:
gu:=1-x-add((t-1)^k*x^(k*m+1)/(k*m+1)!*
DKfk(k,m,a,b),k=1..trunc((N-1)/m)+2):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:


#DK1324cnl(n,l): Formula (6) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK1324cnl(10,3);
DK1324cnl:=proc(n,l) local k:
option remember:
if n=1 then
 if l=0 then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:

if n=2 or n=3 then
 RETURN(0):
fi:

add(binomial(2*k,k)/(k+1)*DK1324cnl(n-2*k-1,l-k),k=1..trunc((n-2)/2)):
end:


#DK1324cn(n,t): Summing the DK1324cnl(n,l) of Formula (6) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK1324cn(10,t);
DK1324cn:=proc(n,t) local l:
expand(add(DK1324cnl(n,l)*(t-1)^l,l=1..n)):
end:


#AnP1324(N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the weight enumerator for
#occurrences of the  consecutive pattern 1324
#taken from Theorem 4 in
#Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For examle, try:
#AnP1324(10,t);
AnP1324:=proc(N,t) local gu,x,i,n:
gu:=1-x-add(x^n/n!*DK1324cn(n,t),n=1..N):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:





#DK1423cnl(n,l): Formula (8) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK1423cnl(10,3);
DK1423cnl:=proc(n,l) local k:
option remember:
if n=1 then
 if l=0 then
  RETURN(1):
 else
  RETURN(0):
 fi:
fi:

if n=2 or n=3 then
 RETURN(0):
fi:

add(binomial(n-k-2,k)*DK1423cnl(n-2*k-1,l-k),k=1..trunc((n-2)/2)):
end:


#DK1423cn(n,t): Summing the DK1324cnl(n,l) of Formula (8) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK1423cn(10,t);
DK1423cn:=proc(n,t) local l:
expand(add(DK1423cnl(n,l)*(t-1)^l,l=1..n)):
end:


#AnP1423(N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the weight enumerator for
#occurrences of the  consecutive pattern 1324
#taken from Theorem 5 in
#Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For examle, try:
#AnP1423(10,t);
AnP1423:=proc(N,t) local gu,x,i,n:
gu:=1-x-add(x^n/n!*DK1423cn(n,t),n=1..N):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:


#AnPincr(a,N): The first N terms of the sequence
#enumerating the permutations of length n
#avoiding the consecutive pattern [1,2,3,.., a]
#using the formula of Sergi Elizalde and Marc Noy. For exaple, try:
#AnPincr(3,20);
AnPincr:=proc(a,N) local gu,x,i,k:
gu:=add(x^(k*a)/(k*a)!,k=0..trunc(N/a)+1)-
add(x^(k*a+1)/(k*a+1)!,k=0..trunc(N/a)+1):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:




#DK2413cnlpq(n,l,p,q): Formula (14) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK2413cnlpq(10,3,1,2);
DK2413cnlpq:=proc(n,l,p,q) local r,s:
option remember:
if n=2 or n=3 then
 RETURN(0):
fi:
if n=4 then

 if l=1 and p=2 and q=4 then
    RETURN(1):
  else
    RETURN(0):
 fi:
fi:

add(add(DK2413cnlpq(n-2,l-1,r,s-1),r=1..p-1),s=p+1..q-1)+
(p-1)*add(add(DK2413cnlpq(n-3,l-1,r-1,s-1),r=p+1..s-1),s=p+2..q-1)+
(p-1)*add(add(DK2413cnlpq(n-3,l-1,r-1,s-2),r=p+1..q-1),s=q+1..n):
end:

#DK2413cnl(n,l): used in Theorem 7 of
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK2413cnl(10,3);
DK2413cnl:=proc(n,l) local p,q:
add(add(DK2413cnlpq(n,l,p,q),p=2..q-2),q=3..n-1): 
end:


#DK2413cn(n,t): Summing the DK2413cnl(n,l) of Formula (14) 
#in Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For example, try:
#DK2413cn(10,t);
DK2413cn:=proc(n,t) local l:
expand(add(DK2413cnl(n,l)*(t-1)^l,l=1..n)):
end:


#AnP2413(N,t): The first N terms of the sequence
#enumerating the weight-enumerators of permutations of length n
#according to the weight enumerator for
#occurrences of the  consecutive pattern 2413
#taken from Theorem 7 in
#Vladimir Dotsenko and Anton Khoroshkin's 
#paper "Anick-type resolutions and consecutive pattern
#avoidance." For examle, try:
#AnP2413(10,t);
AnP2413:=proc(N,t) local gu,x,i,n:
print(`Wrong!`):
gu:=1-x-add(x^n/n!*DK2413cn(n,t),n=1..N):
gu:=1/gu:
gu:=taylor(gu,x=0,N+1):
[seq(expand(i!*coeff(gu,x,i)),i=1..N)]:
end:






#khaverim(p): all the friends of the pattern p.
#For example, try: khaverim([1,2,3,4]);
khaverim:=proc(p) local i,n:
n:=nops(p):
{p,[seq(p[n+1-i],i=1..nops(p))],
[seq(n+1-p[i],i=1..n)],
[seq(n+1-p[n+1-i],i=1..nops(p))]}:
end:

#Khaverim(S): all the friends of the set of patterns S.
#of the same length
#For example, try: Khaverim({[1,2,3,4]});
Khaverim:=proc(S) local i,n,p:
n:={seq(nops(p), p in S)}:
if nops(n)<>1 then
print(`All patterns in `, S, `should be the same length`):
RETURN(FAIL):
fi:
n:=n[1]:
{S,
{seq([seq(p[n+1-i],i=1..nops(p))], p in S)},
{seq([seq(n+1-p[i],i=1..n)],p in S)},
{seq([seq(n+1-p[n+1-i],i=1..nops(p))],p in S)}
}:
end:







#Theorem(FP): The theorem that enables fast enumeration
#of permutations avoiding the set of patterns FP
#of the same length.
#For example, try:
#Theorem({[1,2,3]});
Theorem:=proc(FP) local x,z,pol,i,lu,pi,S,k,i1,g,mu:
option remember:
if FP={} then
 print(`Set of patterns must be non-empty`):
fi:

S:={seq(nops(pi), pi in FP)}:

if nops(S)>1 then
 print(`the case where all patterns do not have the same length is  not implemented`):
 RETURN(FAIL):
fi:

k:=S[1]:

lu:=Um(x,z,FP):

mu:=ApplyUmbralMatrix(lu,[seq(g[i](seq(x[i1],i1=1..k-1),z),i=1..(k-1)!)],
[seq(x[i1],i1=1..k-1),z]):

pol:=Init(k-1,x,z):

print(`The (ordinary) generating function`):
print(`whose coeff. of z^n tells you the number of permutations of`):
print(`length n avoiding the set of patterns`, FP):
print(` equals `):
print(Sum(g[i](1$(k-1),z),i=1..(k-1)!)):
print(`where the vector of formal power series of length`, (k-1)!):
print(`in the variables`, seq(x[i],i=1..(k-1)!),z):
print(`is the unique vector that`):
print(`satisfies the following system of functional equations`):

for i from 1 to (k-1)! do
print(g[i](seq(x[i1],i1=1..k-1),z)=pol[i]+mu[i]):
od:

end:









#####stuff with general t[q]########
#UmNicet(x,z,k,t,f,g): A nice rendition of 
#an Umbral matrix operator with the input function
#being called f[i] and the output function called g[i]
#for all general patterns of length k
#For example, try:
#UmNicet(x,z,3,t,f,g);
UmNicet:=proc(x,z,k,t,f,g) local gu,i1,j1,lu,mu,lu1:

gu:=UmSt(x,z,k,t):

for i1 from 1 to nops(gu) do
mu:=0:
 for j1 from 1 to nops(gu[i1]) do
 lu:=gu[i1][j1]:
 mu:=mu+add(lu1[1]*f[j1](op(lu1[2])),lu1 in lu):
 od:
 print(g[i1]=mu):
od:

end:


#UmSt(x,z,k,t): A simplified version of the global
#umbral evolution matrix operator, for enumerating
#permutations avoiding the patterns (of the same length) in FP
#in terms of an umbral matrix and catalytic variables
#x[1],x[2],..x[k] and  the variable z (for the size of the permutation)
#for general t
#For example, try:
#UmSt(x,z,3,t);
UmSt:=proc(x,z,k,t) local gu,i1,j1:
 gu:=UAtV(k-1,x,z,t):
[seq([seq(SimpU(gu[i1][j1]),j1=1..nops(gu[i1]))],i1=1..nops(gu))]:
end:


Putz1:=proc(ju,z): 
[z*ju[1],ju[2],[op(1..nops(ju[3])-1,ju[3]),ju[3][nops(ju[3])]*z]]:end:

Putz:=proc(S,z) local ju:{seq(Putz1(ju,z),ju in S)}:end:
#Theoremt(k,t): The theorem that enables fast 
#weighted enumeration according to the weight
#keep track of all consecutive patterns of length k
#t is the variable for indexing the consecutive patters.
#For example, try:
#Theoremt(3,t);
Theoremt:=proc(k,t) local x,z,pol,i,lu,i1,g,mu,j1:
option remember:



lu:=UAtV(k-1,x,z,t):

lu:=
[seq(
[seq(Putz(lu[i1][j1],z),j1=1..nops(lu[i1]))],i1=1..nops(lu))]:

lu:=CleanUpU2(lu):
mu:=ApplyUmbralMatrix(lu,[seq(g[i](seq(x[i1],i1=1..k-1),z),i=1..(k-1)!)],
[seq(x[i1],i1=1..k-1),z]):


pol:=Init(k-1,x,z):

print(`The (ordinary) generating function`):
print(`whose coeff. of z^n is the weight-enumerator of permutations`):
print(`of length n where the weight of a permutation pi`):
print(` is the product of t[q] `):
print(`to the power of the number of times the consecutive pattern`):
print(`q shows up in pi, for all patterns q of length`, k):
print(` equals `):
print(Sum(g[i](1$(k-1),z),i=1..(k-1)!)):
print(`where the vector of formal power series of length`, (k-1)!):
print(`in the variables`, seq(x[i],i=1..(k-1)!),z):
print(`is the unique vector that`):
print(`satisfies the following system of functional equations`):

for i from 1 to (k-1)! do
print(g[i](seq(x[i1],i1=1..k-1),z)=pol[i]+mu[i]):
od:

end:

#####end stuff with general t[q]########