##SymUGJ: Save this file as SymUGJ. To use it, stay in the           #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read SymUGJ : <Enter>                              #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Temple University ,                   #
#zeilberg@math.temple.edu.                                           # 
#######################################################################

#Created: April 2001
#This version: April 2001
#SymUGJ: A Maple package that finds Umbral Schemes
#for enumerating words that avoid as factors
#an infinite family of parametrized "mistakes"
#where the set of mistakes is symmetric under the
#symmetric group on the letters
#
#It is one of the packages that accompany Doron Zeilberger's article:
#The Umbral Transfer Matrix Method. V. The Goulden-Jackson
#Cluster Method for Infinitely Many Mistakes
#Please report bugs to zeilberg@math.temple.edu


print(`SymUGJ: A Maple package that finds Umbral Schemes`):
print(`for enumerating words that avoid as factors`):
print(` an infinite family of parametrized "mistakes"`):
print(` where the set of mistakes is symmetric w.r.t. permuting letters`):
print(`It is one of the packages that accompany Doron Zeilberger's article:`):
print(`The Umbral Transfer Matrix Method. V. The Goulden-Jackson`):
print(`Cluster Method for Infinitely Many Mistakes`):
print(`Please report bugs to zeilberg@math.temple.edu`):

lprint(``):
print(`Created: May 2001`):
print(`This version: May 2001`):
lprint(``):
print(`Written by Doron Zeilberger, zeilberg@math.temple.edu`):
lprint(``):
print(`Please report bugs to zeilberg@math.temple.edu`):
lprint(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.temple.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()
print(``):
if nops([args])=0  then
print(` The Main Procedures are : `):
print(`  CheckU, UGJseries, WtUmSc, YafeUmSc `): 
 print(`For help with`):
 print(`a specific procedure, type ezra(procedure_name)`):

fi:

if nops([args])=1 and args[1]=CheckU  then
print(`CheckU(SymbMisSet,Alphabet,L): Uses the Classical Goulden-`):
print(`Jackson to find the first L terms of the set of words in`):
print(`the alphabet Alphabet that avoid factors from SymbMisSet`):
print(`where each symbolic mistake has the form [[letter1,expression1], ...]`):
print(`and expression1... are affine linear expressions in the variables`):
print(`of DisVars`):
fi:

if nops([args])=1 and args[1]=FindUmb  then
print(`FindUmb(U,V,x,t,d): Given two symbolic words, U, and V,`):
print(`and a variable x (x an index variable, for catalytic`):
print(`variables), and a variable t, to keep track of length`):
print(`computes the umbral edge connecting U and V`):
print(` also taking account of the images of V under permutations`):
print(` on d letters`):
print(`For example try: FindUmb([[[1,1],[2,a],[3,a],[1,1]],[a]],`):
print(`[[[1,1],[2,a],[3,a],[1,1]],[a]],x,t,d):`):
fi:

if nops([args])=1 and args[1]=FindUmbInt  then
print(`FindUmbInt(U,V,i,x,t): Given two symbolic letters U, V`):
print(`expressed in symbolic frequency notation`):
print(`{[[let1,AffLinExp1],[let2,AffLinExp2], ...], {set_of_symbols}}`):
print(`and a location i, one of the legal interfaces of`):
print(`U before V in a symbolic cluster, and a sumbol x, `):
print(`for the catalytic variabes and a symbol t, find the umbra that`):
print(`sends x_1^a_1*x_2^a_2 ... to the generating function of all possible`):
print(`such V's `):
print(`For example try: FindUmbInt([[[1,1],[2,a],[3,a],[1,1]],[a]],`):
print(`[[[1,1],[2,a],[3,a],[1,1]],[a]],4,x,t):`):
fi:

if nops([args])=1 and args[1]=UGJseries  then
print(`UGJseries(SetSW,NumLetters,L): Given a set of symbolic`):
print(`letters, SetSW, and the number of letters, NumLetters,`):
print(`an integer, and an integer L, uses the Umbral Scheme`):
print(`to find the first L terms of the series expansion`):
print(`of words that avoid the mistakes of SetSW`):
print(`and all their images  under the symmetric group`):
print(`For example, type`):
print(`UGJseries([ [[[1,a+1],[2,a+b+1],[1,b+1]],[a,b]]],2,20):`):
fi:

if nops([args])=1 and args[1]=WtUmSc   then
print(`WtUmSc(ListSW,x,t,d): Given a list of symbolic words`):
print(`ListSW, constructs a Weighted Umbral Scheme for the Cluster`):
print(`generating function `):
print(`where the mistakes are the members of ListSW and all their images`):
print(` under the symmetric group on d letters`):
print(`For example try: `):
print(`WtUmSc( [ [[[1,a+1],[2,a+b+1],[1,b+1]],[a,b]]],x,t,2):`):
fi:

if nops([args])=1 and args[1]=YafeUmSc then
print(`YafeUmSc (UmbSc,F,D): Given a Weighted Umbral `):
print(`Scheme UmbSc, in`):
print(`Maple notation, converts it to human notation`):
print(`using the symbols F[1], F[2] for the functions and D for `):
print(`differntiation of respective variables`):
fi:

end:

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

if not type(poly2,`+`) then
RETURN(SimplifyUmbra({UmbralTerm(poly2,xlist,alist)})):
fi:

gu:={}:

for i from 1 to nops(poly2) do
 gu:=gu union {UmbralTerm(op(i,poly2),xlist,alist)}:
od:

SimplifyUmbra(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,xlist,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,xlist,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,xlist,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,mono1):
 mu2:=op(2,mono1):
  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:


#SimplifyUmbra(Umb): Given an Umbra, collects like terms
#and returns a simplified version
SimplifyUmbra:=proc(Umb)
local Umb1,gu,kv,i,mu,evar,F,pu:
if Umb=0 then
ERROR(FAIL):
fi:

Umb1:=ConvertToUmbra2(Umb):

gu:=0:
kv:={}:

for i from 1 to nops(Umb1) do
 gu:=gu+Umb1[i][1]*F(Umb1[i][2]):
 kv:=kv union {Umb1[i][2]}:
od:

gu:=expand(gu):

mu:={}:

for i from 1 to nops(kv) do
evar:=op(i,kv):
pu:=normal(coeff(gu,F(evar))):
if pu<>0 then
mu:=mu union {[factor(pu),evar[1],evar[2]]}:
fi:

od:

mu:

end:

#ConvertToUmbra2(Umb): Given an umbra Umb in the 3-list-format
#i.e. as a set of 3-lists of the form
#[FRONT,Orders_Of_Derivaties,Substitutions]
#converts this to a set of 2-lists of the form 
#[FRONT,[Orders_Of_Derivaties,Substitutions]]
ConvertToUmbra2:=proc(Umb)
local Umb1,i,mu:
Umb1:={}:

for i from 1 to nops(Umb) do
 mu:=op(i,Umb):
 Umb1:=Umb1 union {[mu[1],[mu[2],mu[3]]]}:
od:

Umb1:
end:






# ApplyUmbralMatrix(UmbraM,Fvec,ylistvec): Given an Umbral matrix
# UmbraM (in terms of a list of lists), a vector of expressions,
# Fvec, and a vector to variable-lists corresponding to the
# components of UmbraM and Fvec
#
ApplyUmbralMatrix:=proc(UmbraM,Fvec,ylistvec)
local gu,mu,i,j,k:

k:=nops(UmbraM):

if k<>nops(Fvec) or k<>nops(ylistvec) then
ERROR(`Bad input`):
fi:

gu:=[]:

for i from 1 to k do
mu:=0:
for j from 1 to k do
mu:=normal(mu+ApplyUmbra(UmbraM[i][j],Fvec[j],ylistvec[j])):
od:
gu:=[op(gu),mu]:
od:
gu:
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:
gu:=0:

for i from 1 to nops(Umbra1) do
 gu:=normal(gu+ApplyUmbraTerm(Umbra1[i],f,ylist)):
od:
gu:
end:


# ApplyUmbraTerm(Umbra1,f,ylist): given an umbral term, Umbra1, applies
# it to the expression f in the variables in ylist
ApplyUmbraTerm:=proc(Umbra1,f,ylist)
local i, FRONT,Diffs,SubsList,gu,bu:

FRONT:=Umbra1[1]:
Diffs:=Umbra1[2]:
SubsList:=Umbra1[3]:

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


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

gu:=f:

for i from 1 to nops(ylist) do
gu:=normal(xdiff1(gu,ylist[i],Diffs[i])):
od:

bu:={}:
for i from 1 to nops(ylist) do
bu:=bu union {ylist[i]=SubsList[i]}:
od:
gu:=subs(bu,gu):

normal(FRONT*gu):

end:


# xdiff1(f,x,i): given an expression f, and a variable x,
# and an integer i, finds (xD_x)^i f
# 
xdiff1:=proc(f,x,i)
local gu,j:

if i=0 then
 RETURN(f):
fi:

gu:=f:

for j from 1 to i do
gu:=x*diff(gu,x):
od:

gu:

end:




#SerToPol(sidra,q,n): given a series, sidra, finds the polynomial
#consisting of the terms up to degree n
SerToPol:=proc(sidra,q,n)
local lu,i,gu:
lu:=taylor(sidra,q=0,n+1):
gu:=0:
for i from 0 to n do
gu:=gu+coeff(lu,q,i)*q^i:
od:
gu:
end:




# ApplyUmSc(UmSch,q,n,var): Given an Umbral Scheme, UmSch, finds
# the generating functions of creatures up to the wt. q^n
#followed by substition of all x and y variables to 1 and summing
#over leftmost letters
ApplyUmSc:=proc(UmSch,q,n,var)
local UM,INI,ylistvec,lu,gu,i,lu1,i1,kv,ku,fu,Sid,j,mu,gc:
if not type(n,integer) or not n>0 then
ERROR(`n>0`):
fi:


kv:=UmSch[1]:
UM:=UmSch[2]:
INI:=UmSch[3]:
ylistvec:=UmSch[4]:

lu:=[]:

for i from 1 to nops(INI) do
lu:=[op(lu),SerToPol(INI[i],q,0)]:
od:

if n=0 then
RETURN(lu):
fi:

Sid:=[]:

for i from 1 to n do

 lu1:=ApplyUmbralMatrix(UM,lu,ylistvec):
 gu:=[]:

for i1 from 1 to nops(lu1) do
 gu:=[op(gu),expand(SerToPol(INI[i1],q,i)+SerToPol(normal(lu1[i1]),q,i))]:
od:
lu:=gu:

fu:=[]:
for i1 from 1 to nops(lu) do
fu:=[op(fu),coeff(expand(lu[i1]),q,i)]:
od:

mu:=0:

for j from 1 to nops(kv) do
 mu:=mu+kv[j]*fu[j]:
od:


ku:={seq(var[gc]=1,gc=1..nops(var))}:
mu:=subs(ku,normal(mu)):

Sid:=[op(Sid),normal(mu)]:
lprint(Sid):

od:

lu,Sid:

end:
# ApplyUmSc1(UmSch,q,n,var): Same as
#ApplyUmSc(UmSch,q,n,var), but no partial results
#displayed
ApplyUmSc1:=proc(UmSch,q,n,var)
local UM,INI,ylistvec,lu,gu,i,lu1,i1,kv,ku,fu,Sid,j,mu,gc:
if not type(n,integer) or not n>0 then
ERROR(`n>0`):
fi:


kv:=UmSch[1]:
UM:=UmSch[2]:
INI:=UmSch[3]:
ylistvec:=UmSch[4]:

lu:=[]:

for i from 1 to nops(INI) do
lu:=[op(lu),SerToPol(INI[i],q,0)]:
od:

if n=0 then
RETURN(lu):
fi:

Sid:=[]:

for i from 1 to n do

 lu1:=ApplyUmbralMatrix(UM,lu,ylistvec):
 gu:=[]:

for i1 from 1 to nops(lu1) do
 gu:=[op(gu),expand(SerToPol(INI[i1],q,i)+SerToPol(normal(lu1[i1]),q,i))]:
od:
lu:=gu:

fu:=[]:
for i1 from 1 to nops(lu) do
fu:=[op(fu),coeff(expand(lu[i1]),q,i)]:
od:

mu:=0:

for j from 1 to nops(kv) do
 mu:=mu+kv[j]*fu[j]:
od:


ku:={seq(var[gc]=1,gc=1..nops(var))}:
mu:=subs(ku,normal(mu)):

Sid:=[op(Sid),normal(mu)]:

od:

lu,Sid:

end:

#YafeUm(Umb1,F): Given an umbral operator, Umb1, in
#Maple notation, converts it to human notation
#using the symbol F for the function and D for differntiation
YafeUm:=proc(Umb1,F,D)
local T,gu,i,mu,lu,gu1,lu1:

gu:={}:

for i from 1 to nops(Umb1) do
mu:=op(i,Umb1):
gu:=gu union {[mu[2],mu[3]]}:
od:

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


for i from 1 to nops(Umb1) do
mu:=op(i,Umb1):
T[[mu[2],mu[3]]]:=normal(T[ [mu[2],mu[3]] ]+mu[1]):
od:

lu:=0:

for i from 1 to nops(gu) do
gu1:=gu[i]:
lu1:=F(op(gu1[2])):

if convert(gu1[1],`+`)<>0 then
 lu1:=lu1*D[op(gu1[1])]:
fi:
lu1:=lu1*factor(T[gu1]):
lu:=lu+lu1:
od:

lu:

end:


#YafeUmSc(UmbSc,F,D): Given a Weighted umbral 
#Scheme UmbSc, in
#Maple notation, converts it to human notation
#using the symbols F[1], F[2] for the functions and D for differntiation
#of respective variables
YafeUmSc:=proc(UmbSc,F,D)
local gu1,gu2,gu3,gu4,eq,eq1,n,i,j:
gu1:=UmbSc[1]:
gu2:=UmbSc[2]:
gu3:=UmbSc[3]:
gu4:=UmbSc[4]:
eq:=[]:
n:=nops(gu2):
for i from 1 to n do
eq1:=gu3[i]:

for j from 1 to n do
 eq1:=eq1+YafeUm(gu2[i][j],F[j],D):
od:
eq:=[op(eq),F[i](op(gu4[i]))=eq1]:
od:
eq,gu1:
end:
#############End Programs Taken from ROTA###############

#############PROGRAMS TAKEN FROM LinDiophantus
#GFxt(a,E1,x,t): Given a list of symbols, a
#denoting generic (symbolic) non-negative integers
#and a list, E1, of affine linear expressions with integer
#coeffs. in them , and symbols x and t,
#finds the generating function:
#sum of  x[1]^a[1]*...*x[n]^a[n]*t[1]^E1[1]*...*t[r]^E1[r]
#(where n=nops(a) and r=nops(E1)), and the sum is over 0<=a[1]<infinity, ...,
#0<=a[n]<infinity, ...,
GFxt:=proc(a,E1,x,t)
local i,n,gu,r,j:
gu:=1:
n:=nops(a):
r:=nops(E1):

for i from 1 to n do
gu:=gu*x[i]^a[i]:
od:

for i from 1 to r do
gu:=gu*t[i]^E1[i]:
od:

for i from 1 to n do
gu:=subs(a[i]=j,gu):
gu:=sum(gu,j=0..infinity):
gu:=normal(gu):
od:
gu:
end:

#WhoAreYout(f,t): Given a Laurent polynomial f in t
#decides it is of the form P(t), where P is a polynomial in t
#or P(1/t)/t and returns 1 and -1 respectively
WhoAreYout:=proc(f,t)
local f1:
f1:=expand(f):
if degree(f1,t) >0 and ldegree(f1,t)>=0 then 1:
elif degree(f1,t) <0 then -1:
else 0:
fi:
end:

#WhoAreYouxt(f,x,t): Given a polynomial f, in the variables
#x and t, looks at the `normalized' polynomial (in x)
#obtained by dividing by the coeff. of x^0
#and seeing whether all the  coeffs. of the powers of x
#(that are Laurent polynomials in t) have consistently
#coeffs that are polynomials in t, or consistently coeffs.
#that are polynomials in 1/t, or none of the above
#and returns 1, -1, 0 respectively

WhoAreYouxt:=proc(f,x,t)
local f1,lu,i1,i:
f1:=expand(f):

if degree(f1,t)=0 then
RETURN(1):
fi:

if degree(f1,t)<>0 and normal(f1-subs(t=1,f1)*t^degree(f1,t))=0 then
RETURN(-1):
fi:


if coeff(f1,x,0)<>0 then
f1:=expand(f1/coeff(f1,x,0)):
else
f1:=f:
fi:

for i1 from 1 while expand(coeff(f1,x,i1))=0 do od:
lu:=WhoAreYout(coeff(f1,x,i1),t):

for i from i1 to degree(f1,x) do
if WhoAreYout(coeff(f1,x,i),t)<>lu then
 RETURN(0):
fi:
od:
lu:
end:


#isbad2(evar,xSet): given an expression in the variables
#in the set xSet, decides whether evar has positive degree
#in any of the variables of xSet

isbad2:=proc(evar,xSet)
local i:

for i from 1 to nops(xSet) do
if degree(evar,xSet[i])>0 then
RETURN(true):
fi:
od:
false:
end:

#CT1old(f,xSet,t): Given a rational
#function f of t and a set, xSet, of  x-variables
#finds the constant term in t (i.e. the coeff. of t^0)
#of the Laurent expansion of t, viewed as a formal Laurent
#series in the x's of xSet
#For example try CT1old(1/(1-x*t)/(y-t),[x,y],t);
CT1old:=proc(f,xSet,t)
local f1,i,lu,g,mu,BAD1:
mu:=denom(normal(f)):
mu:={solve(mu,t)}:

BAD1:={}:

for i from 1 to nops(mu) do
if isbad2(mu[i],xSet) or mu[i]=0 then
 BAD1:=BAD1 union {mu[i]}:
fi:
od:

f1:=convert(f,parfrac,t):

g:=f1:
for i from 1 to nops(f1) do
lu:=op(i,f1):
mu:=denom(lu):
if isbad3(mu,BAD1,t) then
g:=g-lu:
fi:

od:


factor(normal(subs(t=0,g))):
end:
 


#isbad3(evar,kv,t): if any of the substitutions of t=kv[i]
#into evar result in 0
isbad3:=proc(evar,kv,t)
local i:

for i from 1 to nops(kv) do

 if expand(subs(t=kv[i],evar))=0 then
  RETURN(true):
 fi:

od:

false:
end:



isbad:=proc(evar,t,xSet)
local i,j,coe:

for i from 0 to degree(evar,t)-1 do
 coe:=coeff(evar,t,i):
  for j from 1 to nops(xSet) do
    if degree(coe,xSet[j])>0 then
     RETURN(true):
    fi:
  od:
od:
false:
end:

#isbada(evar,BAD1,t): Given a denominator term (coming
#from a single term of a prtial-fraction decomposition w.r.t t)
#decides whether it is of the form BAD1[i]^(some power) for
#some bad term of BAD1
isbada:=proc(evar,BAD1,t)
local j,gu,d1:
if normal(diff(evar,t))=0 then
 RETURN(false):
fi:
for j from 1 to nops(BAD1) do
gu:=BAD1[j]:
d1:=degree(evar,t)/degree(gu,t):
if type(d1,integer) then
  gu:=normal(evar/gu^d1):
  if normal(diff(gu,t))=0 then
   RETURN(true):
 fi:
fi:
od:
false:
end:

#CT1(f,xSet,t): Given a rational
#function f of t and a set, xSet, of  x-variables
#finds the constant term in t (i.e. the coeff. of t^0)
#of the Laurent expansion of t, viewed as a formal Laurent
#series in the x's of xSet
#For example try CT1(1/(1-x*t)/(y-t),[x,y],t);
CT1:=proc(f,xSet,t)
local f1,i,lu,g,mu,BAD1,mu1,mu11,mu11a:
mu:=denom(normal(f)):

BAD1:={}:
if type(mu,`*`) then
for i from 1 to nops(mu) do
mu1:=op(i,mu):
if type(mu1, `^`) then
 mu11:=op(1,mu1):
else
mu11:=mu1:
fi:

mu11a:=expand(mu11/coeff(mu11,t,degree(mu11,t))):

if mu11a=t or isbad(mu11a,t,xSet) then
 BAD1:=BAD1 union {mu11}:
fi:
od:

else

mu1:=mu:
if type(mu1, `^`) then
 mu11:=op(1,mu1):
else
mu11:=mu1:
fi:

mu11a:=expand(mu11/coeff(mu11,t,degree(mu11,t))):

if mu11a=t or isbad(mu11a,t,xSet) then
 BAD1:=BAD1 union {mu11}:
fi:
fi:


f1:=convert(f,parfrac,t):

g:=f1:
for i from 1 to nops(f1) do
lu:=op(i,f1):
mu:=denom(lu):
if isbada(mu,BAD1,t) then
g:=g-lu:
fi:

od:


factor(normal(subs(t=0,g))):
end:
 

#GFx(a,Elist,x): Given a list of variables, a, 
#of type non-negative integer, and a list
#of linear-diophantine equations in these variables
#and a continuous indexed-variable x, outputs the
#generating function of all solutions 
#where x[1]^a[1]*x[2]^a[2]*... shows up iff (a[1],a[2], ...)
#is a solution
GFx:=proc(a,Eqs,x)
local f,i,t,xSet,g,fbot:
f:=GFxt(a,Eqs,x,t):
xSet:=[seq(x[i],i=1..nops(a))]:
for i from 1 to nops(Eqs) do
f:=CT1(f,xSet,t[i]):
od:

fbot:=denom(f):

for i from 1 to nops(xSet) do
 fbot:=subs(xSet[i]=0,fbot):
od:

if fbot=0 then
 RETURN(0):
fi:
g:=f:
for i from 1 to nops(a) do
g:=subs(x[i]=x[a[i]],g):
od:
g:
end:



#numfacs1(mekh,x): given a product, mekh, of terms
#counts how many of them depend on x
numfacs1:=proc(mekh,x)
local co,i:
co:=0:
for i from 1 to nops(mekh) do
if degree(op(i,mekh),x)>0 then
 co:=co+1:
fi:
od:
co:
end:


#ToUm1(R,x): Given a rational function R and a variable x
#tries to finds an umbra that sends f(x) to Constant
#term of R(x)f(1/x), it only treats the case where
#the denominator of R(x) factors completely over the 
#rationals and there are no multiplicity and
#the degree of the numerator is less than the degree of
#of the denominator
#if this is not the case, it returns 0
ToUm1:=proc(R,x)
local mekh,gu,i,mu,p:
mekh:=factor(denom(R)):

if degree(numer(R),x)>=degree(denom(R),x) then
 RETURN(0):
fi:

if type(mekh, `*`) and numfacs1(mekh,x)<>degree(mekh,x) then
 RETURN(0):
fi:

gu:={solve(mekh,x)}:

mu:={}:
readlib(residue):
for i from 1 to nops(gu) do
p:=normal(-residue(R,x=gu[i])/gu[i]):
mu:=mu union {[p,[0],[1/gu[i]]]}:
od:
mu:
end:




ToUm:=proc(R,xList)
local gu,mu,i,xList1,mui,mui1,mui2,mui3,x,lu,j:

if R=0 then
RETURN({}):
fi:

if nops(xList)=0 then
RETURN( {[R,[],[]]} ):
fi:

if nops(xList)=1 then
RETURN(ToUm1(R,xList[1])):
fi:

x:=xList[nops(xList)]:
xList1:=[op(1..nops(xList)-1,xList)]:
mu:=ToUm(R,xList1):

gu:={}:

for i from 1 to nops(mu) do
mui:=mu[i]:
mui1:=mui[1]:
mui2:=mui[2]:
mui3:=mui[3]:

for j from 1 to nops(mui3) do
if diff(mui3[j],x)<>0 then
 RETURN(0):
fi:
od:

lu:=ToUm1(mui1,x):

if lu=0 then
RETURN(0):
fi:

for j from 1 to nops(lu) do
gu:=gu union {[lu[j][1],[op(mui2),0],[op(mui3),op(lu[j][3])]]} :
od:
od:
gu:
end:

#######End programs taken from LinDiophantus


#CoreW(SymbWord): Given a symbolic word, SymbWord
#in terms of a list of pairs, [[l1,nu1],[l2,nu2], ...] extrats
#the first component, in other words, it contracts all repetitions
CoreW:=proc(SymbWord) local i: [seq(SymbWord[i][1],i=1..nops(SymbWord))]:end:

#CompatibleInterfaces1(u,v): given two (ordinary) words
#finds all the places in u such that the suffix of u starting
#at the i^th letter of u equals a prefix of v
CompatibleInterfaces1:=proc(u,v)
local gu,i,u1,v1:
gu:={}:
for i from 1 to nops(u) do
u1:=[op(i..nops(u),u)]:
if nops(v)>=nops(u1) then
 v1:=[op(1..nops(u1),v)]:
 if u1=v1 then
  gu:=gu union {i}:
 fi:
fi:
od:
gu:
end:

#CompatibleInterfaces(u,v,d): given two (ordinary) words
#finds all the places in u such that the suffix of u starting
#at the i^th letter of u equals a prefix of some image of v
#under the permutation of d letters
#it returns the set of pairs [place, image of v]
CompatibleInterfaces:=proc(u,v,d)
local gu,i,mu,lu,j,T,gu1:
gu:={}:
mu:=khaverim(u,d):

for i from 1 to nops(mu) do
 gu:=[op(gu),op(CompatibleInterfaces1(mu[i],v))]:
od:
gu:
gu1:=convert(gu,set):
for i from 1 to nops(gu1) do
T[gu1[i]]:=0:
od:
for i from 1 to nops(gu) do
T[gu[i]]:=T[gu[i]]+1:
od:
{seq([gu1[i],T[gu1[i]]],i=1..nops(gu1))};
end:


#FindUmbInt(U,V,i,x,t): Given two symbolic letters U, V
#expressed in symbolic frequency notation
#{[[let1,AffLinExp1],[let2,AffLinExp2], ...], {set_of_symbols}}
#and a location i, one of the legal interfaces of
#U before V in a symbolic cluster, and a sumbol x, 
#for the catalytic variabes and a symbol t, find the umbra that
#sends x_1^a_1*x_2^a_2 ... to the generating function of all possible
#such V's
#For example try: FindUmbInt([[[1,1],[2,a],[3,a],[1,1]],[a]],
#[[[1,1],[2,a],[3,a],[1,1]],[a]],4,x,t):
FindUmbInt:=proc(U,V,i,x,t)
local k1,k2,a,b,j,U1,V1,var,eq,lu,mu,m1,zanav,gu,mu1,c,j1,gu1:
eq:=[]:
var:=[]:
U1:=U:
for j from 1 to nops(U[2]) do
 U1:=subs(U[2][j]=a[j],U1):
od:

V1:=V:
for j from 1 to nops(V[2]) do
 V1:=subs(V[2][j]=b[j],V1):
od:

var:=[op(var),op(U1[2]),op(V1[2]),k1,k2]:

if i=nops(U1[1]) then
 var:=[op(var),m1]:
   if nops(V1[1])=1 then
    eq:=[op(eq),U1[1][i][2]-k1-m1-1, V1[1][1][2]-k2-m1-2]:
   else
    eq:=[op(eq),U1[1][i][2]-k1-m1-1 , V1[1][1][2]-k2-m1-1]:
   fi:

else

if i=1 then
eq:=[op(eq),U1[1][i][2]-(V1[1][1][2]+k1+1)]:
else
eq:=[op(eq),U1[1][i][2]-(V1[1][1][2]+k1)]:
fi:

for j from i+1 to nops(U1[1])-1 do
eq:=[op(eq), U1[1][j][2]-V1[1][j-i+1][2]]:
od:
j:=nops(U1[1]):
if j-i+1=nops(V1[1]) then
eq:=[op(eq), U1[1][j][2]+k2+1-V1[1][j-i+1][2]]:
else
eq:=[op(eq), U1[1][j][2]+k2-V1[1][j-i+1][2]]:
fi:
fi:



zanav:=0:
j:=nops(U1[1]):
for j1 from j-i+2 to nops(V1[1]) do
zanav:=zanav+V1[1][j1][2]:
od:

if j-i+1=nops(V1[1]) then
zanav:=zanav+1:
fi:

lu:=GFx(var,convert(eq,set),x):

mu:=[seq(x[op(j,U1[2])],j=1..nops(U1[2]))]:

gu:=ToUm(lu,mu):
gu:=SimplifyUmbra(gu):

if gu=0 then
 RETURN({}):
fi:
gu:=subs({x[k1]=1,x[m1]=1,x[k2]=t},gu):

mu1:=[seq(x[op(j,V1[2])],j=1..nops(V1[2]))]:
for j1 from 1 to nops(V1[2]) do
c:=coeff(zanav,op(j1,V1[2]),1):
zanav:=expand(zanav-c*op(j1,V1[2])):

gu:=subs(x[op(j1,V1[2])]=t^c*x[op(j1,V1[2])],gu):
od:
if not type(zanav,integer) then
 ERROR(`Something is wrong`):
fi:

gu:={seq([-t^zanav*gu[j1][1],gu[j1][2],gu[j1][3]],j1=1..nops(gu))}:

gu1:=gu:
for j1 from 1 to nops(mu1) do
gu1:=subs(mu1[j1]=x[j1],gu1):
od:
SimplifyUmbra(gu1):
end:


#FindUmb(U,V,x,t,d): Given two symbolic words, U, and V,
#and a variable x (x an index variable, for catalytic
#variables), and a variable t, to keep track of length
#computes the umbral edge connecting U and V
#including interfaces with the images of V
#d is the total number of letters
FindUmb:=proc(U,V,x,t,d)
local gu,lu,i,lu1:

lu:=CompatibleInterfaces(CoreW(U[1]),CoreW(V[1]),d):
gu:={}:
for i from 1 to nops(lu) do
lu1:=FindUmbInt(U,V,lu[i][1],x,t):
if lu1=0 then
 print(`Can't find an umbra`):
 RETURN(FAIL):
fi:

lu1:={seq([lu1[j][1]*lu[i][2],lu1[j][2],lu1[j][3]],j=1..nops(lu1))}:
 gu:=gu union lu1:
od:

SimplifyUmbra(gu):
end:

#IniWord(SymWord,x,t): Given a symbolic word, SymWord,
#and two variable symbols x (for the catalytic variables,
#x[1],x[2], ...) and t (to keep track of the length)
#finds the generating function for clusters consisting
#just of  SymWord
IniWord:=proc(SymWord,x,t)
local lu1,lu2,i,gu,mu,c:
lu1:=SymWord[1]:
lu2:=SymWord[2]:
gu:=0:

for i from 1 to nops(lu1) do
gu:=gu+lu1[i][2]:
od:
gu:=expand(gu):
mu:=1:

for i from 1 to nops(lu2) do
mu:=mu/(1-x[i]):
od:

for i from 1 to nops(lu2) do
c:=coeff(gu,lu2[i],1):
gu:=gu-c*lu2[i]:
mu:=subs(x[i]=t^c*x[i],mu):
od:

if not ( type(gu,integer) and gu>=0 ) then
 ERROR(`Something is wrong`):
fi:
-t^gu*mu:
end:


#WtUmSc(ListSW,x,t,d): Given a list of symbolic words
#ListSW, constructs a Weighted Umbral Scheme for the Cluster
#generating function where the mistakes are
#the members of ListSW and all their images under
#the symmetric group on d letters
#Recall that a Weighted Umbral Scheme is like an Umbral
#Scheme, but the first component is a list of
#weights rather then a set of 
WtUmSc:=proc(SetSW,x,t,d)
local U,V,j,gu1,gu2,gu3,gu4,i,gu2a:
gu1:=[seq(nops(khaverim(CoreW(SetSW[i][1]),d)) ,i=1..nops(SetSW))]:
gu2:=[]:

for i from 1 to nops(SetSW) do
U:=SetSW[i]:
gu2a:=[]:
 for j from 1 to nops(SetSW) do
 V:=SetSW[j]:
  gu2a:=[op(gu2a),FindUmb(U,V,x,t,d)]:
 od:
gu2:=[op(gu2),gu2a]:
od:

gu3:=[]:
gu4:=[]:

for i from 1 to nops(SetSW) do
 U:=SetSW[i]:
 gu3:=[op(gu3),IniWord(U,x,t)]:
 gu4:=[op(gu4),[seq(x[j],j=1..nops(U[2]))]]:
od:

gu2:=[seq([seq(gu2[j][i],j=1..nops(gu2[i]))],i=1..nops(gu2))]:

[gu1,gu2,gu3,gu4]:

end:


#UGJseries(SetSW,NumLetters,L): Given a set of symbolic
#letters, SetSW, and the number of letters, NumLetters,
#an integer, and an integer L, uses the Umbral Scheme
#to find the first L terms of the series expansion
#of words that avoid the mistakes of SetSW 
#and all their images  under the symmetric group
UGJseries:=proc(SetSW,NumLetters,L)
local gu,lu,x,t,Sid,i,mu:
gu:=WtUmSc(SetSW,x,t,NumLetters):
lu:={}:

for i from 1 to nops(gu[4]) do
lu:=lu union {op(gu[4][i])}:
od:

Sid:=ApplyUmSc1(gu,t,L,lu)[2]:
lu:=0:

for i from 1 to nops(Sid) do
lu:=lu+Sid[i]*t^i:
od:

lu:=1/(1-NumLetters*t-lu):
lu:=taylor(lu,t=0,L+1):
mu:=[]:
for i from 0 to L do
mu:=[op(mu),coeff(lu,t,i)]:
od:
mu:
end:



#######Stuff for checking using the finite Goulden-Jackson
##(transported from SYMGJ)
with(combinat):
SGJseries:=proc(NUMLETTERS,MISTAKES1,L)
local v,eq, var,i,lu,mu,pols,ku,MISTAKES,MISTAKESFULL:
eq:=[]:
var:=[]:
pols:=[]:
 
MISTAKES:={}:

for i from 1 to nops(MISTAKES1) do
MISTAKES:=MISTAKES union {canon(op(i,MISTAKES1))}:
od:



MISTAKESFULL:={}:

for i from 1 to nops(MISTAKES) do
MISTAKESFULL:=MISTAKESFULL union khaverim(op(i,MISTAKES),NUMLETTERS):
od:


for i from 1 to nops(MISTAKES) do
 v:=op(i,MISTAKES):
 mu:=findeqz_fast(v,MISTAKESFULL):
 eq:= [op(eq),mu[2]]:
 pols:= [op(pols),mu[1]]:
 var:=[op(var),C[op(v)]]:
od:
 

lu:=Ssolser(var,pols,eq,NUMLETTERS,L+2):


ku:=1-NUMLETTERS*s:

for i from 1 to L+2 do
 ku:=ku-op(i,lu)*s^i:
od:

ku:=taylor(1/ku,s=0,L+1):

lu:=[]:

for i from 0 to L do
lu:=[op(lu),coeff(ku,s,i)]:
od:

lu:
end:
 



#canon(mila):Given a word in an alphabet of integers, in terms of a list
#finds an image under the action of the symmetric group that
#is least lexicographically

canon:=proc(mila)
local i,gu,mu,lu:

gu:=[]:

for i from 1 to nops(mila) do
 if nops({op(1..i,mila)})>nops({op(1..i-1,mila)}) then
    gu:=[op(gu),op(i,mila)]:
 fi:
od:

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

lu:=[]:

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

lu:
end:


#khaverim(mila,n):Given a word mila finds the set of all its n(n-1)..(n-j+1)
#images under the symmetric group on {1,2, ..n},where j is the number of
#distinct letters in mila

khaverim:=proc(mila1,n)
local i,gu,r,lu,mila:
r:=nops(convert(mila1,set)):
mila:=canon(mila1):

if n<r then
 ERROR(`Not enough letters`):
fi:

gu:=permute(n,r):

lu:={}:

for i from 1 to nops(gu) do
lu:=lu union {peula(mila,op(i,gu))}:
od:

lu:

end:


#peula(mila,perm): Given a word in the alphabet {1,2,..,r} applies
#to it the (partial) permutation perm i.e. i is replaced by perm[i]

peula:=proc(mila,perm)
local gu,i:
gu:=[]:

for i from 1 to nops(mila) do
gu:=[op(gu),op(op(i,mila),perm)]:
od:

gu:
end:

 
#findeqz_fast sets up the equ C[v]= s+t*Sum_u overlap(u,v) *C[u]
#just like findeqz but returns seperately the polynomial part
#and the equation part
 
findeqz_fast:=proc(v,MISTAKES)
local eq,PO,i,u:

 PO:=(-1)*s^nops(v):

eq:=0: 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq-overlapz(u,v)*C[op(canon(u))]:
od:
 
PO,eq:
 
end:
 
 
#findeqz_fast sets up the equ C[v]= s+t*Sum_u overlap(u,v) *C[u]
#just like findeqz but returns seperately the polynomial part
#and the equation part
 
findeqz_fast:=proc(v,MISTAKES)
local eq,PO,i,u:

 PO:=(-1)*s^nops(v):

eq:=0: 
for i from 1 to nops(MISTAKES) do
 u:=op(i,MISTAKES):
 eq:=eq-overlapz(u,v)*C[op(canon(u))]:
od:
 
PO,eq:
 
end:
 

#overlapz is a procedure that given two words u and v
#computes the weight-enumerator of all v\suffix(u),
#for all suffixes of u that are prefixes of v, but with uniform weight s
overlapz:=proc(u,v)
local i,j,k,lu,gug:
 
lu:=0:
 
for i from 2 to nops(u) do
  for j from i to nops(u) while (j-i+1<=nops(v) and op(j,u)=op(j-i+1,v))
         do
   :
  od:   
 
if j-i=nops(v) and u<>v then
 ERROR(v,`is a subword of`,u,`illegal input`):
 fi:
 
 
 
 if j=nops(u)+1 and (i>1 or j>2) then
  gug:=1:
 
  for k from j-i+1 to nops(v) do
   gug:=gug*s:
  od:
 
  lu:=lu+gug:
 
fi:
 
od:
 
lu:
 
end:




Ssolser:=proc(vars,pols,equs,NUM,L)
local Mem,i,aj,maxi,X,r,lu,zm,gu, j, lu1, n, pij, v:

#Mem is a list of the memories that each a[j] needs


Mem:=[]:
for j from 1 to nops(vars) do
aj:=op(j,vars):

maxi:=0:

for i from 1 to nops(equs) do
 if degree(coeff(op(i,equs),aj,1),s)>maxi then
     maxi:=degree(coeff(op(i,equs),aj,1),s):
 fi:
od:

Mem:=[op(Mem),maxi]:
od:


for i from 1 to nops(vars) do
X[i]:=[seq(0,j=1..op(i,Mem))]:
od:

lu:=[]:

for n from 1 to L do
 for i from 1 to nops(equs) do
   gu[i]:=coeff(op(i,pols),s,n):

  for j from 1 to nops(vars) do
      aj:=op(j,vars):
      pij:=coeff(op(i,equs),aj,1):

       for r from 1 to degree(pij,s) do
         gu[i]:=gu[i]+coeff(pij,s,r)*op(nops(X[j])-r+1,X[j]):
       od:
  od:
 od:

lu1:=0:


for i from 1 to nops(vars) do
 v:=[op(op(i,vars))]:
 zm:=nops(convert(v,set)):
    lu1:=lu1+NUM!/(NUM-zm)!*gu[i]:
od:

lu:=[op(lu),lu1]:

for i from 1 to nops(vars) do
 X[i]:=[op(2..nops(X[i]),X[i]),gu[i]]:
od:

od:
lu:       
   
end:





#CheckU(SymbMisSet,Alphabet,L): Uses the Classical Goulden-
#Jackson to find the first L terms of the set of words in
#the alphabet Alphabet that avoid factors from SymbMisSet
#where each symbolic mistake has the form [[letter1,expression1], ...]
#and expression1... are affine linear expressions in the variables
#of DisVars
CheckU:=proc(SymbMisSet,Alphabet,L)
local MISTAKES1,i:
MISTAKES1:={}:
for i from 1 to nops(SymbMisSet) do
MISTAKES1:=MISTAKES1 union MistSetL(SymbMisSet[i][1],SymbMisSet[i][2],L):
od:
SGJseries(nops(Alphabet),MISTAKES1,L):
end:


#MistSetL(SymbMis,DisVars,L): Given a symbolic Mistake
#using the discrete variables of DisVars (denoting generic
#positive integers), finds all its manifestations
#of length up to L
MistSetL:=proc(SymbMis,DisVars,L)
local LinExp,i,mu,gu:
LinExp:=0:

for i from 1 to nops(SymbMis) do
 LinExp:= LinExp+SymbMis[i][2]:
od:

mu:=FeasVects(LinExp,DisVars,L):
gu:={}:
for i from 1 to nops(mu) do
 gu:=gu union {PlugSymbMis(SymbMis,DisVars,mu[i])}:
od:
gu:
end:



#PlugSymbMis(SymbMis,DisVars,PlugList): plugs the values
#in PlugList for the variables of DisVars, respectivelty
#into the symbolic mistake SymbMis
PlugSymbMis:=proc(SymbMis,DisVars,PlugList)
local gu,i,j,r,lu:
gu:=[]:

for i from 1 to nops(SymbMis) do
lu:=SymbMis[i][2]:
for j from 1 to nops(DisVars) do
lu:=subs(DisVars[j]=PlugList[j],lu):
od:
gu:=[op(gu),seq(SymbMis[i][1],r=1..lu)]:
od:
gu:
end:


#FeasVects(LinExp,DisVars,L): Given a linear expression
#in the list of discrete variables (denoting generic
#non-negative integers), finds all solutions
#for which LinExp at DisVars is <=L
FeasVects:=proc(LinExp,DisVars,L)
local gu,a,i,mu,L1,DisVars1,gu1,LinExp1,j:
gu:={}:

if nops(DisVars)=0 then
RETURN({[]}):
fi:

if nops(DisVars)=1 then
a:=DisVars[1]:
if coeff(LinExp,a,1)=0 then
ERROR(`LinExp does not involve `,a):
fi:
for i from 0 while subs(a=i,LinExp)<=L do
gu:=gu union {[i]}:
od:
RETURN(gu):
fi:

a:=DisVars[nops(DisVars)]:
mu:=coeff(LinExp,a,1)*a:
gu:={}:
for i from 0 while subs(a=i,mu)<=L do
L1:=L-subs(a=i,mu):
DisVars1:=[op(1..nops(DisVars)-1,DisVars)]:
LinExp1:=LinExp-mu:
gu1:=FeasVects(LinExp1,DisVars1,L1):

for j from 1 to nops(gu1) do
gu:=gu union {[op(gu1[j]),i]}:
od:

od:

gu:

end:

#####End part for checking that was transported from SymGJ