#######################################################################
## ROTA: Save this file as ROTA. To use it, stay in the               #
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read ROTA : <Enter>                                #
## Then follow the instructions given there                           #
##                                                                    #
## Written by Doron Zeilberger, Temple University ,                   #
##  zeilberg@math.temple.edu.                                         # 
#######################################################################
#This accompanies Doron Zeilberger's series of articles:
#The Umbral Transfer-Matrix Method. Most specifically, part I
#Created: March 15, 2000
#This version: March 30, 2000
#ROTA: A Maple package to handle Umbral Schemes
#Please report bugs to zeilberg@math.temple.edu
print(`Created:  March 15, 2000.`):
print(`This version: March 30,  1999`):
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 MAIN procedures type ezra(), for help with`):
 print(`a specific procedure, type ezra(procedure_name)`):
 
 print(``):
 
ezra1:=proc()
if args=NULL then
 print(`Contains the following procedures:`):
 print(` ApplyUmbra, ApplyUmbralMatrix`):
 print(`  ApplyUmSc, DataSet, MatrixUmbra, `):
 print(` SolveUmSch2, ToUmbra, Umbra, UmSchEq, YafeUm `):
fi:
end:
 
ezra:=proc()
if args=NULL then
 print(` The Main  procedure is: `):
 print(`  ApplyUmSc `):
 lprint(` `):
 print(`For help with it type: ezra(ApplyUmSc) ; `):
 print(``):
 print(` To see the names of the included examples of Umbral Schemes,`):
 print(` and their descriptions, type ezra(DataSet); `):
 print(``):
 print(`For a list of the other procedures type: ezra1();`):
fi:
 
 
if nops([args])=1 and op(1,[args])=DataSet  then
print(`We have the following examples Umbral Schemes:`):
print(` UmSch1, UmSch2, UMW, UMP1, UMP2 `):
lprint(``):
print(`UmSch1 is an Umbral Scheme for 1-board polyominoes`): 
print(` i.e. vertically-convecx polyominoes `):
print(` (using the catalyst-variable x[1] and the static variable q `):
print(` to weight-enumerate the number of cells ) `):
lprint(``):
print(`UmSch2 is an Umbral Scheme for 2-board polyominoes`): 
print(` (using the catalyst-variables x[1],y[1],x[2], and the `):
print(` static variable q `):
print(` (to weight-enumerate the number of cells ) `):
lprint(``):
print(`UMP1 is the Umbral Scheme for ordinary partitions (allowing 0)`):
print(` using the catalyst-variable x (corresponding to largest part) and `):
print(` the static variables t (number of parts) and q (sum of parts) `):
lprint(``):
print(`UMP2 is the Umbral Scheme for 2-rowed plane-partitions (allowing 0s)`):
print(` using the catalyst-variable x1,x2 `):
print(`corresponding to rightmost column and `):
print(` the static variables t (number of columns) and q (sum of parts) `):
lprint(``):
print(` UMW is the Umbral Scheme given as the last example in the article:`):
print(` The Umbral Transfer-Matrix Method: I. Foundations `):
 
 
fi:
 
if nops([args])=1 and op(1,[args])=UmSchEqEx  then
print(`UmSchEqEx(UmSch,F): Given an umbral scheme, UmSch, and a vector`):
print(`F of expressions ,finds`):
print(`the set of equations in terms of that make it annihilated`):
fi:
 
 
if nops([args])=1 and op(1,[args])=UmSchEq  then
print(` UmSchEq(UmSch,F,a): Given an umbral scheme, UmSch, a lettter `):
print(` F to denote the functions and a dummy letter a, finds `):
print(` the set of functional equations in terms of F[1](args[1]),  `):
print(` F[2](args[2]), corresponding to the letters `):
fi:
 
if nops([args])=1 and op(1,[args])=ApplyUmSc  then
print(`ApplyUmSc(UmSch,q,n,vars): Given an Umbral Scheme, UmSch, finds`):
print(`the generating functions of creatures up to weight n`):
print(`followed by substition of all the variables of var to 1 and summing`):
print(`over leftmost letters`):
fi:
 
if nops([args])=1 and op(1,[args])=SerToPol then
print(`SerToPol(sidra,q,n): given a series, sidra, finds the polynomial`):
print(`consisting of the terms up to degree n`):
fi:
 
if nops([args])=1 and op(1,[args])=SimplifyUmbra then
print(` SimplifyUmbra(Umb): Given an Umbra, collects like terms `):
print( ` and returns a simplified version `):
fi:
 
if nops([args])=1 and op(1,[args])=ConvertToUmbra2 then
print(`ConvertToUmbra2(Umb): Given an umbra Umb in the 3-list-format`):
print(`i.e. as a set of 3-lists of the form`):
print(`[FRONT,Orders_Of_Derivaties,Substitutions]`):
print(`converts this to a set of 2-lists of the form `):
print(`[FRONT,[Orders_Of_Derivaties,Substitutions]]`):
fi:
 
 
if nops([args])=1 and op(1,[args])=ApplyUmbralMatrix then
print(` ApplyUmbralMatrix(UmbraM,Fvec,ylistvec): Given an Umbral matrix`):
print(` UmbraM (in terms of a list of lists), a vector of expressions,`):
print(` Fvec, and a vector to variable-lists corresponding to the`):
print(` components of UmbraM and Fvec`):
fi:
 
if 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`):
fi:
 
if nops([args])=1 and op(1,[args])=ToUmbra then
print(` ToUmbra(poly1,xlist,alist): given a polynomial, poly1, in the`):
print(` variables xlist with exponents that are affine-linear`):
print(` expressions in the discrete variables in the`):
print(` list alist, and in the variables of alist themselves`):
print(` outputs the corresponding Umbra, such that`):
print(` poly1 is the image, under that umbra, of the`):
print(` generic monomial y[1]^alist[1]*...*y[k]^alist[k],`):
print(` (where k=nops(alist), and y[1], ..., y[k] are generic`):
print(` continuous variables that correspond to the disrete variables`):
print(` in alist`):
print(` The format of the output is a set, each element of which`):
print(` is a list of 3-elements`):
print(` [FRONT, Diffs,SubsList], where the FRONT`):
print(` is a rational function in whatever, Diffs is a list of`):
print(` integers of the same length as alist, and SubsList is`):
print(` the list of substitutions that the continuous variables`):
print(` that correspond to alist have to be substituted by`):
print(` For example ToUmbra(a*x^b+b*y^a,[x,y],[a,b]) should yield`):
print(`  {[1, [1, 0], [1, x]], [1, [0, 1], [y, 1]]}`):
fi:
 
 
 
if nops([args])=1 and op(1,[args])=PreUmbra1 then
print(`PreUmbra1(Interface,m,n,VarList,WtList): given an interface`):
print(`Interface of two dynasties of length m and n resp.`):
print(`(the m dynasty standing to the left of the n dysnasty)`):
print(` and the left dynasty m, is held fixed, and the description of`):
print(`the interface is in the form [[a1,b1],[a2,b2],...,[ak,bk]],`):
print(`where the rth pair [ar,br] desribe the lowest (ar) and highest`):
print(` (br) kings of the right dynasty that it has a non-empty`):
print(` overlap with)`):
print(`(where m and n are positive integers), VarList is a list of`):
print(`symbols, of length n-2 and WtList is a list of weights of`):
print(`length m-2`):
print(``):
print(`For example, PreUmbra1([[1,2],[2,2],[2,3]],3,3,[a],[q*x]);`):
fi:
 
if nops([args])=1 and op(1,[args])=schum then
print(`schum(a,reshX): Given a symbol a, and a list of epressions`):
print(`reshX=[X1,...,Xr], finds the expression, featuring a,`):
print(`of the sum of X1^a1*...*Xr^ar over all tuples`):
print(`(a1,...,ar) of POSITIVE integers that sum to a,`):
print(`i.e. a1+...+ar=a`):
fi:
 
if nops([args])=1 and op(1,[args])=YafeUm then
print(`YafeUm(Umb1,F,D): Given an umbral operator, Umb1, in`):
print(`Maple notation, converts it to human notation`):
print(`using the symbol F for the function and D for differntiation`):
fi:
 
end:
 
###Begining of Data
 
UMP1:=
[ {1}, [[{[1/(1-q*x), [0], [q*x]]}]], [q*x/(1-q*x)], [[x]]]:
 
UMP2:=
[ {1},
[[
{[t/(1-q*x1)/(1-q*x2), [0,0], [q*x1,q*x2]],
[-t/(1-q*x1)/(1-q*x2), [0,0], [q^2*x1*x2,1]],
[t/(1-q*x1)/(1-q^2*x1*x2), [0,0], [q^2*x1*x2,1]]}
]], 
[t/(1-q*x1)/(1-q^2*x1*x2)], [[x1,x2]]]:
 
UMW:=[{1,2},
[
[ 
{[1/(1-q*x),[0],[q*x]]},
{[1/(1-q*x),[0],[q*x]]}
],
[{[q*x/(1-q*x),[0],[1]],[-1/(1-q*x),[0],[q*x]]},{}]
]
,
[q*x/(1-q*x),0],[[x],[x]]]:
 
UmSch1:=
[  {1}, [[{[q^2*x[1]^2/(q*x[1]-1)^2, [0], [1]], [-1/(q*x[1]-1)*q*x[
1], [1], [1]]}]], [q*x[1]/(1-q*x[1])], [[x[1]]]]:
 
UmSch2:=
[{1, 2}, [[{[-q*x[1]/(q*x[1]-1), [1], [1
]], [q^2*x[1]^2/(q*x[1]-1)^2, [0], [1]]}, {[-q*x[1]/(q*x[1]-1), [1, 0, 0], [1,
1, 1]], [-q*x[1]/(q*x[1]-1), [0, 0, 1], [1, 1, 1]], [-q^2*x[1]^2/(q*x[1]-1)^2,
[0, 0, 0], [1, q*x[1], 1]], [2*q^2*x[1]^2/(q*x[1]-1)^2, [0, 0, 0], [1, 1, 1]]},
{[q^2*x[1]^2/(q*x[1]-1)^2, [0, 0, 0], [1, q*x[1], 1]]}], [{[-q^2*x[1]*x[2]*y[1]
/(-1+q*x[2])/(-1+y[1])/(q*x[1]-1), [1], [1]], [q^2*x[1]*x[2]/(q*x[1]-1)/(-1+y[1
])^2/(-1+q*x[2]), [0], [y[1]]], [q^2*x[1]*x[2]*y[1]*(y[1]-2)/(q*x[1]-1)/(-1+y[1
])^2/(-1+q*x[2]), [0], [1]]}, {[-y[1]*x[2]^2*x[1]*q^2*(-q*x[2]+q*x[1]-1+y[1])/(
q*x[1]-1)/(-1+y[1])/(q*x[2]-y[1])/(-x[2]+x[1])/(-1+q*x[2]), [0, 0, 0], [1, q*x[
2], 1]], [2*q^2*x[1]*x[2]*y[1]*(y[1]-2)/(q*x[1]-1)/(-1+y[1])^2/(-1+q*x[2]), [0,
0, 0], [1, 1, 1]], [q^2*x[1]*x[2]/(q*x[1]-1)/(-1+y[1])^2/(-1+q*x[2]), [0, 0, 0]
, [y[1], y[1], y[1]]], [-x[2]*q^2*x[1]^2*y[1]*(-y[1]-q*x[2]+q*x[1]+1)/(-y[1]+q*
x[1])/(q*x[1]-1)/(-1+y[1])/(-x[2]+x[1])/(-1+q*x[2]), [0, 0, 0], [1, q*x[1], 1]]
, [-q^2*x[1]*x[2]*y[1]/(-1+q*x[2])/(-1+y[1])/(q*x[1]-1), [1, 0, 0], [1, 1, 1]],
[q^3*x[1]*x[2]^2/(q*x[1]-1)/(-1+y[1])/(q*x[2]-y[1])/(-1+q*x[2]), [0, 0, 0], [y[
1], q*x[2], 1]], [q^2*x[1]*x[2]/(q*x[1]-1)/(-1+y[1])^2/(-1+q*x[2]), [0, 0, 0],
[y[1], 1, 1]], [-q^2*x[1]*x[2]*y[1]/(-y[1]+q*x[1])/(-1+y[1])^2/(-1+q*x[2]), [0,
0, 0], [1, y[1], y[1]]], [q^2*x[1]*x[2]/(q*x[1]-1)/(-1+y[1])^2/(-1+q*x[2]), [0,
0, 0], [1, 1, y[1]]], [q^3*x[1]^2*x[2]/(-y[1]+q*x[1])/(-1+q*x[2])/(-1+y[1])/(q*
x[1]-1), [0, 0, 0], [1, q*x[1], y[1]]], [-q^2*x[1]*x[2]*y[1]/(-1+q*x[2])/(-1+y[
1])/(q*x[1]-1), [0, 0, 1], [1, 1, 1]], [y[1]^2*x[1]*q^2*x[2]/(-y[1]+q*x[1])/(-1
+y[1])^2/(q*x[2]-y[1]), [0, 0, 0], [1, y[1], 1]], [-q^2*x[1]*x[2]*y[1]/(q*x[1]-\
1)/(-1+y[1])^2/(q*x[2]-y[1]), [0, 0, 0], [y[1], y[1], 1]]}, {[q^2*x[1]*x[2]*y[1
]/(q*x[1]-1)/(q*x[2]-y[1])/(-1+q*x[2])^2, [0, 0, 0], [q*x[2], q*x[2], 1]], [-q^
3*x[1]^2*x[2]/(-y[1]+q*x[1])/(-1+q*x[2])/(-1+y[1])/(q*x[1]-1), [0, 0, 0], [1, q
*x[1], y[1]]], [-q^3*x[1]*x[2]^2/(q*x[1]-1)/(-1+y[1])/(q*x[2]-y[1])/(-1+q*x[2])
, [0, 0, 0], [y[1], q*x[2], 1]], [q^2*x[1]*x[2]*y[1]/(q*x[1]-1)^2/(-y[1]+q*x[1]
)/(-1+q*x[2]), [0, 0, 0], [1, q*x[1], q*x[1]]], [q^3*x[1]*y[1]*x[2]^2/(q*x[1]-1
)/(-1+y[1])/(-1+q*x[2])^2, [0, 0, 0], [1, q*x[2], 1]], [q^3*x[1]^2*y[1]*x[2]/(q
*x[1]-1)^2/(-1+y[1])/(-1+q*x[2]), [0, 0, 0], [1, q*x[1], 1]]}], [{[q^2*x[1]*x[2
]*y[1]*(x[2]*y[1]*q+4-3*q*x[2]-3*q*x[1]-2*y[1]+2*q^2*x[1]*x[2]+y[1]*q*x[1])/(-1
+y[1])^2/(q*x[1]-1)^2/(-1+q*x[2])^2, [0], [1]], [-2*q^2*x[1]*x[2]/(q*x[1]-1)/(-\
1+y[1])^2/(-1+q*x[2]), [0], [y[1]]]}, {[2*q^2*x[1]*x[2]*y[1]*(x[2]*y[1]*q+4-3*q
*x[2]-3*q*x[1]-2*y[1]+2*q^2*x[1]*x[2]+y[1]*q*x[1])/(-1+y[1])^2/(q*x[1]-1)^2/(-1
+q*x[2])^2, [0, 0, 0], [1, 1, 1]], [-2*y[1]^2*x[1]*q^2*x[2]/(-y[1]+q*x[1])/(-1+
y[1])^2/(q*x[2]-y[1]), [0, 0, 0], [1, y[1], 1]], [-2*q^3*x[1]^2*x[2]/(-y[1]+q*x
[1])/(-1+q*x[2])/(-1+y[1])/(q*x[1]-1), [0, 0, 0], [1, q*x[1], y[1]]], [-2*q^3*x
[1]*x[2]^2/(q*x[1]-1)/(-1+y[1])/(q*x[2]-y[1])/(-1+q*x[2]), [0, 0, 0], [y[1], q*
x[2], 1]], [(-q^2*x[2]^2+q^2*x[1]*x[2]+y[1]*q*x[1]-2*q*x[1]+x[2]*y[1]*q+2-2*y[1
])*q^2*x[1]*x[2]^2*y[1]/(-1+q*x[2])^2/(-x[2]+x[1])/(q*x[2]-y[1])/(q*x[1]-1)/(-1
+y[1]), [0, 0, 0], [1, q*x[2], 1]], [-2*q^2*x[1]*x[2]/(q*x[1]-1)/(-1+y[1])^2/(-\
1+q*x[2]), [0, 0, 0], [1, 1, y[1]]], [-2*q^2*x[1]*x[2]/(q*x[1]-1)/(-1+y[1])^2/(
-1+q*x[2]), [0, 0, 0], [y[1], y[1], y[1]]], [-2*q^2*x[1]*x[2]/(q*x[1]-1)/(-1+y[
1])^2/(-1+q*x[2]), [0, 0, 0], [y[1], 1, 1]], [2*q^2*x[1]*x[2]*y[1]/(-y[1]+q*x[1
])/(-1+y[1])^2/(-1+q*x[2]), [0, 0, 0], [1, y[1], y[1]]], [(-q^2*x[1]*x[2]+x[1]^
2*q^2+2*q*x[2]-x[2]*y[1]*q-y[1]*q*x[1]-2+2*y[1])*q^2*x[1]^2*x[2]*y[1]/(-1+q*x[2
])/(-x[2]+x[1])/(q*x[1]-1)^2/(-y[1]+q*x[1])/(-1+y[1]), [0, 0, 0], [1, q*x[1], 1
]], [2*q^2*x[1]*x[2]*y[1]/(q*x[1]-1)/(-1+y[1])^2/(q*x[2]-y[1]), [0, 0, 0], [y[1
], y[1], 1]]}, {[-y[1]*x[2]^2*x[1]*q^2*(-q*x[2]+q*x[1]-1+y[1])/(q*x[1]-1)/(-1+y
[1])/(q*x[2]-y[1])/(-x[2]+x[1])/(-1+q*x[2]), [0, 0, 0], [1, q*x[2], 1]], [q^2*x
[1]*x[2]/(q*x[1]-1)/(-1+y[1])^2/(-1+q*x[2]), [0, 0, 0], [y[1], y[1], y[1]]], [-
q^2*x[1]*x[2]*y[1]/(q*x[1]-1)^2/(-y[1]+q*x[1])/(-1+q*x[2]), [0, 0, 0], [1, q*x[
1], q*x[1]]], [-x[2]*q^2*x[1]^2*y[1]*(-y[1]-q*x[2]+q*x[1]+1)/(-y[1]+q*x[1])/(q*
x[1]-1)/(-1+y[1])/(-x[2]+x[1])/(-1+q*x[2]), [0, 0, 0], [1, q*x[1], 1]], [2*q^3*
x[1]*x[2]^2/(q*x[1]-1)/(-1+y[1])/(q*x[2]-y[1])/(-1+q*x[2]), [0, 0, 0], [y[1], q
*x[2], 1]], [2*q^3*x[1]^2*x[2]/(-y[1]+q*x[1])/(-1+q*x[2])/(-1+y[1])/(q*x[1]-1),
[0, 0, 0], [1, q*x[1], y[1]]], [-q^2*x[1]*x[2]*y[1]/(-y[1]+q*x[1])/(-1+y[1])^2/
(-1+q*x[2]), [0, 0, 0], [1, y[1], y[1]]], [-q^2*x[1]*x[2]*y[1]/(q*x[1]-1)/(q*x[
2]-y[1])/(-1+q*x[2])^2, [0, 0, 0], [q*x[2], q*x[2], 1]], [y[1]^2*x[1]*q^2*x[2]/
(-y[1]+q*x[1])/(-1+y[1])^2/(q*x[2]-y[1]), [0, 0, 0], [1, y[1], 1]], [-q^2*x[1]*
x[2]*y[1]/(q*x[1]-1)/(-1+y[1])^2/(q*x[2]-y[1]), [0, 0, 0], [y[1], y[1], 1]]}]],
[q*x[1]/(1-q*x[1]), 0, q^2*x[1]/(1-q*x[1])*y[1]/(1-y[1])*x[2]/(1-q*x[2])], [[x[
1]], [x[1], y[1], x[2]], [x[1], y[1], x[2]]]]:
 
###End of Data
 
 
 
 
# 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:
 
 
# 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:
 
 
# 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:
 
 
# 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(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+fu[kv[j]]:
od:
 
 
ku:={seq(var[gc]=1,gc=1..nops(var))}:
mu:=subs(ku,mu):
 
Sid:=[op(Sid),factor(normal(mu))]:
 
lprint(Sid):
 
od:
 
RETURN(lu,Sid):
 
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:
 
#This version of grab only works for Maple 3
#grab3(mono,xlist,a): given a monomial, mono, in the variables 
#that are in the list of variables xlist, and a single 
# exponent (discrete) variable, a, returns
#the largest monomial, lu, such that lu^a is a factor of mono
grab3:=proc(mono,xlist,a)
local mono1,mono2,i,z,lu:
 
mono1:=expand(mono):
mono2:=mono1:
mono2:=simplify(mono2,symbolic):
for i from 1 to nops(xlist) do
 mono2:=subs(xlist[i]^a=z[i],mono2):
od:
 
lu:=1:
for i from 1 to nops(xlist) do
lu:=lu*xlist[i]^degree(mono2,z[i]):
od:
 
lu, normal(expand(mono1/lu^a)):
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:
 
 
#schum(a,reshX): Given a symbol a, and a list of epressions
#reshX=[X1,...,Xr], finds the expression, featuring a,
#of the sum of X1^a1*...*Xr^ar over all tuples
#(a1,...,ar) of POSITIVE integers that sum to a,
#i.e. a1+...+ar=a
schum:=proc(a,reshX)
local ak,mu,b,k,Xk:
 
k:=nops(reshX):
 
if k<1 then
ERROR(`List must have at least one element`):
fi:
 
if k=1 then
RETURN(reshX[1]^a):
fi:
 
Xk:=reshX[k]:
mu:=schum(b,[op(1..k-1,reshX)]):
normal(expand(simplify(sum(Xk^ak*subs(b=a-ak,mu),ak=1..a-(k-1))))):
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:
 
#SimplifyUmbra(Umb): Given an Umbra, collects like terms
#and returns a simplified version
SimplifyUmbra:=proc(Umb)
local Umb1,gu,kv,i,mu,evar,F:
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):
mu:=mu union {[normal(coeff(gu,F(evar))),evar[1],evar[2]]}:
od:
 
mu:
 
end:
 
# 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({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,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:
 
 
#UmSchEqEx(UmSch,F): Given an umbral scheme, UmSch, and a vector
#F of expressions ,finds
#the set of equations in terms of that make it annihilated
 
UmSchEqEx:=proc(UmSch,F)
local MAT,INI,LiVar,lu,i:
 
MAT:=UmSch[2]:
INI:=UmSch[3]:
LiVar:=UmSch[4]:
 
lu:=ApplyUmbralMatrix(MAT,F,LiVar):
 
[seq(numer(normal(F[i]-lu[i]-INI[i])),i=1..nops(LiVar))]:
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:
 
 
#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:
 
#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: