print("Welcome to mwords, a Maple program by Lara Pudwell"):
print("for finding the multivariate generating function to count pattern-avoiding words"):
print(""):
print("This version was last updated on 19 March 2006"):
print("Try Help(); for a list of available functions"):
print("or Help(function_name); for more information about a particular function"):

Help:=proc()
if args=NULL then
print(`red(Li), pats(k,Li2), genfun(k,Li,x), gwords(n,k,Li), bfwords(n,k,Li), test(n,k,Li)`):

elif nargs=1 and args[1]=red then
print(`red(Li): inputs a multiset permutation, and returns the corresponding reduced multiset permutation`):
print(`for example: red([1,4,5,4]) returns [1,2,3,2]`):

elif nargs=1 and args[1]=pats then
print(`pats(k,Li2): inputs an integer k and a set of multiset permutations Li2`):
print(`and returns all ways of writing patterns in Li2 from an alphabet of k letters`):
print(`for example: pats(4,{[1,2,3],[1,3,2]}) returns {[1, 2, 3], [1, 2, 4], [1, 3, 4], [2, 3, 4], [1, 3, 2], [2, 4, 3], [1, 4, 3], [1, 4, 2]}`):

elif nargs=1 and args[1]=genfun then
print(`genfun(k,Li,x): returns the multivariable generating function in x[1],...,x[k] for the words on k letters that avoid all permutations in the set Li`):
print(`so that the coefficient of mul(x[i]^a[i],i=1..k) is the number of words on a[1] 1s, a[2] 2s,..., a[k] ks avoiding all permutations in Li`):
print(`for example: genfun(3,{[1,2,3,4]},x) returns 1/(1-x[1]-x[2]-x[3]) since all words on 3 letters avoid the pattern [1,2,3,4]`):

elif nargs=1 and args[1]=gwords then
print(`gwords(n,k,Li): returns all words of length n on k letters that avoid all patterns in the set Li`):
print(`for example: gwords(3,2,{[1,2]}); returns {[1, 1, 1], [2, 2, 2], [2, 2, 1], [2, 1, 1]}`):

elif nargs=1 and args[1]=bfwords then
print(`bfwords(n,k,Li): finds the list of all words of length n on k letters that avoid all patterns in Li`):
print(`then converts them into products of variables`):
print(`for example: gwords(3,2,{[1,2]}) returns {[1, 1, 1], [2, 2, 2], [2, 2, 1], [2, 1, 1]}`):
print(`so bfwords(n,k,Li): returns x[1]^3 + x[2]^3 + x[1] *x[2]^2 + x[1]^2 *x[2]`):

elif nargs=1 and args[1]=test then
print(`test(n,k,Li): compares the result of genfun(n,k,Li) and of add(bfwords(j,k,Li),j=0..n)`):
print(`returns true if they are equal, and false otherwise`):

else print(args[1], " is not a valid function name"):
fi:
end:

with(combinat):

test:=proc(n,k,Li)
evalb(mtaylor(genfun(k,Li,x),{seq(x[j],j=1..k)},n+1)-add(bfwords(j,k,Li),j=0..n)=0);
end:

genfun:=proc(k,Li,x) local i, m:
option remember:
#take care of case where there is a singleton pattern
for i from 1 to nops(Li) do
 if nops(Li[i])=1 then return 1: fi:
od:

#take care of case when there are no patterns
if nops(Li)=0 then
return 1/(1-add(x[i],i=1..k)): 
fi:

#m is the smallest number of distinct letters in a single pattern
m:=nops({op(Li[1])}):
for i from 1 to nops(Li) do
if nops({op(Li[i])})<m then m:=nops({op(Li[i])}): fi:
od:

#if there are less letters in the alphabet than in *any* pattern
#return generating function for ALL words on k letters
if k<m then return 1/(1-add(x[i],i=1..k)):
fi:

#if in a non-trivial case, return genfun for alphabet {1,...,k}
#and the set of patterns from before
func2({seq(i, i=1..k)},pats(k, Li),x):

end:

#func2 inputs "literal" patterns, not reduced and works recursively
func2:=proc(alphabe, Li,x) local i, tot, j, temppat, m, temppat2, finis, singles, newalphabe, newLi, newLi2, otherpat, sumalph:
option remember:
tot:=0:

otherpat:={}:
for j from 1 to nops(alphabe) do

temppat:={}:
 for m from 1 to nops(Li) do
  if Li[m][1]=alphabe[j] then temppat:=temppat union{Li[m]}: else 
otherpat:=otherpat union {Li[m]}: fi:
 od:

temppat2:={}:
for m from 1 to nops(temppat) do:
 temppat2:=temppat2 union {[op(2..nops(temppat[m]), temppat[m])]}:
od:

#temppat is the set of patterns beginning with j
#temppat2 is the set of truncated patterns beginning with j

singles:={}:
for i from 1 to nops(temppat2) do
 if nops(temppat2[i])=1 then singles:=singles union {op(temppat2[i])}: fi:
od:
#singles is the list of singleton patterns

temppat2:=temppat2 minus {[]}:
for i from 1 to nops(singles) do
 temppat2:=temppat2 minus {[singles[i]]}:
od:

newalphabe:=alphabe minus singles:
newLi:=temppat2 union Li:
newLi2:={}:

#newalphabe is the letters allowed to use without creating a forbidden pattern
#newLi is the original set of literal patterns along with the truncated ones

for i from 1 to nops(newLi) do
 if {op(newLi[i])} subset newalphabe then
  newLi2:={op(newLi2), newLi[i]}:
 fi:
od:

#we just removed all patterns not consisting of letters in the current alphabet

newLi:=newLi2:

#temppat2 is a set of chopped permutations beginning with letter j
# so the set of all words avoiding Li starting with j must avoid temppat2
#tot:=tot+func2(newalphabe, newLi):
#this should sum over all j to give a generating function

sumalph:=0:
for i from 1 to nops(newalphabe) do
 sumalph:=sumalph+x[newalphabe[i]]:
od:
#sumalph is the sum of all x[i] for i in newalphabe

if nops(newalphabe)=0 then
	tot:=tot+x[alphabe[j]]:
elif nops(newLi)=0 then
	tot:=tot+ (x[alphabe[j]]/(1-sumalph)):
elif newalphabe = alphabe and newLi = Li then
	tot:=x[alphabe[j]]*y+tot:
else
	tot:=tot+x[alphabe[j]]*func2(newalphabe, newLi,x):
fi:
od:

tot:=solve(y=tot+1,y):
tot:
end:

#(red, reduces Li)
red:=proc(Li) local i, temp, word, temp2, k, curr, rests:
temp:=sort(Li): temp2:=[1]: k:=1:
for i from 2 to nops(Li) do
 if temp[i]>temp[i-1] then k:=k+1: fi:
temp2:=[op(temp2), k]:
od:
word:=sort(Li):
rests:={}:
for i from 1 to nops(Li) do
 rests:=rests union {word[i]=temp2[i]}:
od:
word:=subs(rests, Li):
end:

pats:=proc(k,Li2) local 
word,i,j,newalpha,newpat,newword,letters,m,allpats,Li:
Li:={}:
for i from 1 to nops(Li2) do
 Li:={op(Li), red(Li2[i])}:
od:

allpats:={}:
for m from 1 to nops(Li) do

	word:=Li[m]:
	for i from 1 to k do
	word:=subs(i=x[i], word):
	od:
	#word is now the pattern in terms of x_i's

	newalpha:=choose(k,nops({op(word)})):
	#(this gives us the set of letters to work with in our new pattern set)

	newpat:={}:
	for i from 1 to nops(newalpha) do
	newword:=word:

	letters:=[seq(x[j]=newalpha[i][j], j =1..nops(newalpha[i]))]: 
	newword:=subs(letters,newword):

	newpat:={op(newpat), newword}:
	od:
	
	#newpat is the set of all possible permutations that are equivalent to Li[k]
allpats:=allpats union newpat:
od:

allpats:
end:

words:=proc(n, k) local alphaLi, i:
alphaLi:=[]:
for i from 1 to k do
alphaLi:=[op(alphaLi), i$n]:
od:
alphaLi:=permute(alphaLi, n):
#alphaLi is the list of possible words of length n with k letters.
end:

gwords:=proc(n,k,Li) local temp, poss, badw,i,j,m, tempword, s,tempw,t, 
Li2:
Li2:={}:
for i from 1 to nops(Li) do
 Li2:=Li2 union{red(Li[i])}:
od:

temp:=words(n,k):
badw:={}:

for i from 1 to nops(Li2) do
 for j from 1 to nops(temp) do
  poss:=choose({seq(m, m=1..nops(temp[j]))}, nops(Li2[i])): 
#poss is the set of all positions in temp of length of pattern i

for t from 1 to nops(poss) do
tempw:=temp[j]:
 tempword:=[]:
 for s from 1 to nops(poss[t]) do
  tempword:=[op(tempword), tempw[poss[t][s]]]:
 od:

 if red(tempword)=Li2[i] then badw:=badw union {temp[j]}: fi:
od:
 od:
od:

{op(temp)} minus badw:

end:

bfwords:=proc(n,k,Li) local temp, i, var, Li2,j,tot:

Li2:=gwords(n,k,Li):
tot:=0:

for j from 1 to nops(Li2) do

temp:=convert(Li2[j],`multiset`):
var:=1:
for i from 1 to nops(temp) do
 var:=var*x[temp[i][1]]^temp[i][2]:
od:
tot:=tot+var:

od:

tot:
end: