######################################################################
##ELIZALDE: Save this file as  ELIZALDE                              #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read ELIZALDE<Enter>                               #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#zeilberg at math dot rutgers dot edu                                #
######################################################################
 
#Created: Nov. 2010-Jan. 2011
 
print(`Created:  Nov. 2010-Jan. 2011`):
print(`Version of Dec. 29, 2010`):
print(` This is ELIZALDE `):
print(`A Maple package that accompanies the article `):
print(`"Automatic Generation of Theorems and Proofs on `):
print(`Enumerating Consecutive-Wilf classes"`):
print(`by Andrew Baxter, Brian Nakamura and Doron Zeilberger.`):
print(`available from the authors' websites`):
print(`The most current version of this Mape package and of the article`):
print(`are available from Zeilberger's 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 Umbral procedures type ezraU();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

 print(`For a list of the DIRECT`):
print(`Umbral procedures type ezraD();, 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(``):

with(combinat):


ezraD:=proc()

if args=NULL then
 print(` The umbral procedures are: ApplyUD,ApplyUD1, ApplyUD11 `):
 print(` ClusterSeqDir,  UpScheme `):
else
ezra(args):
fi:

end:

ezraU:=proc()

if args=NULL then
 print(` The umbral procedures are: `):
 print(` ApplyUmbra,ApplyUmbraMat, GenSeq, GroupU,`):
 print(` PtoU, PtoBU, `):
 print(` schumT, schumT1, SpellOutUmbra, SPtoUM, TemplateUP, ToUmbra `):

else
ezra(args):
fi:

end:



ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: `):
 print(` Actualize, AsyMoments, CleaunUp, ClusterSeqFg `):
 print(` ClusterSeqS, FastestFriend, `):
 print(` FEx1, FillIn, Foll, `):
 print(` gTemplate, Hakol, khaverim, Moments`):
 print(` Pred, PtoU, redu, RedVar, `):
 print(` ReduceU ,ReduceU1, reps`):
 print(` Scheme,SeqFromDE, SeqS, SeqT, `):
 print(`SeqTu, SeqTuS, , SipurF, Template, Val1 `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: `):
 print(` AllSeqs, AllSeqsF,`):
 print(` ClusterSeq, ClusterSeqF, Hochakha`):
 print(`   Mishpat, PtoU, PtoBU, SeferP, SeferT, Seq, SeqF, SeqW, Sipur `):
 print(` UtoD `):


elif nops([args])=1 and op(1,[args])=Actualize then
print(`Actualize(Temp,i,i0,n0): given a template in terms of symbols`):
print(`i[1],i[2],i[3],... and  a specfic numerical increasing vector`):
print(`of integers i0, and a pos. integer n0, finds all the`):
print(`possible vectors of integers with those specifications.`):
print(`For example, try:`):
print(`Actualize([0,0,i[2]-1],i,[5,7,8],11);`):


elif nops([args])=1 and op(1,[args])=AllSeqs then
print(`AllSeqs(k,N): All possible sequences enumerating`):
print(`single 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,30);`):

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

elif nops([args])=1 and op(1,[args])=ApplyUD then
print(`ApplyUD(U,poly, xvar): Given an upscheme U`):
print(`and a polynomial, poly, in the list of variables xvar,`):
print(`applies the implied umbral operator U to the polynomial pol.`):
print(`For example, try:`):
print(`ApplyUD1([{[j1,0,j2,n+2]},[j1,j2,j3,n]],3*x1*x2^5*x3^11*z^23,[x1,x2,x3,z]);`):


elif nops([args])=1 and op(1,[args])=ApplyUD1 then
print(`ApplyUD1(U,mono, xvar): Given an upscheme U`):
print(`and a monomial in the list of variables xvar,`):
print(`applies the implied umbral operator U to the monomial mono`):
print(`For example, try:`):
print(`ApplyUD1([{[j1,0,j2,n+2]},[j1,j2,j3,n]],3*x1*x2^5*x3^11*z^23,[x1,x2,x3,z]);`):


elif nops([args])=1 and op(1,[args])=ApplyUD11 then
print(`ApplyUD11(T,mono,avar, xvar): Given a template T `):
print(`a monomial mono, a list of discrete variables avar`):
print(`and a list of continuous variables xvar`):
print(`(with nops(T)=nops(avar)=nops(xvar))`):
print(`Applies the implied umbral operator given by T.`):
print(`For example, try:`):
print(`ApplyUD11([j1,0,j2,n+2],3*x1*x2^5*x3^11*z^23,[j1,j2,j3,n],[x1,x2,x3,z]);`):

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])=ApplyUmbraMat then
print(`ApplyUmbraMat(Umbra1Mat,f,ylist,A): Given a list of`):
print(`umbral operators Umbra1Mat phrased in terms of A[1],A[2]`):
print(`and a linear combination of A[1],A[2], ... with coefficients`):
print(`that are polynomials in ylist applies Umbra1Mat to f`):
print(`For example try:`):
print(`ApplyUmbraMat(SPtoUM([[1,3,2]],t,x,z,A),x[1]*x[2]^2*x[3]^3*z^3*A[1],[x[1],x[2],x[3],z]],A);`):



elif nops([args])=1 and op(1,[args])=AsyMoments then
print(`AsyMoments(p,s,n): The  conjectured formula, to order s, for the`):
print(`normalized asymptotic moments (both even and odd)`):
print(`of the random variable #number of occurrences of pattern`):
print(`p in uniformly-at-random n-permutation`):
print(`For example try: `):
print(`AsyMoments([1,2,3],1,n);`):

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

elif nops([args])=1 and op(1,[args])=ClusterSeq then
print(`ClusterSeq(p,t,N): Given a consecutive pattern`):
print(`uses the Umbral scheme to find the first N terms`):
print(`of the cluster-series for counting permutations`):
print(`according to the number of occurrences of p.`):
print(`For example, try:`):
print(`ClusterSeq([1,2,3],t,10);`):

elif nops([args])=1 and op(1,[args])=ClusterSeqFg then
print(`ClusterSeqFg(p,t,N,x): Like ClusterSeqF(p,t,N)`):
print(`but retaining the catalytic variables`):
print(`For example, try:`):
print(`ClusterSeqFg([1,2,3],t,10,x);`):


elif nops([args])=1 and op(1,[args])=ClusterSeqF then
print(`ClusterSeqF(p,t,N): a (sometimes) faster version of ClusterSeq.`):
print(`Given a consecutive pattern`):
print(`uses the Umbral scheme to find the first N terms`):
print(`of the cluster-series for counting permutations`):
print(`according to the number of occurrences of p.`):
print(`For example, try:`):
print(`ClusterSeqF([1,2,3],t,10);`):

elif nops([args])=1 and op(1,[args])=ClusterSeqDir then
print(`ClusterSeqDir(p,t,N): like ClusterSeq(p,t,N)`):
print(`but using the Direct Umbral approach`):
print(`Given a consecutive pattern`):
print(`uses the Umbral scheme to find the first N terms`):
print(`of the cluster-series for counting permutations`):
print(`according to the number of occurrences of p.`):
print(`For example, try:`):
print(`ClusterSeqDir([1,2,3],t,10);`):

elif nops([args])=1 and op(1,[args])=ClusterSeqS then
print(`ClusterSeqS(P,t,N): Given a set of consecutive patterns, P,`):
print(`uses the Umbral matrix-scheme to find the first N terms`):
print(`of the cluster-series for counting permutations`):
print(`according to the number of occurrences of patterns in P.`):
print(`WARNING: VERY SLOW, NOT OPTIMIZED!`):
print(`For example, try:`):
print(`ClusterSeqS({[1,2,3]},t,10);`):


elif nops([args])=1 and op(1,[args])=FastestFriend then
print(`FastestFriend(p,N): The fastest friend of the pattern p`):
print(`for producing the N first terms of the avoiding sequence`):

elif nops([args])=1 and op(1,[args])=FEx1 then
print(`FEx1(U,var): inputs an umbral operator in the set of variables var`):
print(`and outputs the set variables that retain their exponent`):
print(`when applied to the monomial `):
print(`var[1]*var[2]^2*...*var[nops(var)-1]^(nops(var)-1)*`):
print(`var[nops(var)^(nops(var)-1)`):
print(`For example, try:`):
print(`FEx1(PtoU([3,2,1],0,x,z),[x[1],x[2],x[3],z]);`):

elif nops([args])=1 and op(1,[args])=FillIn then
print(`FillIn(N): given a vecor with 0's and positive integers`):
print(`outputs all the increasing sequences such that the`):
print(`non-zero entries are the same. For example, try:`):
print(`FillIn([0,3,0,0,9]);`):

elif nops([args])=1 and op(1,[args])=Foll then
print(`Foll(sig,P): given a consecutive pattern sig and a set of patterns P`):
print(`such that sig belongs to P, finds all triples [sigma,a,pi]`):
print(`such that pi is in P and a is an integer (less than`):
print(`the lenghts of pi and sigma) such that`):
print(`the last a entries of sigma reduce to the first a entries`):
print(`of pi. For example try Foll([1,2,3,4],{[1,2,3,4]});`):

elif nops([args])=1 and op(1,[args])=GenSeq then
print(`GenSeq(p,t,x,N): Given a consecutive pattern`):
print(`uses the Umbral scheme to find the first N terms`):
print(`of the cluster-series for counting permutations`):
print(`according to the number of occurrences of p, and the tails.`):
print(`For example, try:`):
print(`GenSeq([1,2,3],t,x,10);`):


elif nops([args])=1 and op(1,[args])=GroupU then
print(`GroupU(U): Given an Umbral operator U, splits it into`):
print(`a union of simpler ones each with a common denominator.`):
print(`The format now is a sequence of simpler umbral operators`):
print(`each with the same denominator`):
print(`For example, try:`):
print(`GroupU(PtoU([1,3,2],t,x,z));`):


elif nops([args])=1 and op(1,[args])=gTemplate then
print(`gTemplate(T,i,n): Given a triplet T=[sig,a,pi] and a symbol`):
print(`i (for indexing i[1],..., i[nops(pi)]) to describe`):
print(`symbolically how the sigma-part can be made`):
print(`0 denotes blank, and also indicating the n`):
print(`For example, try:`):
print(`gTemplate([[1,2,3],1,[1,2,3]],i,n);`):

elif nops([args])=1 and op(1,[args])=GuessPol then
print(`GuessPol(L,n,s0): guesses a polynomial of degree d in n for`):
print(`the list L, such that L[i]=P[i] for i>=s0 for example, try: `):
print(`GuessPol([(seq(i,i=1..10),n,1);`):

elif nops([args])=1 and op(1,[args])=Hakol then
print(`Hakol(x,r,n): The sum of x[1]^a1*...*x[r]^ar over`):
print(`all increasing sequences 1<=a1<a2<...<ar<=n.`):
print(`For example, try:`):
print(`Hakol(x,3,5);`):


elif nops([args])=1 and op(1,[args])=Hochakha then
print(`Hochakha(p,t): inputs a pattern, and outputs the theorem`):
print(`and its humanly-readable PROOF, `):
print(`that enables fast (poly-time, hopefully optimal)`):
print(`computation of the enumerating sequence for the`):
print(`number of n-permutations avoiding the consecutive Wilf`):
print(`pattern p. For example, try:`):
print(`Hochakha([1,2,3],t);`):

elif nops([args])=1 and op(1,[args])=khaverim then
print(`khaverim(p): all the friends of the pattern p.`):
print(`For example, try: khaverim([1,2,3,4]);`):

elif nops([args])=1 and op(1,[args])=Mishpat then
print(`Mishpat(p,t): inputs a pattern, and a variable t`):
print(`or t=0 `):
print(`and outputs the theorem`):
print(`that enables fast (poly-time, hopefully optimal)`):
print(`computation of the weight-enumerating sequence for`):
print(`n-permutations `):
print(`in z and t, such the coeff. of`, z^n*t^i):
print(`is the number of n-permutations with exactly i occurrences of`):
print(`of the consecutive Wilf pattern p`):
print(`Putting t=0, gives the exponential generating function, in z`):
print(`for the sequence enumerating the number of n-permutations`):
print(`avoiding p`):
print(`For example, try:`):
print(`Mishpat([1,2,3],t);`):
print(`Mishpat([1,2,3],0);`):

elif nops([args])=1 and op(1,[args])=Moments then
print(`Moments(p,R,n): The  first R moments about the mean`):
print(`of the random variable #number of occurrences of pattern`):
print(`p in uniformly-at-random n-permutation.`):
print(`For example try: `):
print(`Moments([1,2,3],4,n);`):

elif nops([args])=1 and op(1,[args])=Pred then
print(`Pred(pi,P): given a pattern pi and a set of patterns P`):
print(`such that pi belongs to P, finds all triples [sigma,a,pi]`):
print(`such that sigma is in P and a is an integer less than `):
print(`the lenghts of pi and sigam, such that`):
print(`the last a entries of sigma reduce to the first a entries`):
print(`of pi. For example try`):
print(` Pred([1,2,3,4],{[1,2,3,4]});`):


elif nops([args])=1 and op(1,[args])=PtoBU then
print(`PtoBU(p,t,x,z): The best Umbral operators, `):
print(`suitably reduced, amongts all images of p`):
print(`under trivial consecutive Wilf-equivalence.`):
print(`The output is a set of lists {[p1,U,mono,var]}`):
print(`where p1 is a pattern, U is the reduced umbral operator`):
print(`mono is the initial monomial, and var is the set`):
print(`of variables.`):
print(`For example, try:`):
print(`PtoBU([1,3,2],t,x,z);`):

elif nops([args])=1 and op(1,[args])=PtoU then
print(`PtoU(p,t,x,z): The Umbral operator, phrased in terms of`):
print(`x[1],x[2], x[nops(p)],z for the cluster generating function`):
print(`the single pattern p with weight (t-1)^#of components`):
print(`For example, try:`):
print(`PtoU([1,3,2],t,x,z);`):


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. For example, try:`):
print(`redu([5,4,2,11]);`):

elif nops([args])=1 and op(1,[args])=ReduceU  then
print(`ReduceU(U ,var,kv,mono): Given an umbral operator {[coeff,subs]} `):
print(`and a list of variables var, and a subset, kv, of [1,nops(var)]`):
print(`and an initial monomial`):
print(`finds the correponding umbral operator, and the correponding`):
print(`set of variables when one does the change of variables`):
print(`f(var)=mul(var[i1]^i1,i1 in kv)*g(var minus var[kv])`):
print(`and the modified initial monomial.`):
print(`the output is a list of length 3.`):
print(`For example, do`):
print(`ReduceU(PtoU([2,3,1],0,x,z),[x[1],x[2],x[3],z],{1,2},-x[1]*x[2]^2*x[3]^3*z^3);`):

elif nops([args])=1 and op(1,[args])=ReduceU1 then
print(`ReduceU1(U1,var,kv): Given an umbral atom U1 [coeff,subs],`):
print(`and a list of variables var, and a subset, kv, of [1,nops(var)]`):
print(`finds the correponding umbral atom, and the correponding`):
print(`set of variables when one does the change of variables`):
print(`f(var)=mul(var[i1]^i1,i1 in kv)*g(var minus var[kv])`):
print(`For example, do`):
print(`ReduceU1(PtoU([2,3,1],0,x,z)[1],[x[1],x[2],x[3],z],{1,2});`):

elif nops([args])=1 and op(1,[args])=RedVar then
print(`RedVar(U,ylist,var1): Given an umbral operator in`):
print(`the list of variables ylist, and an integer i between`):
print(`1 and nops(ylist),`):
print(`decides whether it is redundant when`):
print(`you plug-in ylist[i]=1. If it is, it returns the simplified`):
print(`umbral operator in the corresponding reduced set of variables.`):
print(`For example try`):
print(`RedVar(PtoU([1,3,2],0,x,z),[x[1],x[2],x[3],z],1);`):



elif nops([args])=1 and op(1,[args])=reps then
print(`reps(k): all the representatives of single patterns of length k`):
print(`For example, try:`):
print(`reps(4);`):


elif nops([args])=1 and op(1,[args])=Scheme then
print(`Scheme(P,i,t): Given a set of patterns P, and a symbol i,`):
print(`and a variable t (or 0)`):
print(`outputs a summation scheme for computing the weight-enumetators`):
print(`of the clusters according to the weight t^(#occurrences of patterns)`):
print(`In particular taking t=0 would give a scheme for counting`):
print(`the number of permutations avoiding the pattern in P.`):
print(`The output consists of a list of length 2.`):
print(`The first element is a list of the number of arguments`):
print(`in the functions (labelled naturally 1,2,3..)`):
print(`The second is a list of sets where the i-th item`):
print(`tells you the expression, in codified form for `):
print(`expressing the i-th function evaluated at`):
print(`(n;i[1],i[2], ...). Each member of  that set`):
print(`is a list consisting (i) the coeff. `):
print(`(ii) the function's numerical name (iii) the negative shift in n`):
print(`(iv)The set of {(n-shift; j1,j2, ...)} depending in i[1],i[2], ...`):
print(`over which the function (ii) should be summed. For exaple, try:`):
print(`Scheme({[1,2,3],[1,3,2]},i,t);`):


elif nops([args])=1 and op(1,[args])=schumT1 then
print(`schumT1(a,b,V): the sum of V[1]^i[1]*...V[r]^i[r] over all`):
print(`increasing integer vectors of length r such that`):
print(`a<i[1]<...<i[r]<b`):
print(`For example, try:`):
print(`schumT1(a,b,[x1]);`):

elif nops([args])=1 and op(1,[args])=schumT then
print(`schumT(L,V): given a list L consisting of 0's and symbols`):
print(`and a list of variables V of the same length (say r), finds`):
print(`the symbolic sum of x[1]^i[1]*...*x[r]^i[r]`):
print(`such that i[j]=L[j] if L[j] is not 0 and otherwise`):
print(`ranging over all increasing sequences compatible with`):
print(` L. For example, try: `):
print(` schumT([a,0,b],[x,y,z]); `):

elif nops([args])=1 and op(1,[args])=SeferP then
print(`SeferP(k,N,accur1,n,N1): `):
print(`Inputs a pos. integer k, a pos. integer N, a pos. integer`):
print(`accur1, a symbol n, and a pos. integer N1, outputs`):
print(`a web-book with the statment of theorems, followed by proof,`):
print(`that enable the optimal polynomial-time enumeration`):
print(`of the terms of the enumerating sequence for n-permutations`):
print(`avoiding single consecutive patterns of length k.`):
print(`N1 is a pos. integer less than N that tells how many`):
print(`terms to display in the book (too many terms could be`):
print(`tiring)`):
print(`For example, try:`):
print(`SeferP(3,50,10,n,15);`):

elif nops([args])=1 and op(1,[args])=SeferT then
print(`SeferT(k,N,accur1,n,N1): `):
print(`Inputs a pos. integer k, a pos. integer N, a pos. integer`):
print(`accur1, a symbol n, and a pos. integer N1, outputs`):
print(`a web-book with the statment of theorems (no proofs)a`):
print(`that enable the optimal polynomial-time enumeration`):
print(`of the terms of the enumerating sequence for n-permutations`):
print(`avoiding single consecutive patterns of length k.`):
print(`N1 is a pos. integer less than N that tells how many`):
print(`terms to display in the book (too many terms could be`):
print(`tiring)`):
print(`For example, try:`):
print(`SeferT(3,50,10,n,15);`):

elif nops([args])=1 and op(1,[args])=Seq then
print(`Seq(p,N): `):
print(`Given a pattern p`):
print(` a pos. integer N, finds the first N terms of the`):
print(`sequence enumerating permutations avoiding P`):
print(`For example, try:`):
print(`Seq([1,2,3],10);`):

elif nops([args])=1 and op(1,[args])=SeqF then
print(`SeqF(p,N):a (sometimes) faster version of Seq(p,N) `):
print(`Given a pattern p`):
print(` a pos. integer N, finds the first N terms of the`):
print(`sequence enumerating permutations avoiding P`):
print(`For example, try:`):
print(`SeqF([1,2,3],10);`):

elif nops([args])=1 and op(1,[args])=SeqFromDE then
print(`SeqFromDE(DE1,D1,z,N):  Given a differential equation`):
print(`[Oper,Mono] in terms of the diff. operator D1, and the`):
print(`variable z, finds the  first N terms of the sequence`):
print(`of coeff.s. For example, try:`):
print(`SeqFromDE(UtoD(op(PtoBU([2,3,1],t,x,z)[1][2]),D1),D1,z,30);`):

elif nops([args])=1 and op(1,[args])=SeqS then
print(`SeqS(P,N): `):
print(`Given a set of consecutive patterns P`):
print(` a pos. integer N, finds the first N terms of the`):
print(`sequence enumerating permutations avoiding the members of P`):
print(`WARNING: VERY SLOW, NOT OPTIMIZED!`):
print(`For example, try:`):
print(`SeqS({[1,2,3]},10);`):

elif nops([args])=1 and op(1,[args])=SeqT then
print(`SeqT(P,t,N): Given a set of patterns P, a variable t, and`):
print(` a pos. integer N, finds the first N terms of the`):
print(`sequence enumerating permutations according to the`):
print(`number of occurrences of the consecutive patterns P. With t=0`):
print(`you get the number of permutations avoiding the consecutive`):
print(`patterns of P. For example, try:`):
print(`SeqT({[1,2,3]},t,10);`):

elif nops([args])=1 and op(1,[args])=SeqTu then
print(`SeqTu(p,t,N): Like SeqT({p},t,N) but only for`):
print(`single pattern, and using the umbral method.`):
print(`Given a set of patterns P, a variable t, and`):
print(` a pos. integer N, finds the first N terms of the`):
print(`sequence enumerating permutations according to the`):
print(`number of occurrences of the consecutive patterns P. With t=0`):
print(`you get the number of permutations avoiding the consecutive`):
print(`patterns of P. For example, try:`):
print(`SeqTu([1,2,3],t,10);`):


elif nops([args])=1 and op(1,[args])=Sipur then
print(`Sipur(k,N,accur1,n): All the sequences enumerating`):
print(`permutations avoiding a single pattern of length k`):
print(`together with their guessed asymptotics to accuracy`):
print(`of accur1  digits. n is the symbol for the length of the perm.`):
print(`For example, try:`):
print(`Sipur(3,50,10,n);`):


elif nops([args])=1 and op(1,[args])=SipurF then
print(`SipurF(k,N,accur1,n): a (sometimes) faster version of Sipur.`):
print(`All the sequences enumerating`):
print(`permutations avoiding a single pattern of length k`):
print(`together with their guessed asymptotics to accuracy`):
print(`of accur1  digits. n is the symbol for the length of the perm.`):
print(`For example, try:`):
print(`SipurF(3,50,10,n);`):

elif nops([args])=1 and op(1,[args])=SpellOutUmbra then
print(`SpellOutUmbra(U,F): Given an umbral operator`):
print(`U in the variable list ylist and a symbol F`):
print(` for the function name, spells out the operator`):
print(`in human language. For example, try:`):
print(`SpellOutUmbra({[1,[q*x]]},F); `):


elif nops([args])=1 and op(1,[args])=SPtoUM then
print(`SPtoUM(P,t,x,z,A): The matrix of Umbral operator, phrased in terms of`):
print(`x[1],x[2], x[nops(p)],z for the cluster generating function`):
print(`the list of patterns P with weight (t-1)^#of components`):
print(`A is a symbol for indexing the next pattern`):
print(`For example, try:`):
print(`SPtoUM([[1,3,2]],t,x,z,A);`):

elif nops([args])=1 and op(1,[args])=Template then
print(`Template(T,i): Given a triplet T=[sig,a,pi] and a symbol`):
print(`i (for indexing i[1],..., i[nops(pi)]) to describe`):
print(`symbolically how the sigma-part can be made`):
print(`0 denotes blank`):
print(`For example, try:`):
print(`Template([[1,2,3],1,[1,2,3]],i);`):

elif nops([args])=1 and op(1,[args])=TemplateUP then
print(`TemplateUP(T,j): Given a triplet T=[sig,a,pi] and a symbol`):
print(`j (for indexing j[1],..., j[nops(sig)]) to describe`):
print(`symbolically how the pi-part can be made from the sig part`):
print(`0 denotes blank`):
print(`For example, try:`):
print(`TemplateUP([[1,2,3],1,[1,2,3]],j);`):

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 try:`):
print(`ToUmbra(a*x^b+b*y^a,[x,y],[a,b]);`):



elif nops([args])=1 and op(1,[args])=UpScheme then
print(`UpScheme(p): Given a SINGLE consecutive pattern p,`):
print(`and a symbol j to denote an indexed variable`):
print(`and a symbol n, produces an upward scheme`):
print(`listing the set of upwards templates`):
print(`that a monomial x[1]^j[1]*x[2]^j[2]*...*z^n can go to`):
print(`For example, try:`):
print(`UpScheme([1,2,3],j,n);`):


elif nops([args])=1 and op(1,[args])=UtoD then
print(`UtoD(U,var,mono,D1): Given a functional operator U, in the set of`):
print(`variables var, and initial monomial mono, tries to`):
print(`find an ordinary differential operator Op(z,D1)`):
print(`such that the solution to the functional equation`):
print(`f(var)=mono+Uf(var) upon plugging the catalytic variables in var`):
print(`to be all 1's except the last, that we call z,`):
print(`you get the differential equation for F(z)=f(1,1,1, ...,z)`):
print(`to be F(z)=Mono+Oper(z,D1), where Mono is mono with all`):
print(`the var's (except z) equal 1. `):
print(`It returns the length-2 list [Oper(z,D1),Mono] if succesful,`):
print(`or FAIL, if not.`):
print(`For example, try:`):
print(`UtoD(op(PtoBU([3,2,1],t,x,z)[1][2]),D1);`):

elif nops([args])=1 and op(1,[args])=Val then
print(`Val(S,i,n0,t): The sum of all values of all functions`):
print(`that show up in the scheme S for n=n0.`):
print(`For example, try:`):
print(`Val(Scheme({[1,2,3]},i,t),i,6,t);`):

elif nops([args])=1 and op(1,[args])=Val1 then
print(`Val1(a,S,i,n0,i0,t): The value of function #a in the scheme S`):
print(`at the numerical vector i0 when n=n0. For example, try:`):
print(`Val1(1,Scheme({[1,2,3]},i,t),i,6,[2,3,5],t);`):


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

ez:

ez1:=proc(): 
print(` Actualize(Temp,i,i0,n0) `):
end:


with(combinat):

#redu(perm): Given a permutation on a set of positive integers
#finds its reduced form. For example, try:
#redu([5,4,2,11]);
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:

#Pred(pi,P): given a pattern pi and a set of patterns P
#such that pi belongs to P, finds all triples [sigma,a,pi]
#such that sigma is in P and a is an integer (less than
#the lenghts of pi and sigma) such that
#the last a entries of sigma reduce to the first a entries
#of pi. For example try Pred([1,2,3,4],{[1,2,3,4]});
Pred:=proc(pi,P) local  a,gu,sig:
gu:={}:

for sig in P do
 for a from 1 to min(nops(pi),nops(sig))-1 do
  if redu([op(nops(sig)-a+1..nops(sig),sig)])=redu([op(1..a,pi)]) then
   gu:=gu union {[sig,a,pi]}:
  fi:
 od:
od:

gu:

end:

#KatanMimcha(i,S): given an integer i and a set of integers S
#finds how many members of S are less than i
#For example, try:
#KatanMimcha(4,{2,3,8,9});
KatanMimcha:=proc(i,S) local co,s:
co:=0:

for s in S do
 if s<i then
   co:=co+1:
 fi:
od:
co:
end:

#Template(T,i): Given a triplet T=[sig,a,pi] and a symbol
#i (for indexing i[1],..., i[nops(pi)]) to describe
#symbolically how the sigma-part can be made
#0 denotes blank
#For example, try:
#Template([[1,2,3],1,[1,2,3]],i);
Template:=proc(T,i)
local sig,a,pi,S,lu,i1:

sig:=T[1]: a:=T[2]: pi:=T[3]:

if a>=min(nops(sig),nops(pi)) then
 RETURN(FAIL):
fi:

if redu([op(nops(sig)-a+1..nops(sig),sig)])<>redu([op(1..a,pi)]) then
 RETURN(FAIL):
fi:

lu:={op(a+1..nops(pi),pi)}:

for i1 from 1 to nops(sig)-a do
 S[sig[i1]]:=0:
od:

for i1 from 1 to a do
 S[sig[nops(sig)-a+i1]]:=i[pi[i1]]-KatanMimcha(pi[i1],lu):
od:

[seq(S[i1],i1=1..nops(sig))]:

end:



#Scheme(P,i,t): Given a set of patterns P, and a symbol i,
#and a variable t (or 0)
#outputs a summation scheme for computing the weight-enumetators
#of the clusters according to the weight t^(#occurrences of patterns)
#In particular taking t=0 would give a scheme for counting
#the number of permutations avoiding the pattern in P.
#The output consists of a list of length 2.
#The first element is a list of the number of arguments
#in the functions (labelled naturally 1,2,3..)
#The second is a list of sets where the i-th item
#tells you the expression, in codified form for 
#expressing the i-th function evaluated at
#(n;i[1],i[2], ...). Each member of the that set
#is a list consisting (i) the coeff. 
#(ii) the function's numerical name (iii) the negative shift in n
#(iv)The set of {(n-shift; j1,j2, ...)} depending in i[1],i[2], ...
#over which the function (ii) should be summed. For exaple, try:
#Scheme({[1,2,3],[1,3,2]},i,t);
Scheme:=proc(P,i,t) local gu,pi,mu,T1, T1i,i1,ku,mu1:

for i1 from 1 to nops(P) do
 T1[P[i1]]:=i1:
 T1i[i1]:=P[i1]:
od:

ku:=[seq(nops(T1i[i1]),i1=1..nops(P))]:

gu:=[]:

for i1 from 1 to nops(P) do
 pi:=T1i[i1]:
 mu:=Pred(pi,P):
gu:=[op(gu),
{seq([t-1,T1[mu1[1]],nops(pi)-mu1[2], Template(mu1,i)],mu1 in mu)}]:
od:

[ku,gu]:
end:

#Actualize(Temp,i,i0,n0): given a template in terms of symbols
#i[1],i[2],i[3],... and  a specfic numerical increasing vector
#of integers i0, and a pos. integer n0, finds all the
#possible vectors of integers with those specifications
#For example, try:
#Actualize([0,0,i[2]-1],i,[5,7,8],11);
Actualize:=proc(Temp,i,i0,n0) local i1,mu,mu1,Temp1,X:
Temp1:=subs(0=X,Temp):
mu:=subs({seq(i[i1]=i0[i1],i1=1..nops(i0))},Temp1):
if member(0,mu) then
 RETURN({}):
fi:
mu:=subs(X=0,mu):
mu:=[op(mu),n0+1]:
mu:=FillIn(mu):
{seq([op(1..nops(mu1)-1,mu1)],mu1 in mu)}:
end:

#FillIn(N): given a vecor with 0's and positive integers
#outputs all the increasing sequences such that the
#non-zero entries are the same. For example, try:
#FillIn([0,3,0,0,9]);
FillIn:=proc(N) local i1,j1,a:
for i1 from 1 to nops(N) while N[i1]<>0 do od:


if i1=nops(N)+1 then
 RETURN({N}):
fi:

for j1 from i1+1 to nops(N) while N[j1]=0 do od:
j1:=j1-1:


if i1=1 then
{seq(op(FillIn([op(1..j1-1,N),a,op(j1+1..nops(N),N)])),a=j1..N[j1+1]-1)}:

else

{seq(op(FillIn([op(1..j1-1,N),a,op(j1+1..nops(N),N)])),a=N[i1-1]+j1-i1+1..N[j1+1]-1)}:
fi:



end:


#Val1(a,S,i,n0,i0,t): The value of function #a in the scheme S
#at the numerical vector i0 when n=n0. For example, try:
#Val1(1,Scheme({[1,2,3]},i,t),i,6,[2,3,5],t);
Val1:=proc(a,S,i,n0,i0,t) local co,mu,mu1,coe,fn,kv,shif,b,i0a:
option remember:

if not (a>=1 and a<=nops(S[1])) then
 RETURN(FAIL):
fi:

if S[1][a]<>nops(i0) then 
 RETURN(FAIL):
fi:

if S[1][a]=n0 then
 if i0=[seq(b,b=1..n0)] then
    RETURN(t-1):
  else
    RETURN(0):
 fi:
fi:

if n0<S[1][a] then
  RETURN(0):
fi:

co:=0:
mu:=S[2][a]:

for mu1 in mu do
 coe:=mu1[1]:
 fn:=mu1[2]:
 shif:=mu1[3]:
 kv:=mu1[4]:

kv:=Actualize(kv,i,i0,n0-shif):

co:=expand(co+coe*add(Val1(fn,S,i,n0-shif,i0a,t),i0a in kv)):

od:

co:

end:

#Val(S,i,n0,t): The sum of all values of all functions
#that show up in the scheme S for n=n0.
#For example, try:
#Val(Scheme({[1,2,3]},i,t),i,6,t);
Val:=proc(S,i,n0,t) local co,i0,kv,a,b:
co:=0:

for a from 1 to nops(S[1]) do
 kv:=choose([seq(b,b=1..n0)],S[1][a]):
 co:=co+add(Val1(a,S,i,n0,i0,t), i0 in kv):
od:

co:

end:


#SeqT(P,t,N): Given a set of patterns P, a variable t, and
# a pos. integer N, finds the first N terms of the
#sequence enumerating permutations according to the
#number of occurrences of the consecutive patterns P. With t=0
#you get the number of permutations avoiding the consecutive
#patterns of P. For example, try:
#SeqT({[1,2,3]},t,10);
SeqT:=proc(P,t,N) local gu,n0,i,i1,z,f:
gu:=Scheme(P,i,t):
gu:=[seq(Val(gu,i,n0,t),n0=1..N)]:
f:=1-z-add(gu[i]*z^i/i!,i=2..N):

f:=taylor(1/f,z=0,N+1):
[seq(expand(coeff(f,z,i1))*i1!,i1=1..N)]:
end:




#Hakol(x,r,n): The sum of x[1]^a1*...*x[r]^ar over
#all increasing sequences 1<=a1<a2<...<ar<=n.
#For example, try:
#Hakol(x,3,5);
Hakol:=proc(x,r,n) local i,gu,z,j:
gu:=mul(x[i]^i,i=1..r)*z^r/
mul(1-mul(x[j],j=i..r)*z,i=1..r)/(1-z):
gu:=taylor(gu,z=0,n+1):
expand(coeff(gu,z,n)):
end:


#gTemplate(T,i,n): Given a triplet T=[sig,a,pi] and a symbol
#i (for indexing i[1],..., i[nops(pi)]) and as symbol n, to describe
#symbolically how the sigma-part can be made
#0 denotes blank
#For example, try:
#gTemplate([[1,2,3],1,[1,2,3]],i,n);
gTemplate:=proc(T,i,n) local gu:
gu:=Template(T,i):
[op(gu),n-(nops(T[3])-T[2])]:
end:



#schumT1(a,b,V): the sum of V[1]^i[1]*...V[r]^i[r] over all
#increasing integer vectors of length r such that
#a<i[1]<...<i[r]<b
#For example, try:
#schumT1(a,b,[x1]);
schumT1:=proc(a,b,V) local r,i1,i2,i3,i4,i5,i6:
r:=nops(V):
if r=1 then
  RETURN(expand(sum(V[1]^i1,i1=a+1..b-1))):

elif r=2 then
sum(sum(V[1]^i1*V[2]^i2,i1=a+1..i2-1),i2=a+2..b-1):
elif r=3 then
sum(sum(sum(V[1]^i1*V[2]^i2*V[3]^i3,i1=a+1..i2-1),i2=a+2..i3-1),i3=a+3..b-1):
elif r=4 then
sum(sum(sum(sum(V[1]^i1*V[2]^i2*V[3]^i3*V[4]^i4
,i1=a+1..i2-1),i2=a+2..i3-1),i3=a+3..i4-1),i4=a+4..b-1):

elif r=5 then
sum(sum(sum(sum(sum(V[1]^i1*V[2]^i2*V[3]^i3*V[4]^i4*V[5]^i5
,i1=a+1..i2-1),i2=a+2..i3-1),i3=a+3..i4-1),i4=a+4..i5-1),i5=a+5..b-1):

elif r=6 then
sum(
sum(sum(sum(sum(sum(V[1]^i1*V[2]^i2*V[3]^i3*V[4]^i4*V[5]^i5*V[6]^i6
,i1=a+1..i2-1),i2=a+2..i3-1),i3=a+3..i4-1),i4=a+4..i5-1),i5=a+5..i6-1),
i6=a+6..b-1):

else
 print(`Not yet implemented`):
 RETURN(FAIL):
fi:

end:



#TemplateUP(T,j): Given a triplet T=[sig,a,pi] and a symbol
#i (for indexing i[1],..., i[nops(pi)]) to describe
#symbolically how the pi-part can be made from the sig part
#0 denotes blank
#For example, try:
#TemplateUP([[1,2,3],1,[1,2,3]],j);
TemplateUP:=proc(T,j)
local sig,a,pi,S,lu,i1:

sig:=T[1]: a:=T[2]: pi:=T[3]:

if a>=min(nops(sig),nops(pi)) then
 RETURN(FAIL):
fi:

if redu([op(nops(sig)-a+1..nops(sig),sig)])<>redu([op(1..a,pi)]) then
 RETURN(FAIL):
fi:

lu:={op(a+1..nops(pi),pi)}:

for i1 from a+1 to nops(pi) do
 S[pi[i1]]:=0:
od:

for i1 from 1 to a do
 S[pi[i1]]:=j[sig[nops(sig)-a+i1]]+KatanMimcha(pi[i1],lu):
od:

[seq(S[i1],i1=1..nops(pi))]:

end:


#schumTst(L,V,st): given a list L consisting of 0's and symbols
#and a list of variables V of the same length (say r), finds
#the symbolic sum of x[1]^i[1]*...*x[r]^i[r]
#such that i[j]=L[j] if L[j] is not 0 and otherwise
#ranging over all increasing sequences compatible with
#L. For example, try:
#schumTst([a,0,b],[x,y,z],0);
schumTst:=proc(L,V,st) local gu,i,L1,V1:

if nops(L)<>nops(V) then
 print(`The first two arguments must be lists of the same size`):
 RETURN(FAIL):
fi:

for i from 1 to nops(L) while L[i]=0 do od:
if i=nops(L)+1 then
 RETURN(FAIL):
fi:



if i>1 then
gu:=schumT1(st,L[i],[op(1..i-1,V)]):
else
gu:=1:
fi:
gu:=gu*V[i]^L[i]:

if i=nops(L) then
 RETURN(gu):
else
L1:=[op(i+1..nops(L),L)]:
V1:=[op(i+1..nops(V),V)]:

RETURN(gu*schumTst(L1,V1,L[i])):
fi:

end:


schumT:=proc(L,V) local gu,i,L1,V1:
schumTst(L,V,0):
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 try:
#ToUmbra(a*x^b+b*y^a,[x,y],[a,b]);

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,SubsList,gu,bu:
 
FRONT:=Umbra1[1]:
SubsList:=Umbra1[2]:
 

 
 
if nops(SubsList)<>nops(ylist) then
 ERROR(`SubsList and ylist should have the same length`):
fi:
 

 
gu:=subs({seq(ylist[i]=SubsList[i],i=1..nops(ylist))},f): 

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################






#Foll(sig,P): given a consecutive pattern sig and a set of patterns P
#such that sig belongs to P, finds all triples [sigma,a,pi]
#such that pi is in P and a is an integer (less than
#the lenghts of pi and sigma) such that
#the last a entries of sigma reduce to the first a entries
#of pi. For example try Foll([1,2,3,4],{[1,2,3,4]});
Foll:=proc(sig,P) local  a,gu,pi:
gu:={}:

for pi in P do
 for a from 1 to min(nops(pi),nops(sig))-1 do
  if redu([op(nops(sig)-a+1..nops(sig),sig)])=redu([op(1..a,pi)]) then
   gu:=gu union {[sig,a,pi]}:
  fi:
 od:
od:

gu:

end:


#CleanUp(U): cleans up the Umbra U from 0
CleanUp:=proc(U) local U1,u1:
U1:={}:
for u1 in U do
if u1[1]<>0 then
 U1:=U1 union {[normal(u1[1]),op(2..nops(u1),u1)]}:
fi:
od:
U1:
end:


#PtoU(p,t,x,z): The Umbral operator, phrased in terms of
#x[1],x[2], x[nops(p)],z for the cluster generating function
#the single pattern p with weight (t-1)^#of components
#For example, try:
#PtoU([1,3,2],t,x,z);
PtoU:=proc(p,t,x,z) local gu,lu,j,lu1,ku1,su1,n,i1,xlist,alist:
option remember:
lu:=Foll(p,{p}):
gu:={}:

for lu1 in lu do
 ku1:=TemplateUP(lu1,j):
 ku1:=[op(ku1),n+nops(p)-lu1[2]+1]:
 xlist:=[seq(x[i1],i1=1..nops(p)),z]:
 su1:=(t-1)*schumT(ku1,xlist)/z:
 su1:=expand(su1):
 alist:=[seq(j[i1],i1=1..nops(p)),n]:
 gu:=gu union  ToUmbra(su1,xlist,alist):
od:

if {seq(gu[i1][2],i1=1..nops(gu))}<>{[0$(nops(p)+1)]} then
print(`Something is wrong`):
print({seq(gu[i1][2],i1=1..nops(gu))}):
RETURN(FAIL):
fi:
{seq([gu[i1][1],gu[i1][3]],i1=1..nops(gu))}:
end:



#GenSeqSlow(p,t,x,N): Given a consecutive pattern
#uses the Umbral scheme to find the first N terms
#of the cluster-series for counting permutations
#according to the number of occurrences of p.
#For example, try:
#GenSeqSlow([1,2,3],t,x,10);

GenSeqSlow:=proc(p,t,x,N) local U, Ini1,z,gu,mu,i1,ylist,i2:
U:=PtoU(p,t,x,z):
ylist:=[seq(x[i1],i1=1..nops(p)),z]:
Ini1:=z^nops(p)*(t-1)*mul(x[i1]^i1,i1=1..nops(p)):
mu:=Ini1:
gu:=Ini1:

for i1 from 1 to N do
 mu:=expand(ApplyUmbra(U,mu,ylist)):
 mu:=add(coeff(mu,z,i2)*z^i2,i2=0..N):
 gu:=gu+mu:
od:
gu:=expand(gu):
[seq(coeff(gu,z,i1),i1=1..N)]:
end:



#SeqTu(P,t,N): Like SeqT(P,t,N) but only for
#single pattern, and using the umbral method.
#Given a set of patterns P, a variable t, and
# a pos. integer N, finds the first N terms of the
#sequence enumerating permutations according to the
#number of occurrences of the consecutive patterns P. With t=0
#you get the number of permutations avoiding the consecutive
#patterns of P. For example, try:
#SeqTu({[1,2,3]},t,10);
SeqTuSlow:=proc(P,t,N) local gu,i1,x,z,f:
if nops(P)<>1 then
 print(`Multiple patterns are not yet implemented`):
 RETURN(FAIL):
fi:
gu:=GenSeq(P[1],t,x,N):
gu:=subs({seq(x[i1]=1,i1=1..nops(P[1]))},gu):
f:=1-z-add(gu[i1]*z^i1/i1!,i1=2..N):

f:=taylor(1/f,z=0,N+1):
[seq(expand(coeff(f,z,i1))*i1!,i1=1..N)]:
end:


#SeqSlow(P,N): 
#Given a set of patterns P (so far with one element), a variable t, and
# a pos. integer N, finds the first N terms of the
#sequence enumerating permutations avoiding P
#For example, try:
#SeqSlow({[1,2,3]},10);
SeqSlow:=proc(P,N):SeqTuSlow(P,0,N):end:



#GroupU(U): Given an Umbral operator U, splits it into
#a union of simpler ones each with a common denominator.
#The format now is a 
GroupU:=proc(U) local mekhanim,T,i,u1,m1:
mekhanim:={seq(denom(U[i][1]),i=1..nops(U))}:

for m1 in mekhanim do
T[m1]:={}:
od:

for u1 in U do
 T[denom(u1[1])]:=T[denom(u1[1])] union {u1}:
od:

[seq(T[m1],m1 in mekhanim)]:

end:



#GenSeq(p,t,x,N): Given a consecutive pattern
#uses the Umbral scheme to find the first N terms
#of the cluster-series for counting permutations
#according to the number of occurrences of p
#and the tails
#For example, try:
#GenSeq([1,2,3],t,x,10);

GenSeq:=proc(p,t,x,N) local U, Ini1,z,gu,mu,i1,ylist,i2,U1:
U:=PtoU(p,t,x,z):
U1:=GroupU(U):
ylist:=[seq(x[i1],i1=1..nops(p)),z]:
Ini1:=z^nops(p)*(t-1)*mul(x[i1]^i1,i1=1..nops(p)):
mu:=Ini1:
gu:=Ini1:

for i1 from 1 to N do
 mu:=expand(
normal(add(normal(ApplyUmbra(U1[i2],mu,ylist)),i2=1..nops(U1)))):
 mu:=add(coeff(mu,z,i2)*z^i2,i2=0..N):
 gu:=gu+mu:
od:
gu:=expand(gu):
[seq(coeff(gu,z,i1),i1=1..N)]:
end:





#ClusterSeqOld(p,t,x,N): Given a consecutive pattern
#uses the Umbral scheme to find the first N terms
#of the cluster-series for counting permutations
#according to the number of occurrences of p.
#For example, try:
#ClusterSeqOld([1,2,3],t,10);
ClusterSeqOld:=proc(p,t,N) local U, Ini1,z,gu,mu,i1,ylist,i2,U1,x:
U:=PtoU(p,t,x,z):
U1:=GroupU(U):
ylist:=[seq(x[i1],i1=1..nops(p)),z]:
Ini1:=z^nops(p)*(t-1)*mul(x[i1]^i1,i1=1..nops(p)):
mu:=Ini1:
gu:=z^nops(p)*(t-1):

for i1 from 1 to N do
 mu:=normal(
add(ApplyUmbra(U1[i2],mu,ylist),i2=1..nops(U1))):
 mu:=add(coeff(mu,z,i2)*z^i2,i2=0..N):
 gu:=gu+subs({seq(x[i2]=1,i2=1..nops(p))},mu):
od:
gu:=expand(gu):
[seq(coeff(gu,z,i1),i1=1..N)]:
end:



#SeqTu(p,t,N): Like SeqT({p},t,N) but only for
#single pattern, and using the umbral method.
#Given a set of patterns P, a variable t, and
# a pos. integer N, finds the first N terms of the
#sequence enumerating permutations according to the
#number of occurrences of the consecutive patterns P. With t=0
#you get the number of permutations avoiding the consecutive
#patterns of P. For example, try:
#SeqTu([1,2,3],t,10);
SeqTu:=proc(p,t,N) local gu,i1,x,z,f:

gu:=ClusterSeq(p,t,N):
f:=1-z-add(gu[i1]*z^i1/i1!,i1=2..N):

f:=taylor(1/f,z=0,N+1):
[seq(expand(coeff(f,z,i1))*i1!,i1=1..N)]:
end:


#SeqTuF(p,t,N): Like SeqT({p},t,N) but only for
#single pattern, and using the umbral method.
#Given a set of patterns P, a variable t, and
# a pos. integer N, finds the first N terms of the
#sequence enumerating permutations according to the
#number of occurrences of the consecutive patterns P. With t=0
#you get the number of permutations avoiding the consecutive
#patterns of P. For example, try:
#SeqTuF([1,2,3],t,10);
SeqTuF:=proc(p,t,N) local gu,i1,x,z,f:

gu:=ClusterSeqF(p,t,N)[1]:
f:=1-z-add(gu[i1]*z^i1/i1!,i1=2..N):

f:=taylor(1/f,z=0,N+1):
[seq(expand(coeff(f,z,i1))*i1!,i1=1..N)]:
end:

#Seq(p,N): 
#Given a  pattern p  and
# a pos. integer N, finds the first N terms of the
#sequence enumerating permutations avoiding the consecutive pattern p
#For example, try:
#Seq([1,2,3],10);
Seq:=proc(p,N) option remember:
SeqTu(p,0,N):end:

#SeqF(p,N): 
#Given a  pattern p  and
# a pos. integer N, finds the first N terms of the
#sequence enumerating permutations avoiding the consecutive pattern p
#For example, try:
#SeqF([1,2,3],10);
SeqF:=proc(p,N) option remember:
SeqTuF(p,0,N):end:


#SpellOutUmbra(U,F): Given an umbral operator
#U in the variable list ylist and a symbol F
#for the function name, spells out the operator
#in human language. For example, try:
#SpellOutUmbra({[1,[q*x]]},[x],F); 
SpellOutUmbra:=proc(U,F) local i:
add(U[i][1]*F(op(U[i][2])),i=1..nops(U)):
end:


#RedVar(U,ylist,var1): Given an umbral operator in
#the list of variables ylist, and an integer i between
#1 and nops(ylist),
#decides whether it is redundant when
#you plug-in ylist[i]=1. If it is, it returns the simplified
#umbral operator in the corresponding reduced set of variables.
#For example try:
#RedVar(PtoU([1,3,2],0,x,z),[x[1],x[2],x[3],z],1);
RedVar:=proc(U,ylist,i) local U1,i1,ku:

ku:={seq(expand(subs(ylist[i]=1,denom(U[i1][1]))),i1=1..nops(U))}:
if member(0,ku) then
 RETURN(FAIL):
fi:

U1:=subs(ylist[i]=1,U):

ku:={seq(U1[i1][2][i],i1=1..nops(U1))}:
if ku<>{1} then
   RETURN(FAIL):
fi:

U1:={seq([U1[i1][1],[op(1..i-1,U1[i1][2]),op(i+1..nops(U1[i1][2]),U1[i1][2])]],
     i1=1..nops(U1))}:


U1,[op(1..i-1,ylist),op(i+1..nops(ylist),ylist)]:
end:



#ClusterSeq(p,t,N): Given a consecutive pattern
#uses the Umbral scheme to find the first N terms
#of the cluster-series for counting permutations
#according to the number of occurrences of p.
#For example, try:
#ClusterSeq([1,2,3],t,10);
ClusterSeq:=proc(p,t,N) local U, Ini1,z,gu,mu,i1,ylist,i2,U1,x:
U:=PtoU(p,t,x,z):
ylist:=[seq(x[i1],i1=1..nops(p)),z]:

U1:=RedVar(U,ylist,p[1]):
if U1=FAIL then
 RETURN(FAIL):
fi:
ylist:=U1[2]:
U1:=U1[1]:




Ini1:=z^nops(p)*(t-1)*normal(mul(x[i1]^i1,i1=1..nops(p))/x[p[1]]^p[1]):
mu:=Ini1:
gu:=z^nops(p)*(t-1):

U1:=GroupU(U1):

for i1 from 1 to N do
 mu:=normal(
add(ApplyUmbra(U1[i2],mu,ylist),i2=1..nops(U1))):
 mu:=add(coeff(mu,z,i2)*z^i2,i2=0..N):
 gu:=gu+subs({seq(x[i2]=1,i2=1..nops(p))},mu):
od:

gu:=expand(gu):
[seq(coeff(gu,z,i1),i1=1..N)]:

end:


#SeqFromDE(DE1,D1,z,N):  Given a differential equation
#[Oper,Mono] in terms of the diff. operator D1, and the
#variable z, finds the  first N terms of the sequence
#of coeff.s. For example, try:
#SeqFromDE(UtoD(op(PtoBU(p,t,x,z)[1][2]),D1),D1,z,30):
SeqFromDE:=proc(DE1,D1,z,N) local mu,gu,ope,i:
ope:=DE1[1]: mu:=DE1[2]:
gu:=mu:

while ldegree(mu,z)<=N do
 mu:=coeff(ope,D1,0)*mu+add(coeff(ope,D1,i)*diff(mu,z$i),i=1..degree(ope,D1)):
 mu:=expand(mu):
 mu:=add(coeff(mu,z,i)*z^i,i=0..N):
 gu:=gu+mu:
od:

[seq(coeff(gu,z,i),i=1..N)]:

end:

#SPtoUM(P,t,x,z,A): The matrix of Umbral operator, phrased in terms of
#x[1],x[2], x[nops(p)],z for the cluster generating function
#the list of patterns P with weight (t-1)^#of components
#A is a symbol for indexing the next pattern
#For example, try:
#SPtoUM([[1,3,2]],t,x,z,A);
SPtoUM:=proc(P,t,x,z,A) local gu,lu,j,lu1,ku1,su1,n,i1,xlist,alist,p,Gu,i,T:

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

Gu:=[]:

for i from 1 to nops(P) do
p:=P[i]:
lu:=Foll(p,P):
gu:={}:

for lu1 in lu do
 ku1:=TemplateUP(lu1,j):
 ku1:=[op(ku1),n+nops(p)-lu1[2]+1]:
 xlist:=[seq(x[i1],i1=1..nops(p)),z]:
 su1:=(t-1)*schumT(ku1,xlist)/z*A[T[lu1[3]]]:
 su1:=expand(su1):
 alist:=[seq(j[i1],i1=1..nops(p)),n]:
 gu:=gu union  ToUmbra(su1,xlist,alist):
od:

if {seq(gu[i1][2],i1=1..nops(gu))}<>{[0$(nops(p)+1)]} then
print(`Something is wrong`):
print({seq(gu[i1][2],i1=1..nops(gu))}):
RETURN(FAIL):
fi:
Gu:=[op(Gu),{seq([gu[i1][1],gu[i1][3]],i1=1..nops(gu))}]:
od:
Gu:
end:



#ApplyUmbraMat(Umbra1Mat,f,ylist,A): Given a list of
#umbral operators Umbra1Mat phrased in terms of A[1],A[2]
#and a linear combination of A[1],A[2], ... with coefficients
#that are polynomials in ylist applies Umbra1Mat to f
#For example try:
#ApplyUmbraMat(SPtoUM([[1,3,2]],t,x,z,A),x[1]*x[2]^2*x[3]^3*z^3*A[1],[x[1],x[2],x[3],z]],A);
ApplyUmbraMat:=proc(Umbra1Mat,f,ylist,A) local i:

add(ApplyUmbra(Umbra1Mat[i],coeff(f,A[i]),ylist),i=1..nops(Umbra1Mat)):
end:






###multiple patterns
#SeqTuS(P,t,N): Like SeqT(P,t,N) but using Umbral methods.
#Given a set of patterns P, a variable t, and
# a pos. integer N, finds the first N terms of the
#sequence enumerating permutations according to the
#number of occurrences of the consecutive patterns P. With t=0
#you get the number of permutations avoiding the consecutive
#patterns of P. For example, try:
#SeqTuS({[1,2,3]},t,10);
SeqTuS:=proc(P,t,N) local gu,i1,x,z,f:


gu:=ClusterSeqS(P,t,N):
f:=1-z-add(gu[i1]*z^i1/i1!,i1=2..N):

f:=taylor(1/f,z=0,N+1):
[seq(expand(coeff(f,z,i1))*i1!,i1=1..N)]:
end:


#ClusterSeqS(P,t,N): Given a set of consecutive patterns
#P,uses the Umbral scheme to find the first N terms
#of the cluster-series for counting permutations
#according to the number of occurrences of members of P
#For example, try:
#ClusterSeqS([1,2,3],t,10);
ClusterSeqS:=proc(P,t,N) local U, Ini1,z,gu,mu,i1,ylist,i2,x,A:

if nops({seq(nops(P[i1]),i1=1..nops(P))})<>1 then
 print(`All members of the set of patterns, P (first argument)`):
 print(`must be of the same lenght`):
 RETURN(FAIL):
fi:

U:=SPtoUM(P,t,x,z,A):


ylist:=[seq(x[i1],i1=1..nops(P[1])),z]:


Ini1:=z^nops(P[1])*(t-1)*
mul(x[i1]^i1,i1=1..nops(P[1]))*
add(A[i1],i1=1..nops(P[1])):

mu:=Ini1:
gu:=z^nops(P[1])*(t-1)*nops(P):


for i1 from 1 to N do

mu:=ApplyUmbraMat(U,mu,ylist,A):
mu:=(add(expand(normal(coeff(mu,A[i2],1)))*A[i2],i2=1..nops(P))):


 mu:=add(coeff(mu,z,i2)*z^i2,i2=0..N):


 gu:=gu+subs({seq(A[i2]=1,i2=1..nops(P)),seq(x[i2]=1,i2=1..nops(P[1]))},mu):
od:
gu:=expand(gu):
[seq(coeff(gu,z,i2),i2=1..N)]:
end:




#SeqS(P,N): 
#Given a set of patterns P,   and
# a pos. integer N, finds the first N terms of the
#sequence enumerating permutations avoiding P
#For example, try:
#SeqS({[1,2,3]},10);
SeqS:=proc(P,N):SeqTuS(P,0,N):end:


#UpScheme(p): Given a SINGLE consecutive pattern p,
#and a symbol j to denote an indexed variable
#and a symbol n, produces an upward scheme
#listing the set of upwards templates
#that a monomial x[1]^j[1]*x[2]^j[2]*...*z^n can go to
#For example, try:
#UpScheme([1,2,3],j,n);
UpScheme:=proc(p,j,n) local gu,lu,lu1,ku1,i1:
lu:=Foll(p,{p}):
gu:={}:

for lu1 in lu do
 ku1:=TemplateUP(lu1,j):
 ku1:=[op(ku1),n+nops(p)-lu1[2]+1]:
gu:=gu union {ku1}:
od:
[gu,[seq(j[i1],i1=1..nops(p)),n]]:
end:



#ApplyUD11(T,mono,avar, xvar): Given a template T 
#a monomial mono, a list of discrete variables avar
#and a list of continuous variables xvar
#(with nops(T)=nops(avar)=nops(xvar))
#Applies the implied umbral operator given by T.
#For example, try:
#ApplyUD11([j1,0,j2,n+2],3*x1*x2^2*x3^3*z^3,[j1,j2,j3,n],[x1,x2,x3,z]);
ApplyUD11:=proc(T,mono,avar, xvar) local gu,i1,degs1,coef1,T1,gu1:

degs1:=[seq(degree(mono,xvar[i1]),i1=1..nops(xvar))]:
coef1:=normal(mono/mul(xvar[i1]^degs1[i1],i1=1..nops(xvar))):


T1:=subs({seq(avar[i1]=degs1[i1],i1=1..nops(avar))},T):


gu:=FillIn(T1):

expand(coef1*add(mul(xvar[i1]^gu1[i1],i1=1..nops(xvar)),gu1 in gu)):

end:


#ApplyUD1(U,mono, xvar): Given an upscheme U
#and a monomial in the list of variables xvar,
#applies the implied umbral operator U to the monomial mono
#For example, try:
#ApplyUD1([{[j1,0,j2,n+2]},[j1,j2,j3,j4],3*x1*x2^2*x3^3*z^3,[x1,x2,x3,z]);
ApplyUD1:=proc(U,mono, xvar) local avar,T1:
avar:=U[2]:
add(ApplyUD11(T1,mono, avar,xvar), T1 in U[1]):
end:



#ApplyUD(U,poly, xvar): Given an upscheme U
#and a polynomial, poly, in the list of variables xvar,
#applies the implied umbral operator U to the polynomial pol
#For example, try:
#ApplyUD([{[j1,0,j2,n+2]},[j1,j2,j3,j4],3*x1*x2^2*x3^3*z^3,[x1,x2,x3,z]);
ApplyUD:=proc(U,poly, xvar) local i1:

if not type(poly,`+`) then
 ApplyUD1(U,poly, xvar):
else
add( ApplyUD1(U,op(i1,poly), xvar),i1=1..nops(poly)):
fi:
end:







#ClusterSeqDir(p,t,N): Given a consecutive pattern
#uses the DIRECT Umbral scheme to find the first N terms
#of the cluster-series for counting permutations
#according to the number of occurrences of p.
#For example, try:
#ClusterSeqDir([1,2,3],t,10);
ClusterSeqDir:=proc(p,t,N) local U,j,x,i1,xvar,z,Ini,i2,n,gu,mu:

U:=UpScheme(p,j,n):

Ini:=mul(x[i1]^i1,i1=1..nops(p))*z^nops(p)*(t-1):

xvar:=[seq(x[i1],i1=1..nops(p)),z]:
mu:=Ini:
gu:=z^nops(p)*(t-1):


while ldegree(mu,z)<=N do 

mu:=expand(normal((t-1)*ApplyUD(U,mu, xvar)/z)):

 gu:=gu+subs({seq(x[i2]=1,i2=1..nops(p))},mu):
od:
gu:=expand(gu):
[seq(coeff(gu,z,i1),i1=1..N)]:
end:

#SeqW(FP,N,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:
#SeqW({[1,2,3]},40,10);
SeqW:=proc(FP,N,accur1)  local gu,mu1,mu2,mu3,A1,A2,A3,rho,A:
Digits:=3*accur1:
gu:=Seq(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):

 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:







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

gu:=convert(reps(k),list):



lu:=[]:

for i from 1 to nops(gu) do
lu:=[op(lu),[gu[i][1],Seq(gu[i][1],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:






#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:


#Reps(k,r): a set of representatives for the 
#sets of sets of r patterns all of length r
#For example, try: 
#Reps(3,2);
Reps:=proc(k,r) local gu,netsigim,gu1:
gu:=convert(permute(k),set):
gu:=choose(gu,r):
netsigim:={}:
while gu<>{} do
gu1:=gu[1]:
netsigim:=netsigim union {gu1}:
gu:=gu minus KHAVERIM(gu1):
od:
netsigim:
end:


#khaverim(p): all the friends of the pattern p.
#For example, try: khaverim([1,2,3,4]);
khaverim:=proc(p) local gu,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:



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

reps:=proc(k) local gu,mu,i:
mu:=convert(permute(k),set):
gu:={}:
while mu<>{} do
gu:=gu union {mu[1]}:
mu:=mu minus khaverim(mu[1]):
od:
gu:
end:



#FastestFriend(p,N): The fastest friend of the pattern p
#for producing the N first terms of the avoiding sequence
FastestFriend:=proc(p,N) local gu,t0,i,aluf, cand1, si:
gu:=khaverim(p):

aluf:=gu[1]:

t0:=time(): Seq(gu[1],N): si:=time()-t0:

for i from 2 to nops(gu) do

t0:=time(): Seq(gu[i],N): cand1:=time()-t0:

if cand1<si then
 si:=cand1:
 aluf:=gu[i]:
fi:
od:
aluf, si:

end:



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

gu:=convert(reps(k),list):



lu:=[]:

for i from 1 to nops(gu) do

p:=FastestFriend(gu[i],15)[1]:


lu:=[op(lu),[p,Seq(p,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:




#SipurOld(k,N,accur1,n): Old version of Sipur
#All the sequences enumerating
#permutations avoiding a single pattern of length k
#together with their guessed asymptotics to accuracy
#of accur1  digits. n is the symbol for the length
#For example, try:
#SipurOld(3,50,10,n);
SipurOld:=proc(k,N,accur1,n) local gu,lu,i:
gu:=AllSeqs(k,N):

print(`This is the story of all enumerating sequences for counting`):
print(`permutations avoiding a single CONSECUTIVE pattern of length`, k):
print(`We list the first`, N , `terms and the conjectured asymptotics`):

for i from 1 to nops(gu) do
print(`The equivalence class of patterns that ranked number`, i):
print(`is `, khaverim(gu[i][1]) ):
print(`The first`, N, `terms, starting at n=1, are`):
print(gu[i][2]):
lu:=SeqW(gu[i][1],N,accur1):
if lu[2]<>FAIL and lu[3]<>FAIL then
 print(`Asympotically it is`, n!*lu[3]*lu[2]^n):
fi:
od:

print(`The whole thing took`, time() , `seconds `):
end:



#Squeeze1(L): Given a list of pairs outputs a shorter list
#where those with equal second argument are lumped together
#For example, try: Squeeze1([[3,1],[5,1],[7,3],[9,3],[11,3],[15,5]]);
Squeeze1:=proc(L) local M,i:
M:=[[{L[1][1]},L[1][2]]]:

for i from 2 to nops(L) do
if L[i][2]=L[i-1][2] then
 M:=[op(1..nops(M)-1,M),[M[nops(M)][1] union {L[i][1]},
   M[nops(M)][2]]  ]:
 else
  M:=[op(M),[{L[i][1]},L[i][2]]]:
fi:
od:
M:
end:

#Sipur(k,N,accur1,n): 
#All the sequences enumerating
#permutations avoiding a single pattern of length k
#together with their guessed asymptotics to accuracy
#of accur1  digits. n is the symbol for the length
#For example, try:
#Sipur(3,50,10,n);
Sipur:=proc(k,N,accur1,n) local gu,lu,i,j:
gu:=AllSeqs(k,N):

print(`This is the story of all enumerating sequences for counting`):
print(`permutations avoiding a single CONSECUTIVE pattern of length`, k):
print(`We list the first`, N , `terms and the conjectured asymptotics`):
print(`There are`, nops(gu), `equivalence classes up to trivial`):
print(`consecutive-Wilf equivalence`):
gu:=Squeeze1(gu):
print(`But there are only `, nops(gu), `distinct sequences that show up`):

for i from 1 to nops(gu) do
if nops(gu[i][1])=1 then
print(`The single pattern that ranked number`, i, `is : `):
print(gu[i][1][1]):
print(`and of course its images`, khaverim(gu[i][1][1])):
else
print(`The set of single patterns that ranked number`, i, ` is: `):
print(gu[i][1]):
print(`and of course their images`):
print( {seq(khaverim(gu[i][1][j]),j=1..nops(gu[i][1])) }):
fi:
print(`The first`, N, `terms, starting at n=1, of the enumerating`):
print(`sequence for permutations that avoid this pattern(s)`):
print(gu[i][2]):
lu:=SeqW(gu[i][1][1],N,accur1):
if lu[2]<>FAIL and lu[3]<>FAIL then
 print(`Asympotically it is`, n!*lu[3]*lu[2]^n):
fi:
od:

print(`The whole thing took`, time() , `seconds `):
end:



#FEx1(U,var): inputs an umbral operator in the set of variables var
#and outputs the set variables that retain their exponent
#when applied to the monomial 
#var[1]*var[2]^2*...*var[nops(var)-1]^(nops(var)-1)*
#ar[nops(var)^(nops(var)-1)
#For example, try:
#FEx1(PtoU([3,2,1],0,x,z),[x[1],x[2],x[3],z]);
FEx1:=proc(U,var) local mono,i1,mu,gu,hafukh,c,h1,g1:
mono:=mul(var[i1]^i1,i1=1..nops(var)-1)*var[nops(var)]^(nops(var)-1):
mu:=expand(ApplyUmbra(U,mono,var)):

gu:={}:

for i1 from 1 to nops(var)-1 do
 if ldegree(mu,var[i1])=degree(mu,var[i1]) and degree(mu,var[i1])=i1 then
    gu:=gu union {i1}:
 fi:
od:
gu:

hafukh:={seq(i1,i1=1..nops(var)-1)} minus gu:

mono:=mul(var[g1]^g1, g1 in gu)*mul(var[h1]^c[h1],h1 in hafukh)*
var[nops(var)]^c[nops(var)]:

if
{seq(
evalb(ldegree(mono,var[g1])=degree(mono,var[g1])), g1 in gu),
seq(evalb(degree(mono,var[g1])=g1), g1 in gu)} ={true} then
RETURN(gu):
else
 RETURN({}):
fi:


end:

#SimplifyUmbra(Umb): Given an Umbra, collects like terms
#and returns a simplified version
SimplifyUmbra:=proc(Umb)
local gu,kv,i,mu,evar,F:
 
gu:=0:
kv:={}:
 
for i from 1 to nops(Umb) do
 gu:=gu+Umb[i][1]*F(Umb[i][2]):
 kv:=kv union {Umb[i][2]}:
od:

 
gu:=expand(gu):


mu:={}:
 
for i from 1 to nops(kv) do
evar:=op(i,kv):
mu:=mu union {[normal(coeff(gu,F(evar))),evar]}:
od:
 
mu:
 
end:

#ReduceU1(U1,var,kv): Given an umbral atom U1 [coeff,subs],
#and a list of variables var, and a subset, kv, of [1,nops(var)]
#finds the correponding umbral atom, and the correponding
#set of variables when one does the change of variables
#f(var)=mul(var[i1]^i1,i1 in kv)*g(var minus var[kv])
#For example, do
#ReduceU1(PtoU([2,3,1],0,x,z)[1],[x[1],x[2],x[3],z],{1,2});
ReduceU1:=proc(U1,var,kv) local i1,gu,mono,U12:
gu:={seq(var[i1]=U1[2][i1],i1=1..nops(U1[2]))}:
mono:=1:
for i1 from 1 to nops(var) do

if member(i1,kv) then
 mono:=mono*var[i1]^i1:
fi:
od:

mono:=normal(subs(gu,mono)/mono)*U1[1]:
U12:=[]:

for i1 from 1 to nops(var) do
 if not member(i1,kv) then
   U12:=[op(U12),U1[2][i1]]:
 fi:
od:

[mono,U12]     :
end:


#ReduceU(U,var,kv,mono): Given an umbral operator U {[coeff,subs]},
#and a list of variables var, and a subset, kv, of [1,nops(var)]
#and an initial monomial
#finds the correponding umbral operator, and the correponding
#set of variables when one does the change of variables
#and the modified initial monomial.
#f(var)=mul(var[i1]^i1,i1 in kv)*g(var minus var[kv])
#For example, do
#ReduceU(PtoU([2,3,1],0,x,z),[x[1],x[2],x[3],z],{1,2},-x[1]*x[2]^2*x[3]^3*z^3);
ReduceU:=proc(U,var,kv,mono) local u1,V,var1,i1:
var1:=[]:
 for i1 from 1 to nops(var) do
  if not member(i1,kv) then
   var1:=[op(var1),var[i1]]:
  fi:
 od:

V:={seq(ReduceU1(u1,var,kv),u1 in U)}:


V:=SimplifyUmbra(V):
[V,var1, subs({seq(var[i1]=1,i1 in kv)},mono)]:

end:




#PtoBUold(p,t,x,z): The best Umbral operators, 
#suitably reduced, amongts all images of p
#under trivial consecutive Wilf-equivalence.
#The output is a set of lists {[p1,U,mono,var]}
#where p1 is a pattern, U is the reduced umbral operator
#mono is the initial monomial, and var is the set
#of variables.
#For example, try:
#PtoBU([1,3,2],t,x,z);
PtoBUold:=proc(p,t,x,z) local alufim,si,mu,U,halevai,i,var,mono,kv,i1,gu,q:
option remember:
var:=[seq(x[i],i=1..nops(p)),z]:
mu:=khaverim(p):

U:=PtoU(mu[1],t,x,z):
si:=nops(FEx1(U,var)):
alufim:={mu[1]}:

for i from 2 to nops(mu) do
U:=PtoU(mu[i],t,x,z):
halevai:=nops(FEx1(U,var)):

 if halevai>si then
   alufim:={mu[i]}:
   si:=halevai:
 elif halevai=si then
  alufim:=alufim union {mu[i]}:
 fi:
od:

gu:={}:

mono:=(t-1)*mul(x[i1]^i1,i1=1..nops(p))*z^nops(p):

for q in alufim do 
 U:=PtoU(q,t,x,z):
 kv:=FEx1(U,var):
 U:=ReduceU(U,var,kv,mono):
 gu:=gu union {[q,U]}:
od:
gu:
end:

#BetterU1(U,ylist,mono): Given a functional operator U in ylist
#and a monomial mono,
#finds a better one with fewer variables
BetterU1:=proc(V) local i,U1,U,ylist,mono:
U:=V[1]: ylist:=V[2]: mono:=V[3]:
for i from 1 to nops(ylist)-1 do
U1:=RedVar(U,ylist,i):
if U1<>FAIL then
RETURN([U1,subs(ylist[i]=1,mono)] ):
fi:
od:
V:
end:


#BetterU(U): Given a functional operator U in ylist
#and a monomial mono representing a functional equation
#finds a better one with fewer variables
BetterU:=proc(U) local U1,U2,U3:
U1:=U:
U2:=BetterU1(U1):

while U1<>U2 do
U3:=BetterU1(U2):
U1:=U2: U2:=U3: 
od:
U2:
end:

#PtoBU(p,t,x,z): The best Umbral operators, 
#suitably reduced, amongts all images of p
#under trivial consecutive Wilf-equivalence.
#The output is a set of lists {[p1,U,mono,var]}
#where p1 is a pattern, U is the reduced umbral operator
#mono is the initial monomial, and var is the set
#of variables.
#For example, try:
#PtoBU([1,3,2],t,x,z);
PtoBU:=proc(p,t,x,z) local alufim,si,mu,U,halevai,i,var,mono,kv,i1,gu,q:
option remember:
var:=[seq(x[i],i=1..nops(p)),z]:
mu:=khaverim(p):

U:=PtoU(mu[1],t,x,z):
si:=nops(FEx1(U,var)):
alufim:={mu[1]}:

for i from 2 to nops(mu) do
U:=PtoU(mu[i],t,x,z):
halevai:=nops(FEx1(U,var)):

 if halevai>si then
   alufim:={mu[i]}:
   si:=halevai:
 elif halevai=si then
  alufim:=alufim union {mu[i]}:
 fi:
od:

gu:={}:

mono:=(t-1)*mul(x[i1]^i1,i1=1..nops(p))*z^nops(p):

for q in alufim do 
 U:=PtoU(q,t,x,z):
 kv:=FEx1(U,var):
 U:=ReduceU(U,var,kv,mono):
 U:=BetterU(U):
 gu:=gu union {[q,U]}:
od:
gu:
end:


   





#ClusterSeqF(p,t,N): Given a consecutive pattern
#uses the Umbral scheme to find the first N terms
#of the cluster-series for counting permutations
#according to the number of occurrences of p.
#For example, try:
#ClusterSeqF([1,2,3],t,10);
ClusterSeqF:=proc(p,t,N) local U, Ini1,z,gu,mu,i1,ylist,i2,U1,x,q:

U:=PtoBU(p,t,x,z)[1]:
q:=U[1]:
ylist:=U[2][2]:
Ini1:=U[2][3]:
U1:=U[2][1]:
mu:=Ini1:
gu:=z^nops(p)*(t-1):

for i1 from 1 to N do
 mu:=ApplyUmbra(U1,mu,ylist):
 mu:=add(coeff(mu,z,i2)*z^i2,i2=0..N):
 gu:=gu+subs({seq(x[i2]=1,i2=1..nops(p))},mu):
od:

gu:=expand(gu):
[seq(coeff(gu,z,i1),i1=1..N)],q:

end:


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

gu:=convert(reps(k),list):



lu:=[]:

for i from 1 to nops(gu) do

p:=FastestFriend(gu[i],15)[1]:


lu:=[op(lu),[p,Seq(p,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:


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

gu:=convert(reps(k),list):



lu:=[]:

for i from 1 to nops(gu) do
p:=PtoBU(gu[i],0,x,z)[1][1]:
lu:=[op(lu),[p,SeqF(p,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:


#SipurF(k,N,accur1,n): 
#All the sequences enumerating
#permutations avoiding a single pattern of length k
#together with their guessed asymptotics to accuracy
#of accur1  digits. n is the symbol for the length
#For example, try:
#SipurF(3,50,10,n);
SipurF:=proc(k,N,accur1,n) local gu,lu,i,j:
gu:=AllSeqsF(k,N):

print(`This is the story of all enumerating sequences for counting`):
print(`permutations avoiding a single CONSECUTIVE pattern of length`, k):
print(`We list the first`, N , `terms and the conjectured asymptotics`):
print(`There are`, nops(gu), `equivalence classes up to trivial`):
print(`consecutive-Wilf equivalence`):
gu:=Squeeze1(gu):
print(`But there are only `, nops(gu), `distinct sequences that show up`):

for i from 1 to nops(gu) do
if nops(gu[i][1])=1 then
print(`The single pattern that ranked number`, i, `is : `):
print(gu[i][1][1]):
print(`and of course its images`, khaverim(gu[i][1][1])):
else
print(`The set of single patterns that ranked number`, i, ` is: `):
print(gu[i][1]):
print(`and of course their images`):
print( {seq(khaverim(gu[i][1][j]),j=1..nops(gu[i][1])) }):
fi:
print(`The first`, N, `terms, starting at n=1, of the enumerating`):
print(`sequence for permutations that avoid this pattern(s)`):
print(gu[i][2]):
lu:=SeqW(gu[i][1][1],N,accur1):
if lu[2]<>FAIL and lu[3]<>FAIL then
 print(`Asympotically it is`, n!*lu[3]*lu[2]^n):
fi:
od:

print(`The whole thing took`, time() , `seconds `):
end:

#GuessPol(L,n,s0): guesses a polynomial of degree d in n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#GuessPol([(seq(i,i=1..10),n,1);
GuessPol:=proc(L,n,s0) local d,gu:

for d from 0 to nops(L)-s0-3 do
 gu:=GuessPol1(L,d,n,s0):
 if gu<>FAIL then
    RETURN(gu):
 fi:
od:

FAIL:

end:

#GuessPol1(L,d,n,s0): guesses a polynomial of degree d in n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#GuessPol1([(seq(i,i=1..10),1,n,1);
GuessPol1:=proc(L,d,n,s0) local P,i,a,eq,var:
if d>=nops(L)-s0-2 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=s0..s0+d+3)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

subs(var,P):

end:



#Moments(p,R,n): The  first R moments about the mean
#of the random variable #number of occurrences of pattern
#p in uniformly-at-random n-permutation
#For example try: 
#Moments([1,2,3],4,n);
Moments:=proc(p,R,n) local mu,lu,gu,r,j,vu,i,t:
r:=nops(p):
mu:=SeqTu(p,t,6*R+10+r):
mu:=[seq(mu[i]/t^((i-r+1)/r!)/i!,i=1..nops(mu))]:

gu:=[]:

for j from 1 to R do
 mu:=[seq(t*diff(mu[i],t),i=1..nops(mu))]:

 lu:=subs(t=1,mu):

 vu:=GuessPol(lu,n,r+3*j):

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

gu:

end:

#asympt1(f,x,s): the asymptotics to order s of f(x)
asympt1:=proc(f,x,s) local f1,i1:
f1:=normal(subs(x=1/x,f)):
f1:=taylor(f1,x=0,s+1):
add(coeff(f1,x,i1)/x^i1,i1=0..s):
end:

#AsyMoments(p,s,n): The  conjectured formula, to order s, for the
#normalized asymptotic moments (both even and odd)
#of the random variable #number of occurrences of pattern
#p in uniformly-at-random n-permutation
#For example try: 
#AsyMoments([1,2,3],s,n);
AsyMoments:=proc(p,s,n) local mu,R,luE,luO,i,s1:

R:=14*s+2:

mu:=Moments(p,R,n):

if nops(mu)<R then
 RETURN(FAIL):
fi:

luE:=[seq(asympt1(mu[2*i]/mu[2]^i/((2*i)!/i!/2^i),n,s+4),i=1..R/2)]:

#luO:=[seq(asympt1(simplify(mu[2*i-1]/mu[2]^(i-1/2))
#/((2*i)!/i!/2^i),n,s+4
#),i=1..R/2)]:



[seq([seq(coeff(luE[i],n,-s1),i=1..nops(luE))],s1=0..s)]:

#[seq([seq(op(s1+1,luO[i])*n^s1,i=1..nops(luO))],s1=0..s)]:


#to be continued
end:





#UtoD(U,var,mono,D1): Given a functional operator U, in the set of
#variables var, and initial monomial mono, tries to
#find an ordinary differential operator Op(z,D1)
#such that the solution to the functional equation
#f(var)=mono+Uf(var) upon plugging the catalytic variables in var
#to be all 1's except the last, that we call z,
#you get the differential equation for F(z)=f(1,1,1, ...,z)
#to be F(z)=Mono+Oper(z,D1), where Mono is mono with all
#the var's (except z) equal 1. 
#It returns the length-2 lisf of
#Oper(z,D1) followed by Mono if succesful,
#or FAIL, if not.
#For example, try:
#UtoD(op(PtoBU([3,2,1],t,x,z)[1][2]),D1);
UtoD:=proc(U,var,mono,D1) local gu,z,f,i,Mono,zub,mu:
z:=var[nops(var)]:
Mono:=mono:
gu:=ApplyUmbra(U,f(op(var)),var):
for i from 1 to nops(var)-1 do
 Mono:=subs(var[i]=1,Mono):
 gu:=limit(gu,var[i]=1):
od:

mu:=0:

zub:=f(1$(nops(var)-1),z):
mu:=coeff(gu,zub,1):
gu:=gu-mu*zub:

for i from 1 to 20 do
mu:=mu+coeff(gu,D[nops(var)$i](f)(1$(nops(var)-1),z),1)*D1^i:
gu:=gu-coeff(gu,D[nops(var)$i](f)(1$(nops(var)-1),z),1)
*D[nops(var)$i](f)(1$(nops(var)-1),z):
od:

gu:=expand(gu):

if gu<>0 then
 RETURN(FAIL):
else
 RETURN([mu,Mono]):
fi:


end:


#Mishpat(p,t): inputs a pattern, and outputs the theorem
#that enables fast (poly-time, hopefully optimal)
#computation of the enumerating sequence for the
#number of n-permutations avoiding the consecutive Wilf
#(or if t is a variable, the weight-enumerator according to number
#of occurrences of)
#pattern p. For example, try:
#Mishpat([1,2,3],t);
#Mishpat([1,2,3],0);
Mishpat:=proc(p,t) local D1, gu,lu,n,i,j,x,z,f,pi,a,c,mu,ku,var,Mono,U:

if not type(t, symbol) and t<>0 then
 print(`The second argument must be either 0 or a symbol`):
 RETURN(FAIL):
fi:

if t=0 then
lu:=khaverim(p):
gu:=PtoBU(p,0,x,z):
gu,lu:
print(`Theorem: Let a(n) be the number of permutations of length n`):
print(`avoiding the consecutive Wilf pattern`, p):
print(`In other words, the number of permutations`, pi):
print(`such that for all i between 1 and`, n-nops(p)+1):
print([seq(pi[i+j],j=0..nops(p)-1)], `does not reduce to `, p):
print(` Then `):
print(Sum(a(n)*z^n/n!,n=0..infinity)=(1-z-Sum(c(n)*z^n/n!,n=2..infinity))^(-1)):

gu:=PtoBU(p,0,x,z)[1]:
gu:=gu[2]:
mu:=UtoD(op(gu),D1):
if mu<>FAIL and degree(mu[1],D1)=0 then

ku:=solve(f-mu[2]-mu[1]*f,f):
print(`Where `):
print(Sum(c(n)*z^n,n=2..infinity)=ku):
elif mu<>FAIL and degree(mu[1],D1)>0 then
print(`Where g(z):=`):
print(Sum(c(n)*z^n,n=2..infinity)):
print(`is the (unique!) formal power series that`):
print(`satisfies the following ordinary differential equation`):
print(g(z)=mu[2]+coeff(mu[1],D1,0)*g(z)+add(diff(g(z),z$i)*coeff(mu[1],D1,i),i=1..degree(mu[1],D1))):
else
U:=gu[1]:
var:=gu[2]:
Mono:=gu[3]:
print(`Where f(z):=`):
print(Sum(c(n)*z^n,n=2..infinity)=g(1$(nops(var)-1),z)):
print(`and `, g(op(var)), `is the unique formal-power-series solution `):
print(`of the functional equation:`):
print(g(op(var))=Mono+ApplyUmbra(U,g(op(var)),var)):
fi:

elif type(t,symbol) then
lu:=khaverim(p):
gu:=PtoBU(p,t,x,z):
gu,lu:
print(`Theorem: Let a(n,t) be the polynomial in`, t, ` such that for any j,`):
print(`its coeff. of t^j  is the number of n-permutations`):
print(`with exacty j occurrences of the consecutive pattern`, p):
print(`In other words, the coeff. of t^j in a(n,t) is `):
print(`the number of n-permutations`, pi):
print(`that have exactly j places`):
print(`i, between 1 and`, n-nops(p)+1, `such that `):
print([seq(pi[i+j],j=0..nops(p)-1)], `reduces to `, p , `. `):
print(` Then: `):
print(Sum(a(n,t)*z^n/n!,n=0..infinity)=(1-z-Sum(c(n,t)*z^n/n!,n=2..infinity))^(-1)):

gu:=PtoBU(p,t,x,z)[1]:
gu:=gu[2]:
mu:=UtoD(op(gu),D1):
if mu<>FAIL and degree(mu[1],D1)=0 then

ku:=solve(f-mu[2]-mu[1]*f,f):
print(`Where `):
print(Sum(c(n,t)*z^n,n=2..infinity)=numer(ku)/sort(collect(denom(ku),z))):
elif mu<>FAIL and degree(mu[1],D1)>0 then
print(`Where g(z,t):=`):
print(Sum(c(n,t)*z^n,n=2..infinity)):
print(`is the (unique!) formal power series that`):
print(`satisfies the following ordinary differential equation`):
print(g(z,t)=mu[2]+coeff(mu[1],D1,0)*g(z,t)+add(diff(g(z,t),z$i)*coeff(mu[1],D1,i),i=1..degree(mu[1],D1))):
else
U:=gu[1]:
var:=gu[2]:
Mono:=gu[3]:
print(`Where f(z):=`):
print(Sum(c(n,t)*z^n,n=2..infinity)=g(1$(nops(var)-1),z,t)):
print(`and `, g(op(var),t), `is the unique formal-power-series solution `):
print(`of the functional equation:`):
print(g(op(var),t)=Mono+ApplyUmbra(U,g(op(var),t),var)):
fi:
fi:

end:



#Hochakha(p,t): inputs a pattern, and outputs the theorem 
#that enables fast (poly-time, hopefully optimal)
#computation of the enumerating sequence for the
#number of n-permutations avoiding the consecutive Wilf
#pattern p, and its PROOF.
#For example, try:
#Hochakha([1,2,3],t);
Hochakha:=proc(p,t) local p1,x,z,j,i,i1,lu,X,r1,k1,lu1,ku1,n,ku1a,gu,alist,
gu1,su1,xlist,mono,F,kv,kv1,D1,f,g,kv1a,ku,mu:

if not type(t, symbol) and t<>0 then
 print(`The second argument must be either 0 or a symbol`):
 RETURN(FAIL):
fi:

if t=0 then

Mishpat(p,t):
print(``):
p1:=PtoBU(p,t,x,z)[1][1]:

print(`Proof: `):
print(`It is easier to consider the equivalent problem of`):
print(`enumerating permutations avoiding the equivalent pattern`, p1):
print(`Recall that the generalized Goulden-Jackson method`):
print(`developed by Vladimir Dotsenko and Anton Khoroshkin, arXiv:1002.2761v1 [math.CO] 15 Feb 2010 `):
print(`(and independently by Brian Nakamura in his PhD thesis) tells us that`):
print(`if a(n) is the number of permutations avoiding the pattern`, p1, `then`):
print(Sum(a(n)*z^n/n!,n=0..infinity)=(1-z-Sum(c(n)*z^n/n!,n=2..infinity))^(-1)):
print(`where c(n) is the the sum of the weights of all clusters of length n`):
print(`where every block reduces to the pattern`, p1, `and the weight is`):
print((t-1), `to the power the number of blocks. `):
print(`So we are interested in the weight-enumerator`):
print(`(ordinary generating function) of c(n).`):
print(`We need to introduce`, nops(p), `catalytic variables`):
print(seq(x[i1],i1=1..nops(p))):
print(``):
print(`Define the weight of a cluster of length n `):
print(`where the sorted entries of the last block are`):
print([seq(j[i1],i1=1..nops(p))], `to be equal to ` ):
print(mul(x[i1]^j[i1],i1=1..nops(p))*z^n):
print(``):
lu:=Foll(p1,{p1}):
print(`Consider such a cluster, and let's see how we can extend it`):
print(`by adding legally another block`):
print(`Of course when we add new entries, we'll have to make room`):
print(`for those newcomers by relabelling the current entries.`):
print(`so that the labels of the extended cluster will be from 1 to n+e .`):
print(`where e is the number of new entries contributed by the added block.`):
print(`The possible overlaps between the penultimate block`):
print(`and the last block are at members of the`):
print(`following set of locations:`, {seq(lu[i][2],i=1..nops(lu))}):
print(`for the new block.`):
print(`Let `, [seq(j[i1],i1=1..nops(p))], ` be the sorted list`):
print(`of the letters of the current righmost block`):
print(`In other words, the last`, nops(p),` letters of the underlying`):
print(`permutation, in that order, are`):
print([seq(j[p1[i1]],i1=1..nops(p1))]):
print(`Let's see how we can add another block to the current last block`):
print(`Let `, [seq(i[i1],i1=1..nops(p))], ` be the sorted list`):
print(`of the letters of the new righmost block`):
print(`In other words, the last`, nops(p),` letters of the extended`):
print(`underlying permutation, in that order, are`):
print([seq(i[p1[i1]],i1=1..nops(p1))]):

gu:={}:
for  lu1 in lu do
print(`----------------------------`):
  k1:=lu1[2]:
 print(`if the new block has its `, k1, ` first entries touch`):
 print(`the last`, k1, ` entries of the existing block then`):
 print(seq(i[p1[r1]]=j[p1[nops(p)-k1+r1]],r1=1..k1)):
  ku1:=TemplateUP(lu1,j):
  ku1:=[op(ku1),n+nops(p)-lu1[2]+1]:
 print(`This means that the new (sorted)entries`, [seq(i[i1],i1=1..nops(p))]):
 print(`follow the template`):
 ku1a:=subs(0=X,ku1):
 print(ku1a):
 print(`where the places that are`):
 print(`not X are already committed, while the places with an X`):
 print(`denote ANY positive integer, provided that the resulting list`):
 print(`will be STRICTLY increasing.`):
 print(`So any specific cluster that ends with`):
 print([seq(j[p1[i1]],i1=1..nops(p1))]):
 print(`and hence with weight`):
 print(mul(x[i1]^j[i1],i1=1..nops(p))*z^n):
 print(`gives rise to the sum of`):
 print((t-1)*mul(x[i1]^i[i1],i1=1..nops(p))*z^(n+nops(p)-k1)):
 print(`over all the feasible increasing vectors`,[seq(i[i1],i1=1..nops(p))]):
 print(`satisfying the template`):
 print(ku1a):
 print(`This is a sum (or iterated sum) of geometrical series that turns out to be`):
 print(`the rational function`):

 xlist:=[seq(x[i1],i1=1..nops(p)),z]:
 su1:=(t-1)*schumT(ku1,xlist)/z:
 su1:=expand(su1):
 alist:=[seq(j[i1],i1=1..nops(p)),n]:
 print(su1):
 print(`By linearity this correponds to the functional operation`):
 alist:=[seq(j[i1],i1=1..nops(p)),n]:
 gu1:=ToUmbra(su1,xlist,alist):
gu1:={seq([gu1[i1][1],gu1[i1][3]],i1=1..nops(gu1))}:
 print(F(seq(x[i1],i1=1..nops(p)),z), `goes to `):
 print(SpellOutUmbra(gu1,F)):
  gu:=gu union  gu1:
od:


print(`Recal that the one-block cluster has weight`):
mono:=(t-1)*mul(x[i]^i,i=1..nops(p))*z^nops(p):
print(mono):

if nops(lu)>1 then
print(`Combining all the contributions from the different interfaces`):
fi:

print(`we get that the weight-enumerator of all clusters according to the pattern`,p):
print(`satisfies the functional equation`):
print(F(seq(x[i1],i1=1..nops(p)),z)=mono+SpellOutUmbra(gu,F)):

 kv:=FEx1(gu,xlist):
kv1:=convert({seq(i1,i1=1..nops(p))} minus kv,list):
kv1:=sort(kv1):

kv1a:=sort(convert(kv1,list)):
kv1a:=seq(x[kv1a[i1]],i1=1..nops(kv1a)):



if kv={seq(i1,i1=1..nops(p))} then
 print(`If the functional operator featured on the Right side is applied`):
 print(`to any formal power series of the form`):
 print(mul(x[i1]^i1,i1 in kv)*g(kv1a,z)):
 print(`it would still be of that form (with a different g, of course)`):
 print(`Hence we can simplify the above functional equation and put`):

 print(F(op(xlist))=mul(x[i1]^i1,i1 in kv)*g(kv1a,z)):

print(`Making this change of dependent variables yields the`):
print(`simple algebraic equation for`, g(kv1a,z)):
 gu:=ReduceU(gu,xlist,kv,mono):
 gu:=BetterU(gu):
print(g(op(gu[2]))=gu[3]+ApplyUmbra(gu[1],g(op(gu[2])),gu[2])):

gu:=PtoBU(p,t,x,z)[1]:
p1:=gu[1]:
gu:=gu[2]:
mu:=UtoD(op(gu),D1):
if not (mu<>FAIL and degree(mu[1],D1)=0) then
print(`Something is wrong`):
RETURN(FAIL):
fi:
ku:=solve(f-mu[2]-mu[1]*f,f):
print(`Solving this simple algebaric equation yields`):
print(g(z)=ku, `QED.`):
fi:

if kv<>{} and kv<>{seq(i1,i1=1..nops(p))} then
 print(`If the functional operator featured on the right side is applied`):
 print(`to any formal power series of the form`):
 print(mul(x[i1]^i1,i1 in kv)*g(kv1a,z)):
 print(`it would still be of that form (with a different g, of course)`):
 print(`Hence we can simplify the above functional equation, by setting:`):
 print(F(op(xlist))=mul(x[i1]^i1,i1 in kv)*g(kv1a,z)):


print(`Making this change of dependent variables yields the`):
print(`simplified functional equation for`, g(kv1a,z)):
 gu:=ReduceU(gu,xlist,kv,mono):
 gu:=BetterU(gu):
print(g(op(gu[2]))=gu[3]+ApplyUmbra(gu[1],g(op(gu[2])),gu[2])):


gu:=PtoBU(p,t,x,z)[1]:
p1:=gu[1]:
gu:=gu[2]:
mu:=UtoD(op(gu),D1):
if mu<>FAIL then
print(`Taking the limit as the variables`, kv1a, `go to 1`):
print(`and using L'Hospital's rule, yields the differential equation`):
print(g(z)=mu[2]+coeff(mu[1],D1,0)*g(z)+add(diff(g(z),z$i)*coeff(mu[1],D1,i),i=1..degree(mu[1],D1)),`QED.`):

fi:
fi:

fi:




if type(t,symbol) then

Mishpat(p,t):
print(``):
p1:=PtoBU(p,t,x,z)[1][1]:

print(`Proof: `):
print(`It is easier to consider the equivalent problem of`):
print(`enumerating permutations avoiding the equivalent pattern`, p1):
print(`Recall that the generalized Goulden-Jackson method`):
print(`developed by Vladimir Dotsenko and Anton Khoroshkin, arXiv:1002.2761v1 [math.CO] 15 Feb 2010 `):
print(`(and independently by Brian Nakamura in his PhD thesis) tells us that`):
print(`if a(n,t) is the weight-enumerator according to the number of`):
print(`occurrences of the pattern`, p1, `then`):
print(Sum(a(n,t)*z^n/n!,n=0..infinity)=(1-z-Sum(c(n,t)*z^n/n!,n=2..infinity))^(-1)):
print(`where c(n,t) is the the sum of the weights of all clusters of length n`):
print(`where every block reduces to the pattern`, p1, `and the weight is`):
print((t-1), `to the power the number of blocks. `):
print(`So we are interested in the weight-enumerator`):
print(`(ordinary generating function) of c(n,t).`):
print(`We need to introduce`, nops(p), `catalytic variables`):
print(seq(x[i1],i1=1..nops(p))):
print(``):
print(`Define the weight of a cluster of length n `):
print(`where the sorted entries of the last block are`):
print([seq(j[i1],i1=1..nops(p))], `to be equal to ` ):
print(mul(x[i1]^j[i1],i1=1..nops(p))*z^n):
print(``):
lu:=Foll(p1,{p1}):
print(`Consider such a cluster, and let's see how we can extend it`):
print(`by adding legally another block`):
print(`Of course when we add new entries, we'll have to make room`):
print(`for those newcomers by relabelling the current entries.`):
print(`so that the labels of the extended cluster will be from 1 to n+e .`):
print(`where e is the number of new entries contributed by the added block.`):
print(`The possible overlaps between the penultimate block`):
print(`and the last block are at members of the`):
print(`following set of locations:`, {seq(lu[i][2],i=1..nops(lu))}):
print(`for the new block.`):
print(`Let `, [seq(j[i1],i1=1..nops(p))], ` be the sorted list`):
print(`of the letters of the current righmost block`):
print(`In other words, the last`, nops(p),` letters of the underlying`):
print(`permutation, in that order, are`):
print([seq(j[p1[i1]],i1=1..nops(p1))]):
print(`Let's see how we can add another block to the current last block`):
print(`Let `, [seq(i[i1],i1=1..nops(p))], ` be the sorted list`):
print(`of the letters of the new righmost block`):
print(`In other words, the last`, nops(p),` letters of the extended`):
print(`underlying permutation, in that order, are`):
print([seq(i[p1[i1]],i1=1..nops(p1))]):

gu:={}:
for  lu1 in lu do
print(`----------------------------`):
  k1:=lu1[2]:
 print(`if the new block has its `, k1, ` first entries touch`):
 print(`the last`, k1, ` entries of the existing block then`):
 print(seq(i[p1[r1]]=j[p1[nops(p)-k1+r1]],r1=1..k1)):
  ku1:=TemplateUP(lu1,j):
  ku1:=[op(ku1),n+nops(p)-lu1[2]+1]:
 print(`This means that the new (sorted)entries`, [seq(i[i1],i1=1..nops(p))]):
 print(`follow the template`):
 ku1a:=subs(0=X,ku1):
 print(ku1a):
 print(`where the places that are`):
 print(`not X are already committed, while the places with an X`):
 print(`denote ANY positive integer, provided that the resulting list`):
 print(`will be STRICTLY increasing.`):
 print(`So any specific cluster that ends with`):
 print([seq(j[p1[i1]],i1=1..nops(p1))]):
 print(`and hence with weight`):
 print(mul(x[i1]^j[i1],i1=1..nops(p))*z^n):
 print(`gives rise to the sum of`):
 print((t-1)*mul(x[i1]^i[i1],i1=1..nops(p))*z^(n+nops(p)-k1)):
 print(`over all the feasible increasing vectors`,[seq(i[i1],i1=1..nops(p))]):
 print(`satisfying the template`):
 print(ku1a):
 print(`This is a sum (or iterated sum) of geometrical series that turns out to be`):
 print(`the rational function`):

 xlist:=[seq(x[i1],i1=1..nops(p)),z]:
 su1:=(t-1)*schumT(ku1,xlist)/z:
 su1:=expand(su1):
 alist:=[seq(j[i1],i1=1..nops(p)),n]:
 print(su1):
 print(`By linearity this correponds to the functional operation`):
 alist:=[seq(j[i1],i1=1..nops(p)),n]:
 gu1:=ToUmbra(su1,xlist,alist):
gu1:={seq([gu1[i1][1],gu1[i1][3]],i1=1..nops(gu1))}:
 print(F(seq(x[i1],i1=1..nops(p)),z), `goes to `):
 print(SpellOutUmbra(gu1,F)):
  gu:=gu union  gu1:
od:


print(`Recal that the one-block cluster has weight`):
mono:=(t-1)*mul(x[i]^i,i=1..nops(p))*z^nops(p):
print(mono):

if nops(lu)>1 then
print(`Combining all the contributions from the different interfaces`):
fi:

print(`we get that the weight-enumerator of all clusters according to the pattern`,p):
print(`satisfies the functional equation`):
print(F(seq(x[i1],i1=1..nops(p)),z)=mono+SpellOutUmbra(gu,F)):

 kv:=FEx1(gu,xlist):
kv1:=convert({seq(i1,i1=1..nops(p))} minus kv,list):
kv1:=sort(kv1):

kv1a:=sort(convert(kv1,list)):
kv1a:=seq(x[kv1a[i1]],i1=1..nops(kv1a)):



if kv={seq(i1,i1=1..nops(p))} then
 print(`If the functional operator featured on the Right side is applied`):
 print(`to any formal power series of the form`):
 print(mul(x[i1]^i1,i1 in kv)*g(kv1a,z,t)):
 print(`it would still be of that form (with a different g, of course)`):
 print(`Hence we can simplify the above functional equation and put`):

 print(F(op(xlist))=mul(x[i1]^i1,i1 in kv)*g(kv1a,z,t)):

print(`Making this change of dependent variables yields the`):
print(`simple algebraic equation for`, g(kv1a,z,t)):
 gu:=ReduceU(gu,xlist,kv,mono):
 gu:=BetterU(gu):
print(g(op(gu[2]),t)=gu[3]+ApplyUmbra(gu[1],g(op(gu[2]),t),gu[2])):

gu:=PtoBU(p,t,x,z)[1]:
p1:=gu[1]:
gu:=gu[2]:
mu:=UtoD(op(gu),D1):
if not (mu<>FAIL and degree(mu[1],D1)=0) then
print(`Something is wrong`):
RETURN(FAIL):
fi:
ku:=solve(f-mu[2]-mu[1]*f,f):
print(`Solving this simple algebaric equation yields`):
print(g(z,t)=ku, `QED.`):
fi:

if kv<>{} and kv<>{seq(i1,i1=1..nops(p))} then
 print(`If the functional operator featured on the right side is applied`):
 print(`to any formal power series of the form`):
 print(mul(x[i1]^i1,i1 in kv)*g(kv1a,z,t)):
 print(`it would still be of that form (with a different g, of course)`):
 print(`Hence we can simplify the above functional equation, by setting:`):
 print(F(op(xlist),t)=mul(x[i1]^i1,i1 in kv)*g(kv1a,z,t)):


print(`Making this change of dependent variables yields the`):
print(`simplified functional equation for`, g(kv1a,z,t)):
 gu:=ReduceU(gu,xlist,kv,mono):
 gu:=BetterU(gu):
print(g(op(gu[2]),t)=gu[3]+ApplyUmbra(gu[1],g(op(gu[2]),t),gu[2])):


gu:=PtoBU(p,t,x,z)[1]:
p1:=gu[1]:
gu:=gu[2]:
mu:=UtoD(op(gu),D1):
if mu<>FAIL then
print(`Taking the limit as the variables`, kv1a, `go to 1`):
print(`and using L'Hospital's rule, yields the differential equation`):
print(g(z,t)=mu[2]+coeff(mu[1],D1,0)*g(z,t)+add(diff(g(z,t),z$i)*coeff(mu[1],D1,i),i=1..degree(mu[1],D1)),`QED.`):

fi:
fi:

fi:



end:






#ClusterSeqFg(p,t,N,x): Given a consecutive pattern
#uses the Umbral scheme to find the first N terms
#of the generalized cluster-series for counting permutations
#according to the number of occurrences of p and the catalytic variables
#For example, try:
#ClusterSeqFg([1,2,3],t,10,x);
ClusterSeqFg:=proc(p,t,N,x) local U, Ini1,z,gu,mu,i1,ylist,i2,U1,q:

U:=PtoBU(p,t,x,z)[1]:
q:=U[1]:
ylist:=U[2][2]:
Ini1:=U[2][3]:
U1:=U[2][1]:
mu:=Ini1:
gu:=z^nops(p)*(t-1):

for i1 from 1 to N do
 mu:=ApplyUmbra(U1,mu,ylist):
 mu:=add(coeff(mu,z,i2)*z^i2,i2=0..N):
 gu:=gu+mu:
od:

gu:=expand(gu):
[seq(coeff(gu,z,i1),i1=1..N)],q:

end:





#SeferT(k,N,accur1,n,N1): 
#Inputs a pos. integer k, a pos. integer N, a pos. integer
#accur1, a symbol n, and a pos. integer N1, outputs
#a web-book with the statment of theorems (no proofs)
#that enable the optimal polynomial-time enumeration
#of the terms of the enumerating sequence for n-permutations
#avoiding single consecutive patterns of length k.
#N1 is a pos. integer less than N that tells how many
#terms to display in the book (too many terms could be
#tiring)
#For example, try:
#SeferT(3,50,10,n,15);
SeferT:=proc(k,N,accur1,n,N1) local gu,lu,i,j,co:
gu:=AllSeqs(k,N):
co:=0:
print(`Theorems on Enumerating Permutations Avoiding Single Consecutive Patterns of Length`, k):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This (computer-generated!) article contains all the theorems enabling`):
print(`efficient (poly-time and believably optimal) enumeration of`):
print(`permutations avoiding a single CONSECUTIVE pattern of length`, k):
print(`for all possible single patterns of that length.`):
print(`We also list the first`, N1 , `terms and the conjectured asymptotics`):
print(`There are`, nops(gu), `equivalence classes up to trivial`):
print(`consecutive-Wilf equivalence`):
gu:=Squeeze1(gu):
print(`But there are only `, nops(gu), `distinct sequences that show up`):

for i from 1 to nops(gu) do
if nops(gu[i][1])=1 then
print(`The single pattern that ranked number`, i, `is : `):
print(gu[i][1][1]):
print(`and of course its images`, khaverim(gu[i][1][1])):
else
print(`The set of single patterns that ranked number`, i, ` is: `):
print(gu[i][1]):
print(`and of course their images`):
print( {seq(khaverim(gu[i][1][j]),j=1..nops(gu[i][1])) }):
fi:
print(`The first`, N1, `terms, starting at n=1, of the enumerating`):
print(`sequence for permutations that avoid this pattern(s) are`):
print(op(1..N1,gu[i][2])):
lu:=SeqW(gu[i][1][1],N,accur1):
if lu[2]<>FAIL and lu[3]<>FAIL then
 print(`Asympotically it is`, n!*lu[3]*lu[2]^n):
fi:
print(`These terms (and as many more as you wish) could be efficiently`):
print(`computed from the following theorem, whose proof we know, but omit,`.``):
Mishpat(gu[i][1][1],0):

if nops(gu[i][1])>1 then
print(`The statements for the other patterns, that are`):
print(`non-trivially-Wilf equivalent to it are similar`):
print(`and imply (rigorously) that they are indeed Wilf-equivalent`):
fi:
print(`-------------------------------------------------`):
od:

print(`The whole thing took`, time() , `seconds `):
end:



#SeferP(k,N,accur1,n,N1): 
#Inputs a pos. integer k, a pos. integer N, a pos. integer
#accur1, a symbol n, and a pos. integer N1, outputs
#a web-book with the statment of theorems and proofs
#that enable the optimal polynomial-time enumeration
#of the terms of the enumerating sequence for n-permutations
#avoiding single consecutive patterns of length k.
#N1 is a pos. integer less than N that tells how many
#terms to display in the book (too many terms could be
#tiring)
#For example, try:
#SeferP(3,50,10,n,15);
SeferP:=proc(k,N,accur1,n,N1) local gu,lu,i,j,co:
gu:=AllSeqs(k,N):
co:=0:
print(`Theorems on Enumerating Permutations Avoiding Single Consecutive Patterns of Length`, k):
print(`By Shalosh B. Ekhad `):
print(``):
print(`This (computer-generated!) article contains all the theorems enabling`):
print(`efficient (poly-time and believably optimal) enumeration of`):
print(`permutations avoiding a single CONSECUTIVE pattern of length`, k):
print(`for all possible single patterns of that length.`):
print(`We also list the first`, N1 , `terms and the conjectured asymptotics`):
print(`There are`, nops(gu), `equivalence classes up to trivial`):
print(`consecutive-Wilf equivalence`):
gu:=Squeeze1(gu):
print(`But there are only `, nops(gu), `distinct sequences that show up`):

for i from 1 to nops(gu) do
if nops(gu[i][1])=1 then
print(`The single pattern that ranked number`, i, `is : `):
print(gu[i][1][1]):
print(`and of course its images`, khaverim(gu[i][1][1])):
else
print(`The set of single patterns that ranked number`, i, ` is: `):
print(gu[i][1]):
print(`and of course their images`):
print( {seq(khaverim(gu[i][1][j]),j=1..nops(gu[i][1])) }):
fi:
print(`The first`, N1, `terms, starting at n=1, of the enumerating`):
print(`sequence for permutations that avoid this pattern(s) are`):
print(op(1..N1,gu[i][2])):
lu:=SeqW(gu[i][1][1],N,accur1):
if lu[2]<>FAIL and lu[3]<>FAIL then
 print(`Asympotically it is`, n!*lu[3]*lu[2]^n):
fi:
print(`These terms (and as many more as you wish) could be efficiently`):
print(`computed from the following theorem, whose proof follows`.``):
Hochakha(gu[i][1][1],0):

if nops(gu[i][1])>1 then
print(`The statements for the other patterns, that are`):
print(`non-trivially-Wilf equivalent to it are similar`):
fi:
print(`-------------------------------------------------`):
od:

print(`The whole thing took`, time() , `seconds `):
end: