###################################################################### ##WILF: Save this file as WILF. To use it, stay in the # ##same directory, get into Maple (by typing: maple ) # ##and then type: read WILF : # ##Then follow the instructions given there # ## # ##Written by Doron Zeilberger, Temple University , # #zeilberg@math.temple.edu. # ####################################################################### #Created: April 23 , 1998 #Prevuius version: May 27, 2003 [Thanks to Vince Vatter!] #This version: June 24, 2004 [Thanks to Vince Vatter!, again!!] #WILF: A Maple package to study Wilf classes #Please report bugs to zeilberg@math.temple.edu with(combinat): print(`Added July 19, 2004: A much more general program is now`): print(`available from Vince Vatter http://www.math.rutgers.edu/~vatter/`): print(`called WILFPLUS `): print(``): print(`Created: April 23, 1998.`): print(`Previous version: May 27, 2003 [thanks to Vince Vatter]`): print(`This version: June 24, 2004 [thanks to Vince Vatter, again!]`): print(`Thanks to VINCE VATTER for pointing out that the previous version`): print(` of SchemeF did not capture the full diehedral group!`): lprint(``): print(`Written by Doron Zeilberger, zeilberg@math.rutgers.edu`): print(`Current address:, zeilberg@math.rutgers.edu`): lprint(``): print(`This is WILF, named in honor of Herb Wilf,`): print(`to study Wilf classes of permutations`): print(`It accompanies the paper `): print(`Enumeration Schemes, and more importantly, their Automaic Generation`): print(`by Doron Zeilberger, published in Annals of Combinatorics`): print(`and also available from his website`): lprint(``): print(`Please report bugs to zeilberg@math.rutgers.edu`): lprint(``): print(`The most current version of this package and paper`): print(` are available from`): print(`http://www.math.rutgers.edu/~zeilberg`): print(`For a list of the procedures type ezra(), for help with`): print(`a specific procedure, type ezra(procedure_name)`): print(``): ezra:=proc() if args=NULL then print(`Contains the following procedures: `): print(` `): print(``): print(` FindiJ, FindSchemes, funyot `): print(`good, good1, Insert, Julian, Miklos, perush, perush1, PScheme `): print(`redu, Refinements, Scheme, SchemeF, SchemeFv, SchemeRecurrence`): print(`SchemeSeries, `): print(` SchemeU,Sim1, Singles, Wilf, Zinn `): fi: if nops([args])=1 and op(1,[args])=FindSchemes then print(`FindSchemes(Resh,a,Depth,L): Given a list of integers Resh, and an`): print(`integers a, Depth and L, finds the set of all permutation`): print(`sets with Resh[i] permutations of length a+i-1 (i=1..nops(Resh))`): print(`that have a scheme of depth<=Depth, followed by their scheme, and`): print(`the first L terms. For example, try:`): print(`FindSchemes([1],3,2,10);`): fi: if nops([args])=1 and op(1,[args])=`SchemeRecurrence` then print(`SchemeRecurrence(schem,n,N,L): Given a scheme schem, finds`): print(`empirically a linear recurrence equation with`): print(`polynomial coefficients (where n is the discrete argument`): print(`and N is the shift operator in n) satisfied by it, by trying`): print(`up to L terms. If unsuccesful, it returns the sequence`): fi: if nops([args])=1 and op(1,[args])=`Zinn` then print(` Zinn's method `): print(`Given an increasing sequence a(n) of positive integers , expressed`): print(`in terms of a list, estimates the theta and mu such that `): print(` a(n) is asympt. to n^(theta)*mu^n. The output is `): print(` theta, mu `): fi: if nops([args])=1 and op(1,[args])=SchemeSeries then print(`SchemeSeries(schem,L): Given a scheme schem, and an integer`): print(`L, computes the first L terms of the`): print(`sequence generated by schem`): fi: if nops([args])=1 and op(1,[args])=SchemeU then print(`SchemeU(Patterns,GVUL): Given a set of permutations `): print(`(of various lengths)`): print(`Patterns, finds an Upper scheme [Redu,Expa,A,B,C] where`): print(`Redu is a set of reducible prefix-permutations with respect to`): print(`Patterns, Expa is the set of expandible ones, A is a table `): print(`where A[sigma] is the place (between 1 and nops(sigma)) that says`): print(`where to delete, for sigma in Redu, `): print(`and C[sigma] is a table that says the Findj of sigma`): print(`[i.e. the restriction, if any on the prefixes]`): print(`and B is a table `): print(`where B[sigma] is the set of refinements of sigma, of length one`): print(`more. For this to be a scheme, the deleted versions of all`): print(`sigma in Redu, and all the B[sigma] must be good`): print(`GVUL is the maximal permitted prefix size on either`): print(`Redu or Expa,`): print(`The Upper Scheme means that the sequence computed by Miklos `): print(`using the scheme is`): print(`an UPPER bound for the actual sequence of cardinalities`): print(`of Wilf(n,Patterns)`): print(`Once again, the syntax is:SchemeU(Patterns,GVUL)`): fi: if nops([args])=1 and op(1,[args])=PScheme then print(`PScheme(Patterns,GVUL): Given a set of permutations`): print(`(of various lengths)`): print(`Patterns, finds a Partial scheme [Redu,Expa,A,B,C] where`): print(`Redu is a set of reducible prefix-permutations with respect to`): print(`Patterns, Expa is the set of expandible ones, A is a table `): print(`where A[sigma] is the place (between 1 and nops(sigma)) that says`): print(`where to delete, for sigma in Redu, `): print(`and C[sigma] is a table that says the Findj of sigma`): print(`[i.e. the restriction, if any on the prefixes]`): print(`and B is a table `): print(`where B[sigma] is the set of refinements of sigma, of length one`): print(`more. For this to be a scheme, the deleted versions of all`): print(`sigma in Redu, and all the B[sigma] must be good`): print(`GVUL is the maximal permitted prefix size on either`): print(`Redu or Expa,`): print(`If a complete scheme is found, it returns it, otherwise`): print(`it returns 0 followed by the partial scheme`): fi: if nops([args])=1 and op(1,[args])=funyot then print(`funyot(tau,Sch,x,L): Given a scheme, Sch, and a prefix-permutation`): print(`tau, a letter x, and an intger L, finds the generating function`): print(`sum(Julian(n,tau,vect,Sch)*x[1]^(vect[tauh[1])*...`): print(`*x[k]^(vect[tauh[k])`): print(`*x[k+1]^n, over all conceivable vectors vect with max(vect)<=n<=L`): fi: if nops([args])=1 and op(1,[args])=Sim1 then print(`Sim1(k,L): All vectors of integers with k components`): print(`[v_1,...,v_k] such that 1<=v_1k or min(op(J))<0 then ERROR(`Bad input`): fi: if member(0,J) and w[1]<>1 then RETURN(false): fi: if member(k,J) and w[k]<>n then RETURN(false): fi: for j from 1 to k-1 do if member(j,J) and w[j+1]-w[j]<>1 then RETURN(false): fi: od: true: end: #IV1(n,k,J) all increasing vectors of length #k of the form 1<=i_1< ...n then RETURN(1): fi: gu:=IV(n,k-2): for i from 1 to nops(gu) do vec:=op(i,gu): subperm:=[op(1,pi)]: for j from 1 to k-2 do subperm:=[op(subperm),pi[vec[j]]]: od: subperm:=[op(subperm),pi[n]]: if redu(subperm)=pattern1 then RETURN(0): fi: od: 1: end: #good(pi,SetPatterns) #Given a permutation pi, and a set of patterns #SetPatterns, decides whether #pattern1 occurs in pi with the first and last letter of #pi coinciding #with the last letter of pattern1. It returns 0 if this is #the case, otherwise 1 (0 means bad) good:=proc(pi,SetPatterns) local i: if member(pi,SetPatterns) then RETURN(0): fi: for i from 1 to nops(SetPatterns) do if good1(pi,op(i,SetPatterns))=0 then RETURN(0): fi: od: 1: end: #Insert(pi,i): Given a permutation pi, and an integer #i between 1 and nops(pi)+1, constructs the permutation #of length nops(pi)+1 whose last entry is i, and such that #reduced form of the first nops(pi) letters is pi Insert:=proc(pi,i) local mu,j: mu:=[]: for j from 1 to nops(pi) do if pi[j]nops(vect) then ERROR(`the 2nd and 3rd args must be of same length`): fi: Sch:=Scheme(Patterns,GVUL): if Sch=0 then RETURN(0): fi: if n<=nops(tau) then if member(tau,Wilf(n,Patterns)) then RETURN(1): else RETURN(0): fi: fi: if redu(vect)<>tau then ERROR(`the reductions of the third argument must be the second arg.`): fi: if vect<>[] then if min(op(vect))<1 or min(op(vect))>n or nops(vect)<>nops(convert(vect,set)) then ERROR(`The third argument is bad`): fi: fi: Redu:=Sch[1]: Expa:=Sch[2]: A:=Sch[3]: B:=Sch[4]: C:=Sch[5]: if member(tau,Redu) then i:=A[tau]: J:=C[tau]: if member(0,J) then if not vect[tauh[1]]=1 then RETURN(0): fi: fi: if member(k,J) then if not vect[tauh[k]]=n then RETURN(0): fi: fi: for j from 1 to k do if member(j,J) then if j>0 and j{} do KV:=KV union kha: kha:=NewRelations(KV): od: KV: end: #kosherj(sigma,tau,i,j): Given a prefix permutation sigma #and a pattern permutation tau, and a set of places in #[1,nops(sigma)], and an integer j between 0 and nops(sigma) # indicating the restrictions #on the permutations x whose first nops(sigma) entries #reduce to sigma, decides whether the event i is possible #in light of the J restrictions kosherj:=proc(sigma,tau,i,j) local gu2,sigmah,k,mu,i1,mu1: k:=nops(sigma): if not(j>=0 and j<=nops(sigma)) then ERRROR(`Last argument, j, is out of range`): fi: sigmah:=hofkhi(sigma): gu2:=perush1(tau,x,y,i): if j=0 then if member(x[sigmah[1]],convert(gu2,set)) and not op(1,gu2)=x[sigmah[1]] then RETURN(false): else RETURN(true): fi: fi: if j=k then if member(x[sigmah[k]],convert(gu2,set)) and not op(nops(gu2),gu2)=x[sigmah[k]] then RETURN(false): else RETURN(true): fi: fi: if member(x[sigmah[j]],convert(gu2,set)) and member(x[sigmah[j+1]],convert(gu2,set)) then mu:={}: for i1 from 1 to nops(gu2)-1 do mu:=mu union {[gu2[i1],gu2[i1+1]]}: od: mu1:=Closure(mu) minus mu: if member([x[sigmah[j]],x[sigmah[j+1]]],mu1) then RETURN(false): else RETURN(true): fi: fi: true: end: #kosher(sigma,tau,i,J): Given a prefix permutation sigma #and a pattern permutation tau, and a set of places i, in #[1,nops(sigma)], and a set of integers J between 0 and nops(sigma) # indicating the restrictions #on the permutations x whose first nops(sigma) entries #reduce to sigma, decides whether the event i is possible #in light of the J restrictions kosher:=proc(sigma,tau,i,J) local j: for j from 1 to nops(J) do if not kosherj(sigma,tau,i,op(j,J)) then RETURN(false): fi: od: true: end: #ImpliedRelations(sigma,tau,i): Given a permutation #sigma( a prefix) of length N, say, and a #permutation tau (a pattern) of length M, say #and a set of places i (nops(i)=I0) then ERROR(`Bad input`): fi: sigma1:=[seq(sigma[i[k1]],k1=1..nops(i))]: tau1:=[op(1..nops(i),tau)]: if not redu(sigma1)=redu(tau1) then ERROR(`Bad input for sigma and tau`): fi: gu1:=perush(sigma,x): gu2:=perush1(tau,x,y,i): kv:={}: for i1 from 1 to nops(gu1) do for j1 from i1+1 to nops(gu1) do kv:=kv union {[op(i1,gu1),op(j1,gu1)]}: od: od: for i1 from 1 to nops(gu2) do for j1 from i1+1 to nops(gu2) do kv:=kv union {[op(i1,gu2),op(j1,gu2)]}: od: od: Closure(kv): end: #iImpliesi1(sigma,tau,i,i1): Given a prefix-permutation #sigma, and a pattern permutation tau, and #two set of places i and i1, decides whether #the existence of the pattern tau where the first #nops(i) participants of the pattern are at #the places indicated by i, implies the #existence of a pattern with the same #phenomenon with i1 iImpliesi1:=proc(sigma,tau,i,i1) local kv,kv1: kv:=ImpliedRelations(sigma,tau,i): kv1:=ImpliedRelations(sigma,tau,i1): if kv1 minus kv={} then RETURN(true): else RETURN(false): fi: end: #Events1J(sigma,tau,J): Given a prefix permutation sigma #and a pattern permutation tau , and a subset J of [0,nops(sigma)] #indicating the mandatory twins #finds all the subsets #of places of [1,nops(sigma)] that produce `events's i.e. #[sigma[i_1], ..., sigma[i_r]] reduces to a prefix of #tau Events1J:=proc(sigma,tau,J) local kv,KV,i,j,Pref_tau,kulam,mu: Pref_tau:={}: for i from 1 to nops(tau)-1 do Pref_tau:=Pref_tau union {redu([op(1..i,tau)])}: od: KV:={}: kulam:=powerset(nops(sigma)) minus {{}}: for i from 1 to nops(kulam) do kv:=op(i,kulam): mu:=sort([op(kv)]): if member(redu([seq(sigma[mu[j]],j=1..nops(mu))]),Pref_tau) then if kosher(sigma,tau,kv,J) then KV:=KV union {kv}: fi: fi: od: KV: end: #EventsJ(sigma,Patterns,J): #Given a prefix permutation sigma, and a set of #pattern permutation, Patterns, finds all #possible events , subsets of [1,nops(sigma)] with #the accompanying tau of Patterns that make it works #also given is a subset J of [0,nops(sigma)] # #The output is a set of ordered pairs #compatible with J: (recall that a pair #is a list of two elements), the first member of each such pair #is the subset of [1,nops(sigma)], the second member is #the successful permutation tau from Patterns EventsJ:=proc(sigma,Patterns,J) local i,tau,kv,KV,j: KV:={}: for i from 1 to nops(Patterns) do tau:=op(i,Patterns): kv:=Events1J(sigma,tau,J): for j from 1 to nops(kv) do KV:=KV union {[op(j,kv),tau]}: od: od: KV: end: #E1ImpliesE2(sigma,E1,E2): Given a prefix-permutation #sigma, and two events E1 and E2 of the form #[i1,tau1], [i2, tau2] #where i and i1 are sets of places and tau1 and tau2 #are pattern permutations decides whether #the existence of the pattern tau1 where the first #nops(i1) participants of the pattern are at #the places indicated by i1, implies the #existence of a pattern tau2 with the same #phenomenon with i2 #At present we assume nops(i1)=nops(i2) E1ImpliesE2:=proc(sigma,E1,E2) local i1,tau1,i2,tau2,kv1,kv2: i1:=E1[1]: tau1:=E1[2]: i2:=E2[1]: tau2:=E2[2]: kv1:=ImpliedRelations(sigma,tau1,i1): kv2:=ImpliedRelations(sigma,tau2,i2): if kv2 minus kv1={} then RETURN(true): else RETURN(false): fi: end: #BailOut(sigma,E1,SetE2): Given a prefix permutation #sigma,an event E1, and a set of events SetE2, decides #whether there exists at least one member of SetE2, let's #call it E2, such that E1ImpliesE2(sigma,E1,E2) is true BailOut:=proc(sigma,E1,SetE2) local E2,i: for i from 1 to nops(SetE2) do E2:=op(i,SetE2): if E1ImpliesE2(sigma,E1,E2) then RETURN(true): fi: od: false: end: #TestiJ:=proc(sigma,Patterns,i,J) #Given a prefix permutation, sigma, and a set of #patterns, Patterns, and an integer in the closed #interval between 1 and nops(sigma), and a set #of integers J (a subset of [0,nops(sigma)], #decides whether inserting #the entry at the i^th place is `harmless', i.e. #any patterns that it might participate in imply #previous patterns where he did not participate TestiJ:=proc(sigma,Patterns,i,J) local EventsWithi, EventsWithouti,j,AllEvents,eve: EventsWithi:={}: EventsWithouti:={}: AllEvents:=EventsJ(sigma,Patterns,J): for j from 1 to nops(AllEvents) do eve:=op(j,AllEvents): if member(i,eve[1]) then EventsWithi:=EventsWithi union {eve}: else EventsWithouti:=EventsWithouti union {eve}: fi: od: for j from 1 to nops(EventsWithi) do eve:=op(j,EventsWithi): if not BailOut(sigma,eve,EventsWithouti) then RETURN(false): fi: od: true: end: #FindiJ(sigma,Patterns,J): #Given a prefix permutation sigma, a set of patterns #Patterns, and a set of restrictions J, finds the set of #removable i's FindiJ:=proc(sigma,Patterns,J) local kv,i: kv:={}: for i from 1 to nops(sigma) do if TestiJ(sigma,Patterns,i,J) then kv:=kv union {i}: fi: od: kv: end: #Refinements(permu,Patterns,J): All the refinements of permu # (i.e. permutations of length one more, s.t. chopping #the last letter reduces to permu) that avoid the #patterns in Patterns, honoring the restrictions given by J Refinements:=proc(permu,Patterns,J) local i,j,lu,gu: gu:={}: for i from 1 to nops(permu)+1 do if not member(i-1,J) then lu:=[]: for j from 1 to nops(permu) do if op(j,permu)0 then print(`The set of permutations`): print(Patterns): print(`has a scheme of depth`,i): print(gu[1],gu[2],gu[3],op(gu[4]),op(gu[5])): RETURN([Patterns,i,gu]): fi: for j from 1 to nops(khaver) do gu:=Scheme(op(j,khaver),i): if gu<>0 then print(`The set of permutations`): print(Patterns): print(`does not have a scheme of depth`,i): print(`but its following image does`): print(op(j,khaver)): print(`The scheme is`): print(gu[1],gu[2],gu[3],op(gu[4]),op(gu[5])): RETURN([op(j,khaver),i,gu]): fi: od: od: 0: end: #SchemeF(Patterns,GVUL): Given a set of permutations #(of various lengths) #Patterns, finds a scheme for Patterns, #of depth <=GVUL, or #a scheme for one of its images under OPPOS , REVERS or OPPOS(REVERS) #of minimal depth. #If such a scheme exist, the program would find it (time and memory #permitting, if there is no such scheme, it returns 0 SchemeF:=proc(Patterns,GVUL) local i,j,gu,khaver: khaver:=KHAVERIM(Patterns) minus {Patterns}: for i from 1 to GVUL do gu:=Scheme(Patterns,i): if gu<>0 then RETURN([Patterns,i,gu]): fi: for j from 1 to nops(khaver) do gu:=Scheme(op(j,khaver),i): if gu<>0 then RETURN([op(j,khaver),i,gu]): fi: od: od: 0: end: #Singles(k,GVUL,L): Tries to find schemes for all Wilf #classes for all single-permutations of k letter #up to the trivial equivalence. #it prints out the succesful permutation from its class, #the scheme, and the list of cardinalities for n<=L #It returns a set of lists where each list has three #entries: the succesful permutation, the scheme, and #the list of cardinalities SinglesV:=proc(k,GVUL,L) local mu,gu,per,ku,lu,i,i1,Sch: gu:=convert(permute(k),set): mu:={}: while gu<>{} do per:=op(1,gu): mu:=mu union {khaverim(per)}: gu:=gu minus khaverim(per): od: ku:={}: for i from 1 to nops(mu) do gu:=SchemeF({op(1,op(i,mu))},GVUL): if gu=0 then print(`No scheme of depth<=`,GVUL,`found for the Wilf classes `): print(`in the trivial equivalence class`): print(op(i,mu)): else print(`One of the member-permutation in the equivalence class`): print(op(i,mu)): print(` in other words `, op(gu[1])): print(`has the following Scheme of depth`, gu[2]): Sch:=[gu[3][1],gu[3][2],gu[3][3],gu[3][4],gu[3][5]]: print([Sch[1],Sch[2],op(Sch[3]),op(Sch[4]),op(Sch[5])]): lu:=[seq(Miklos(i1,Sch),i1=1..L)]: print(`The cardinalities are`): print(lu): ku:=ku union {[op(1,gu[1]),gu[2],gu[3],lu]}: fi: od: ku: end: Singles:=proc(k,GVUL,L) local mu,gu,per,ku,lu,i,i1,Sch: gu:=convert(permute(k),set): mu:={}: while gu<>{} do per:=op(1,gu): mu:=mu union {per}: ku:=KHAVERIM({per}): ku:={seq(ku[i][1],i=1..nops(ku))}: gu:=gu minus ku: od: ku:={}: for i from 1 to nops(mu) do gu:=SchemeF({mu[i]},GVUL): if gu<>0 then Sch:=[gu[3][1],gu[3][2],gu[3][3],gu[3][4],gu[3][5]]: lu:=[seq(Miklos(i1,Sch),i1=1..L)]: ku:=ku union {[op(1,gu[1]),gu[2],gu[3],lu]}: fi: od: ku: end: #Scheme(Patterns,GVUL): Given a set of permutations (of various lengths) #Patterns, finds a scheme [Redu,Expa,A,B,C] where #Redu is a set of reducible prefix-permutations with respect to #Patterns, Expa is the set of expandible ones, A is a table #where A[sigma] is the place (between 1 and nops(sigma)) that says #where to delete, for sigma in Redu, #and C[sigma] is a table that says the Findj of sigma #[i.e. the restriction, if any on the prefixes] #and B is a table #where B[sigma] is the set of refinements of sigma, of length one #more. For this to be a scheme, the deleted versions of all #sigma in Redu, and all the B[sigma] must be "good" #GVUL is the maximal permitted prefix size on either #Redu or Expa, #If unsucessul, it returns 0 Scheme:=proc(Patterns,GVUL) local Redu,Expa,A,B,C,i,KVU,hakhi,permu,permud,I1,J: option remember: if member([],Patterns) then ERROR(`[] can't be a member of Patterns`): fi: hakhi:=0: KVU:={[]}: Redu:={}: Expa:={}: while KVU<>{} and hakhi<=GVUL do permu:=op(1,KVU): if nops(permu)>hakhi then hakhi:=nops(permu): fi: if hakhi>GVUL then RETURN(0): fi: J:=Findj(permu,Patterns): I1:=FindiJ(permu,Patterns,J): if I1={} then Expa:=Expa union {permu}: B[permu]:=Refinements(permu,Patterns,J): C[permu]:=J: KVU:=KVU minus {permu}: KVU:=KVU union (B[permu] minus (Expa union Redu)): fi: if I1<>{} then Redu:=Redu union {permu}: i:=op(1,I1): A[permu]:=i: C[permu]:=J: permud:=[op(1..i-1,permu),op(i+1..nops(permu),permu)]: permud:=redu(permud): KVU:=KVU minus {permu}: if not member(permud,Redu union Expa) then KVU:=KVU union {permud}: fi: fi: od: [Redu,Expa,A,B,C]: end: #Julian(n,tau,vect,Sch): Given an integer n, a #permutation tau, a vector, vect, of distinct integers #<=n, and a Scheme uses the Scheme #to compute the number of permutations of #length n that start with vect, reduce to tau, Julian:=proc(n,tau,vect,Sch) local gu,Redu,Expa,A,B,C,i,j,k,vect1,vect1n,vi,j1,Refs,x,tau1,tauh,J: option remember: if Sch=0 then ERROR(`Scheme is 0`): fi: k:=nops(tau): tauh:=hofkhi(tau): if nops(tau)<>nops(vect) then ERROR(`the 2nd and 3rd args must be of same length`): fi: if redu(vect)<>tau then ERROR(`the reductions of the third argument must be the second arg.`): fi: if vect<>[] then if min(op(vect))<1 or min(op(vect))>n or nops(vect)<>nops(convert(vect,set)) then ERROR(`The third argument is bad`): fi: fi: if n<=nops(tau) then RETURN(1): fi: Redu:=Sch[1]: Expa:=Sch[2]: A:=Sch[3]: B:=Sch[4]: C:=Sch[5]: if member(tau,Redu) then i:=A[tau]: J:=C[tau]: if member(0,J) then if not vect[tauh[1]]=1 then RETURN(0): fi: fi: if member(k,J) then if not vect[tauh[k]]=n then RETURN(0): fi: fi: for j from 1 to k do if member(j,J) then if j>0 and j{} and hakhi<=GVUL do permu:=op(1,KVU): if nops(permu)>hakhi then hakhi:=nops(permu): fi: if hakhi>GVUL then RETURN([0,KVU,[Redu,Expa,A,B,C]]): fi: J:=Findj(permu,Patterns): I1:=FindiJ(permu,Patterns,J): if I1={} then Expa:=Expa union {permu}: B[permu]:=Refinements(permu,Patterns,J): C[permu]:=J: KVU:=KVU minus {permu}: KVU:=KVU union (B[permu] minus (Expa union Redu)): fi: if I1<>{} then Redu:=Redu union {permu}: i:=op(1,I1): A[permu]:=i: C[permu]:=J: permud:=[op(1..i-1,permu),op(i+1..nops(permu),permu)]: permud:=redu(permud): KVU:=KVU minus {permu}: if not member(permud,Redu union Expa) then KVU:=KVU union {permud}: fi: fi: od: [Redu,Expa,A,B,C]: end: #CUBE(DIM,SIZE): All vectors of length DIM, with non- #negative components <=SIZE CUBE:=proc(DIM,SIZE) local gu,i,mu,akha,lu: if DIM<0 or SIZE<0 then ERROR(`Bad input`): fi: if DIM=0 then RETURN({[]}): fi: lu:={}: gu:=CUBE(DIM-1,SIZE): for i from 1 to nops(gu) do mu:=op(i,gu): for akha from 0 to SIZE do lu:=lu union {[op(mu),akha]}: od: od: lu: end: #Given a vector v=[v_1,v_2,...,v_n] satisfying #v_1-v_2>=d, v_2-v_3>=d, ..., v_{d-1}-v_d>=d #finds the next vector in lex. order nextvect:=proc(v,d) local i,j,n: n:=nops(v): for j from n by -1 to 2 while op(j-1,v)-op(j,v)=d do od: [op(1..j-1,v),op(j,v)+1,seq(d*(n-j-i+1)+1,i=1..n-j)]: end: #applymonoJul(mono,Sch,tau,vec,A): applies a single monomial in the #shift operators A[1], ..., A[n] to f at vec applymonoJul:=proc(mono,Sch,tau,vec,A) local i,vec1,mu,vec2: if type(mu,`+`) then ERROR(`must be a product`): fi: mu:=mono: vec1:=[]: for i from 1 to nops(vec) do vec1:=[op(vec1),op(i,vec)-degree(mono,A[i])]: mu:=normal(mu/A[i]^degree(mono,A[i])): od: vec2:=[seq(vec1[i]+nops(vec1)-i+1,i=2..nops(vec1))]: mu*JulianSort(vec1[1]+nops(vec1)-1,tau,vec2,Sch): end: #applyopeJul(ope,Sch,tau,vec,A) Given an operator ope, expressed in #terms of A[i], finds ope JulianSort(vec) applyopeJul:=proc(ope,Sch,tau,vec,A) local i, gu: if not type(ope,`+`) then RETURN(applymonoJul(ope,Sch,tau,vec,A)): fi: gu:=0: for i from 1 to nops(ope) do gu:=gu+applymonoJul(op(i,ope),Sch,tau,vec,A): od: expand(gu): end: #Sim2(k,L): All vectors of integers with k components #[v_1,...,v_k] such that 1<=v_1{} and hakhi<=GVUL do permu:=op(1,KVU): if nops(permu)>hakhi then hakhi:=nops(permu): fi: if hakhi>GVUL then RETURN(0,[Redu,Expa,A,B,C]): fi: J:=Findj(permu,Patterns): I1:=FindiJ(permu,Patterns,J): if I1={} then Expa:=Expa union {permu}: B[permu]:=Refinements(permu,Patterns,J): C[permu]:=J: KVU:=KVU minus {permu}: KVU:=KVU union (B[permu] minus (Expa union Redu)): fi: if I1<>{} then Redu:=Redu union {permu}: i:=op(1,I1): A[permu]:=i: C[permu]:=J: permud:=[op(1..i-1,permu),op(i+1..nops(permu),permu)]: permud:=redu(permud): KVU:=KVU minus {permu}: if not member(permud,Redu union Expa) then KVU:=KVU union {permud}: fi: fi: od: [Redu,Expa,A,B,C]: end: #SchemeU(Patterns,GVUL): Given a set of permutations (of various lengths) #Patterns, finds an Upper scheme [Redu,Expa,A,B,C] where #Redu is a set of reducible prefix-permutations with respect to #Patterns, Expa is the set of expandible ones, A is a table #where A[sigma] is the place (between 1 and nops(sigma)) that says #where to delete, for sigma in Redu, #and C[sigma] is a table that says the Findj of sigma #[i.e. the restriction, if any on the prefixes] #and B is a table #where B[sigma] is the set of refinements of sigma, of length one #more. For this to be a scheme, the deleted versions of all #sigma in Redu, and all the B[sigma] must be "good" #GVUL is the maximal permitted prefix size on either #Redu or Expa, #The Upper Scheme means that the sequence computed by Miklos #using the scheme is #an UPPER bound for the actual sequence of cardinalities #of Wilf(n,Patterns) SchemeU:=proc(Patterns,GVUL) local Redu,Expa,A,B,C,i,KVU,permu,permud,I1,J: option remember: if member([],Patterns) then ERROR(`[] can't be a member of Patterns`): fi: KVU:={[]}: Redu:={}: Expa:={}: while KVU<>{} do permu:=op(1,KVU): J:=Findj(permu,Patterns): I1:=FindiJ(permu,Patterns,J): if I1={} then if nops(permu){} then Redu:=Redu union {permu}: i:=op(1,I1): A[permu]:=i: C[permu]:=J: permud:=[op(1..i-1,permu),op(i+1..nops(permu),permu)]: permud:=redu(permud): KVU:=KVU minus {permu}: if not member(permud,Redu union Expa) then KVU:=KVU union {permud}: fi: fi: od: [Redu,Expa,A,B,C]: end: #SchemeSeries(schem,L): Given a scheme schem, and an integer #L, computes the first L terms of the #sequence generated by schem SchemeSeries:=proc(schem,L) local gu,i: gu:=[]: for i from 0 to L do gu:=[op(gu),Miklos(i,schem)]: #print(gu); od: gu: end: Zinn:=proc(resh,n) local s1,s2: s1:=sn(resh,n): s2:=sn(resh,n-1): evalf(2*(s1+s2)/(s1-s2)^2), evalf(sqrt(op(n+1,resh)/op(n-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))): end: sn:=proc(resh,n): -1/log(op(n+1,resh)*op(n-1,resh)/op(n,resh)^2): #evalf(%): end: #checkrec(ope,n,N,Seq1) applies the operator ope(N,n) to the sequence gu checkrec:=proc(ope,n,N,Seq1) local i,j,lu,ru: lu:=[]: for i from 1 to nops(Seq1)-degree(ope,N) do ru:=0: for j from 0 to degree(ope,N) do ru:=ru+subs(n=i,coeff(ope,N,j))*op(i+j,Seq1): od: lu:=[op(lu),ru]: od: lu: end: findrec:=proc(f,DEGREE,ORDER,n,N) local ope,var,eq,a,i,j,n0,kv,var1,eq1,mu: if (1+DEGREE)*(1+ORDER)+2+ORDER>nops(f) then ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE): fi: ope:=0: var:={}: for i from 0 to ORDER do for j from 0 to DEGREE do ope:=ope+a[i,j]*n^j*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 do eq1:=0: for i from 0 to ORDER do eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f): od: eq:= eq union {eq1}: od: var1:=solve(eq,var): kv:={}: for i from 1 to nops(var1) do mu:=op(i,var1): if op(1,mu)=op(2,mu) then kv:= kv union {op(1,mu)}: fi: od: ope:=subs(var1,ope): if nops(kv)>1 then RETURN(0): fi: for i from 1 to nops(kv) do ope:=subs(op(i,kv)=1,ope): od: ope: end: Findrec:=proc(f,n,N) local i1,kal,DEGREE, ORDER,L,gu,gu1,i: for L from 1 to nops(f)-5 do for ORDER from 1 to min(L,trunc((nops(f)-3)/2)-1) do for DEGREE from 0 to min(L-ORDER, trunc((nops(f)-3)/2/(1+ORDER))-2) do gu:=findrec(f,DEGREE,ORDER,n,N): if gu<>0 then gu1:=0: for i from 0 to degree(gu,N) do gu1:=gu1+factor(coeff(gu,N,i))*N^i: od: kal:=checkrec(gu1,n,N,f): if kal=[seq(0,i1=1..nops(kal))] then RETURN(gu1): fi: fi: od: od: od: 0: end: #SchemeRecurrence(schem,n,N,L): Given a scheme schem, finds #empirically a linear recurrence equation with #polynomial coefficients (where n is the discrete argument #and N is the shift operator in n) satisfied by it, by trying #up to L terms. If unsuccesful, it returns the sequence SchemeRecurrence:=proc(schem,n,N,L) local gu,i,mu: gu:=[]: for i from 0 to L do gu:=[op(gu),Miklos(i,schem)]: if i>10 and i mod 5 =0 then mu:=Findrec(gu,n,N): if mu<>0 then RETURN(mu): fi: fi: od: gu: end: #PermSets1(Resh,a) : Given a list of integers Resh, and #an integer a, finds the set of all permutations #that contail Resh[i] permutations of length a+i-1 #For example, try: PermSets1([2,1],3); PermSets1:=proc(Resh,a) local gu,Resh1,mu,ku,j,i: if Resh=[] then RETURN({{}}): fi: if Resh[nops(Resh)]=1 then Resh1:=[op(1..nops(Resh)-1,Resh)]: else Resh1:=[op(1..nops(Resh)-1,Resh),Resh[nops(Resh)]-1]: fi: mu:=PermSets1(Resh1,a): ku:=permute(nops(Resh)+a-1): gu:={}: for i from 1 to nops(mu) do for j from 1 to nops(ku) do if SubPatts(ku[j]) intersect mu[i]={} then gu:=gu union {mu[i] union {ku[j]}}: fi: od: od: gu: end: #PermSets(Resh,a): Like PermSets1(Resh,a) but only with #one representative for each equivalence class under #the dihedral group PermSets:=proc(Resh,a) local gu,lu: lu:=PermSets1(Resh,a): gu:={}: while lu<>{} do gu:=gu union {lu[1]}: lu:=lu minus KHAVERIM(lu[1]): od: gu: end: SubPatts:=proc(pi) local gu,n,i: n:=nops(pi): gu:=powerset({seq(i,i=1..n)}) minus {{},seq({i},i=1..n)}: {seq(redu([seq(pi[gu[i][j]],j=1..nops(gu[i]))]),i=1..nops(gu))}: end: #FindSchemes(Resh,a,Depth,L): Given a list of integers Resh, and an #integers a, Depth and L, finds the set of all permutation #sets with Resh[i] permutations of length a+i-1 (i=1..nops(Resh)) #that have a scheme of depth<=Depth, followed by their scheme, and #the first L terms. For example, try: #FindSchemes([1],3,2,10); FindSchemes:=proc(Resh,a,Depth,L) local gu,lu,mu,ku,i: gu:={}: lu:=PermSets(Resh,a): print(`There are `, nops(lu), `equivalence classes `): for i from 1 to nops(lu) do mu:=lu[i]: ku:=SchemeF(mu,Depth): if ku<>0 then ku:=[op(ku),SchemeSeries(ku[3],L)]: gu:=gu union {ku}: fi: od: print(`There are`, nops(gu), `successul ones with depth `, Depth): gu: end: