The main purpose of this package is to provide tools for manipulating combinatorial (Hopf) algebras. To setup the stage, we start this guided tour by presenting a few sample computations with two examples of such algebras. Then, we proceed by illustrating with many examples the predefined combinatorial objects and how to define new ones and the predefined combinatorial algebras and how to define new ones. We conclude this tour by a summary of the current features.
We load the MuPAD-Combinat package, and define a shortcut for the algebra of symmetric functions [MacDonald.SF.1995]:
>> package("Combinat", Quiet):
S := examples::SymmetricFunctions():
We consider the three first elementary symmetric polynomials in
the variables {x1,...,x6}
:
>> alphabet := [ x1, x2, x3, x4, x5, x6]:
e1 := expand(S::e([1])(alphabet));
x1 + x2 + x3 + x4 + x5 + x6
>> e2 := expand(S::e([2])(alphabet));
x1 x2 + x1 x3 + x1 x4 + x2 x3 + x1 x5 + x2 x4 + x1 x6 + x2 x5 + x3 x4 + x2 x6 + x3 x5 + x3 x6 + x4 x5 + x4 x6 + x5 x6
>> e3 := expand(S::e([3])(alphabet))
x1 x2 x3 + x1 x2 x4 + x1 x2 x5 + x1 x3 x4 + x1 x2 x6 + x1 x3 x5 + x2 x3 x4 + x1 x3 x6 + x1 x4 x5 + x2 x3 x5 + x1 x4 x6 + x2 x3 x6 + x2 x4 x5 + x1 x5 x6 + x2 x4 x6 + x3 x4 x5 + x2 x5 x6 + x3 x4 x6 + x3 x5 x6 + x4 x5 x6
Computing the product of two such polynomials yields a huge polynomial which is not quite practical to manipulate:
>> expand(e2*e3)
10 x1 x2 x3 x4 x5 + 10 x1 x2 x3 x4 x6 + 10 x1 x2 x3 x5 x6 + 10 x1 x2 x4 x5 x6 + 10 x1 x3 x4 x5 x6 + 10 x2 x3 x4 x5 x6 + 2 3 x1 x2 x3 x4 + ... (one page of output) 2 2 x1 x2 x3 + ... (another page of output)
Instead, if we use the symmetries, the previous product can be expressed as compactly as:
>> S::m( S::e([2]) * S::e([3]) );
10 m[1, 1, 1, 1, 1] + 3 m[2, 1, 1, 1] + m[2, 2, 1]
Here, m[2, 1, 1, 1]
denotes the monomial symmetric function
m2,1,1,1 obtained by summing all the monomials with one
variable elevated to the power 2 and three variables to the power
1.
Now, we are working in an algebra whose basis is indexed by partitions, and we want to compute efficiently in this algebra.
As another typical example, we are currently working on the so-called Loday-Ronco algebra [Loday_Ronco.1998], which is in particular of interest for theoretical physicists [Brouder_Frabetti.2001,Brouder,Foissy.thesis]. It is implemented as a combinatorial algebra having binary trees as basis.
>> LRA:= (examples::LodayRoncoAlgebra())::p:
For example, take the two following trees:
>> t1 := LRA(combinat::binaryTrees::unrank(6, 4))
o / \ \
>> t2 := LRA(combinat::binaryTrees::unrank(26, 5))
o / \ / \
You can make a formal linear combination of them:
>> t3 := 2*t2 + 3/4*t1
o o 2 / \ + 3/4 / \ / \ \
or take their product:
>> t2*t1
o o o o o / \ / \ / \ / \ / \ / \ \ + / \ / \ + / \ \ + / \ / \ + / \ \ + / \ \ \ / \ \ / \ / \ / \ o / \ / \ / \ / \
Here is a more complicated product in this algebra:
>> t3*t2
o o o / \ / \ / \ \ \ / \ 3/4 \ + 3/4 / \ + 3/4 \ + / \ \ \ / \ / \ / \ o o o / \ / \ / \ / \ \ / \ 3/4 \ + 3/4 / \ + 3/4 \ + \ / \ / \ / \ \ \ o o o / \ / \ / \ / \ / \ / \ 3/4 \ + 3/4 / \ + 3/4 / \ + / \ \ \ \ \ \ o o o / \ / \ / \ / \ \ / \ 3/4 / \ + 3/4 / \ + 3/4 \ + \ / \ / \ \ / / o o o / \ / \ / \ / \ / \ / \ 3/4 \ + 3/4 / \ + 3/4 / \ + / \ \ \ / / / o o o / \ / \ / \ / \ / \ / \ 3/4 / \ + 3/4 / \ + 3/4 / \ + \ / / / \ \ o o / \ / \ o / \ / \ / \ 3/4 / \ + 3/4 / + 2 / \ \ + / / \ / \ \ \ / \ o o o o / \ / \ / \ / \ 2 / \ / \ + 2 / \ + 2 / \ / \ + 2 / \ + \ / \ \ / \ / \ / \ / \ / \ \ \ o o o o / \ / \ / \ / \ 2 / \ + 2 / \ / \ + 2 / \ + 2 / \ + / \ / \ / \ / \ / \ / \ \ / / / \ / o / \ / \ 2 / / \ / \
We now describe in more details how all of this works. In the following, we assume that the package MuPAD-Combinat has been loaded into MuPAD and, for shortening the notations, that the library combinat has been exported. We also assume that the reader is somewhat familiar with the MuPAD syntax, and some of the exercises requires acquaintance with more advanced features of MuPAD like domain programming;we refer to the MuPAD tutorial for details. Technicalities can be safely ignored in a first reading; they will be better understood after the explanations in the design notes
.>> package("Combinat", Quiet):
export(combinat):
We start by some sample applications at random. We compute the
first terms of the famous Catalan sequence, we generate the
Cartesian product of three lists, we compute all permutations of the
numbers 1, 2, 3, and we ask for all sub-words of the word [a,
b, c, d]
:
>> catalan(i) $ i = 0..10
1, 1, 2, 5, 14, 42, 132, 429, 1430, 4862, 16796
>> cartesianProduct::list([1,2,3],[a,b],[i,ii,iii])
[[1, a, i], [1, a, ii], [1, a, iii], [1, b, i], [1, b, ii], [1, b, iii], [2, a, i], [2, a, ii], [2, a, iii], [2, b, i], [2, b, ii], [2, b, iii], [3, a, i], [3, a, ii], [3, a, iii], [3, b, i], [3, b, ii], [3, b, iii]]
>> permutations::list([1, 2, 3])
[[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]]
>> subwords::list([a,b,c,d])
[[], [a], [b], [c], [d], [a, b], [a, c], [a, d], [b, c], [b, d], [c, d], [a, b, c], [a, b, d], [a, c, d], [b, c, d], [a, b, c, d]]
We turn now to various combinatorial classes; in short a
combinatorial class is a set of related combinatorial objects,
like the set of all integer partitions. For each such class, there
is a sub-library of combinat
. We can use the library
combinat::partitions
to list all the integer partitions of
5:
>> partitions::list(5)
[[5], [4, 1], [3, 2], [3, 1, 1], [2, 2, 1], [2, 1, 1, 1], [1, 1, 1, 1, 1]]
Let us draw the partition [3,2]
using boxes
(French or Cartesian notation):
>> partitions::printPretty([3, 2])
+---+---+ | | | +---+---+---+ | | | | +---+---+---+
Now, we can fill those boxes with the numbers 1,2,3,4,5
so
that the numbers are increasing along rows and columns to obtain
so-called standard tableaux. Here are all the standard
tableaux of shape [3,2]
:
>> map(tableaux::list([3, 2]), tableaux::printPretty)
-- +---+---+ +---+---+ +---+---+ +---+---+ | | 4 | 5 | | 3 | 5 | | 2 | 5 | | 3 | 4 | | +---+---+---+ +---+---+---+ +---+---+---+ +---+---+---+ | | 1 | 2 | 3 |, | 1 | 2 | 4 |, | 1 | 3 | 4 |, | 1 | 2 | 5 |, -- +---+---+---+ +---+---+---+ +---+---+---+ +---+---+---+ +---+---+ -- | 2 | 4 | | +---+---+---+ | | 1 | 3 | 5 | | +---+---+---+ --
Ordered trees are another typical combinatorial class. Here are all the trees on four nodes:
>> trees::list(4)
-- o -- | o o o o | | | /|\, / \, / \, | , | | -- | | / \ | --
and here are a lot of trees:
>> trees::list(6)
-- | o o | o o o o /|\ o o o /|\ | //|\\, // \\, // \\, /| \ , |, // \\, /|\, / | \, | , | | | / \ | | || / \ | -- o o o o o / \ / \ / \ / \ o o o o / \ , / \, / \, | , |, // \\, /|\, /|\, / \ , /|\ | | / \ | | | | || | / \ | o o o o o o / \ o o /|\ / \ o / \ / \ / \ | |, / |\, / \, | , | |, / \, / \ , / \ , | , | / \ / \ | | | /|\ | | / \ o o o o o o o o o / \ o | | | | | | | | | , | , /|\, /|\, / \ , / \, /|\, / \, / \, / \, | // \\ | | / \ | | | | / \ | | | | o -- o o o o | | | | | | | | | , | , | , | , | | /|\ / \ / \ | | | | | / \ | --
All the sub-libraries of combinat
share a standardized
interface. Let us look in more detail at the library
combinat::partitions
. We can count partitions:
>> partitions::count(i) $ i = 0..10
1, 1, 2, 3, 5, 7, 11, 15, 22, 30, 42
list them under some extra conditions (here we list the partitions of 5 whose length is between 2 and 3):
>> partitions::list(5, MinLength = 2, MaxLength = 3);
[[4, 1], [3, 2], [3, 1, 1], [2, 2, 1]]
or compare them (lexicographically):
>> bool(partitions::_less([3, 1], [2, 2]))
FALSE
An important feature of MuPAD-Combinat are the so-called
generators, which allow to run through huge lists of
combinatorial objects without expanding the full lists into memory.
Technically, a generator is a function g
such that each call
g()
returns either a new object, or FAIL if no more objects
are available. Let us build a generator for the partitions of 4:
>> g := partitions::generator(4):
Here is the first partition of 4:
>> g()
[4]
Here is the second partition of 4:
>> g()
[3, 1]
And here are the remaining ones:
>> g(), g(), g(), g()
[2, 2], [2, 1, 1], [1, 1, 1, 1], FAIL
Generators become handy when you want to work with the 53174 partitions of 42:
>> g := partitions::generator(42):
g(), g(), g(), g(), g(), g()
[42], [41, 1], [40, 2], [40, 1, 1], [39, 3], [39, 2, 1]
Most of the sub-libraries of combinat
provide such generators.
Whenever possible (i.e. when it does not harm the computational complexity), we focus on providing the user with generic tools that cover many kinds of applications. For example, the libraries for partitions, integer vectors, and compositions share a very similar interface:
>> integerVectors::list(10, 3, MinPart = 2, MaxPart = 5, Inner = [2, 4, 2])
(Note: Inner = [2, 4, 2]
means that the three parts should
be respectively at least 2
, 4
and 2
).
[[4, 4, 2], [3, 5, 2], [3, 4, 3], [2, 5, 3], [2, 4, 4]]
>> compositions::list(5, MaxPart = 3, MinPart = 2, MinLength = 2, MaxLength = 3)
[[3, 2], [2, 3]]
>> partitions::list(5, MaxSlope = -1)
[[5], [4, 1], [3, 2]]
Those libraries actually use internally the same computational
engine combinat::integerListsLexTools
:
>> partitions::list(9, MinPart = 2, MaxPart = 5)
[[5, 4], [5, 2, 2], [4, 3, 2], [3, 3, 3], [3, 2, 2, 2]]
>> integerListsLexTools::list(9, 0, infinity, 2, 5, -infinity, 0)
[[5, 4], [5, 2, 2], [4, 3, 2], [3, 3, 3], [3, 2, 2, 2]]
In fact, combinat::integerListsLexTools
could also be used
to generate Motzkin and Dyck words, etc.
In the same spirit, instead of implementing a specific generator for standard tableaux, we implemented a generator for the linear extensions of a poset. We already reused this generator internally for generating standard binary search trees, and it could be reused as well for generating standard skew tableaux, standard ribbons, and so on (any volunteers?).
We also incorporated and extended the former CS
library by
S. Corteel, A. Denise, I. Dutour, and P. Zimmermann. This library
allows for manipulating any combinatorial classes that can be
defined by a deterministic grammar. Here we consider words of A's
and B's without two consecutive B's:
>> fiWords := decomposableObjects(
[FiWords = Union(Epsilon,
Atom(B),
Prod(FiWords, Atom(A)),
Prod(FiWords, Atom(A), Atom(B)))
]):
(Note: an
Epsilon
is an object of size 0 while an Atom
is
an object of size 1).
>> fiWords::list(4)
[Prod(Prod(Prod(B, A), A), A), Prod( Prod(Prod(Prod(Epsilon, A), A), A), A), Prod(Prod(Prod(Epsilon, A, B), A), A), Prod(Prod(B, A, B), A), Prod(Prod(Prod(Epsilon, A), A, B), A ), Prod(Prod(B, A), A, B), Prod(Prod(Prod(Epsilon, A), A), A, B), Prod(Prod(Epsilon, A, B), A, B)]
The result is not very readable, but this can be fixed by a quick substitution:
>> map(fiWords::list(4), p -> [eval(subs(p, Prod = id, Epsilon = null()))])
[[B, A, A, A], [A, A, A, A], [A, B, A, A], [B, A, B, A], [A, A, B, A], [B, A, A, B], [A, A, A, B], [A, B, A, B]]
Alternatively, we could have provided some rewriting rules in the grammar:
>> fiWords := decomposableObjects(
[FiWords = Alias(FiWordsRec, DOM_LIST),
FiWordsRec = Union(Epsilon(),
Atom(B),
Alias(Prod(FiWordsRec, Atom(A)), op),
Alias(Prod(FiWordsRec, Atom(A), Atom(B)), op))
]):
fiWords::list(4);
[[B, A, A, A], [A, A, A, A], [A, B, A, A], [B, A, B, A], [A, A, B, A], [B, A, A, B], [A, A, A, B], [A, B, A, B]]
This seems to work nicely. Let us count those words:
>> fiWords::count(i) $ i = 0..10
1, 2, 3, 5, 8, 13, 21, 34, 55, 89, 144
You will certainly recognize the Fibonnaci sequence. Not quite a surprise, the recurrence relation can be seen right away from the grammar. Actually, this recurrence relation is automatically determined by the library, and used for counting efficiently:
>> fiWords::recurrenceRelation() = 0
u(n - 1) - u(n) + u(n - 2) = 0
This also applies for several libraries which are based on
combinat::decomposableObjects
. For example, here is the
recurrence relation for binary trees:
>> collect(binaryTrees::grammar::recurrenceRelation(),
[u(n), u(n-1)], factor) = 0
u(n) (n + 1) - 2 u(n - 1) (2 n - 1) = 0
Let us define one of the most trivial combinatorial classes:
>> domain nonNegativeIntegers
// This is a domain (not a library):
inherits Dom::BaseDomain;
// This is a combinatorial class:
category Cat::CombinatorialClass;
// This is a facade domain:
axiom Ax::systemRep;
info := "Domain 'nonNegativeIntegers': the class of non negative integers";
// The domain of the elements of this class:
domtype := DOM_INT;
// Test if an object belongs to this class
isA := (obj) -> testtype(obj, Type::NonNegInt);
// The size of an integer is itself:
size := n -> n;
// There is exactly one integer of size n:
count := n -> 1;
list := n -> [n];
// No need to define generator; it is defined via list by default
end_domain:
>> testtype(x, nonNegativeIntegers), testtype(-3, nonNegativeIntegers),
testtype(3, nonNegativeIntegers);
FALSE, FALSE, TRUE
>> nonNegativeIntegers::count(4);
1
>> nonNegativeIntegers::list(4)
[4]
In a first approximation, the three lines inherits
,
category
, and axiom
may be safely ignored and kept
verbatim. For a deeper understanding, we strongly recommend to read
the detailed explanations about the implementation of combinatorial
classes in the design
notes.
Exercise: Improve the definition above so that
nonNegativeIntegers::count(-1)
raises an error instead of
returning -1.
It is often practical to define a sub-class of an existing class.
Here we show how to define the class of the permutations of
[1,2,3]
:
>> domain permutationsOf123
// Inherits all the methods from combinat::permutations
inherits combinat::permutations;
// This is a combinatorial class
category Cat::CombinatorialClass;
// This is a facade domain
axiom Ax::systemRep;
info := "Domain 'permutationsOf123': the class of the permutations of [1,2,3]";
// Redefinition of isA
isA := (p) -> permutations::isA(p, [1,2,3]);
// Redefinition of count
count := () -> permutations::count([1,2,3]);
// Redefinition of generator
generator := () -> permutations::generator([1,2,3]);
// No need to redefine list, since it is defined via generator by default
end_domain:
Let us use this new combinatorial class:
>> testtype(x, permutationsOf123),
testtype([1, 2, 3, 4], permutationsOf123),
testtype([1, 2, 2], permutationsOf123),
testtype([1, 3, 2], permutationsOf123);
FALSE, FALSE, FALSE, TRUE
>> permutationsOf123::count();
6
>> permutationsOf123::list()
[[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]]
Note: instead of implementing permutationsOf123 by hand, we could
have alternatively used the generic utility
combinat::subClass
; it allows to automatically define a
sub-class of an existing combinatorial class by providing extra
parameters to be passed down to all the methods count
,
list
, etc.:
>> permutationsOf123 := subClass(permutations, 3):
permutationsOf123::list()
[[1, 2, 3], [1, 3, 2], [2, 1, 3], [2, 3, 1], [3, 1, 2], [3, 2, 1]][2]
Exercises:
permutationsOf(m)
for the
permutations of the list m
(or of [1..m]
if m
is
an integer);
classOfSet({a,b,c})
that
creates the trivial combinatorial class composed of the elements
a
, b
, and c
. You may assume that all the
elements of the set are identifiers, and define the domtype
entry as DOM_IDENT
.
To conclude, we define the combinatorial class of Fibonacci
words. Essentially, we reuse the definition of fiWords
above,
and wrap it into a domain to add type checking:
>> domain FibonacciWords
inherits Dom::BaseDomain;
// The objects of this class are defined by a grammar
category Cat::DecomposableClass;
// This is a facade domain
axiom Ax::systemRep;
info := "Domain 'FibonacciWords': the class of Fibonacci words";
// The domain of the elements of this class
domtype := DOM_LIST;
// The type of the elements of this class
isA := proc(w : Type::AnyType, // a procedure that tests if w is a Fibonacci word
sz = 0-1 : Type::NonNegInt) // of optional size sz.
local i;
begin
if domtype(w) <>DOM_LIST then return(FALSE); end_if;
for i from 1 to nops(w) do
if (w[i]<>A and w[i]<>B) or
(i<nops(w) and w[i]=B and w[i+1]=B) then
return(FALSE);
end_if;
end_for;
if args(0) = 1 then TRUE else
bool(nops(w) = sz);
end_if;
end_proc;
// The size of a Fibonacci word is its length
size := nops;
// The grammar which defines the objects of this class
grammar := decomposableObjects(
[FiWords = Alias(FiWordsRec, DOM_LIST),
FiWordsRec = Union(Epsilon(),
Atom(B),
Alias(Prod(FiWordsRec, Atom(A)), op),
Alias(Prod(FiWordsRec, Atom(A), Atom(B)), op))
]);
end_domain:
Now, we can do type checking with this domain:
>> testtype(x, FibonacciWords),
testtype([A, B, C], FibonacciWords),
testtype([A, B, B], FibonacciWords),
testtype([A, B, A], FibonacciWords)
FALSE, FALSE, FALSE, TRUEAnd of course, we can still use all the previous functionalities of fiWords:
>> FibonacciWords::list(4);
[[B, A, A, A], [A, A, A, A], [A, B, A, A], [B, A, B, A], [A, A, B, A], [B, A, A, B], [A, A, A, B], [A, B, A, B]]
>> FibonacciWords::count(4);
8
>> FibonacciWords::unrank(2, 4)
[A, A, A, A]
Exercises:
toString
and fromString
which
converts a Fibonacci word to and from a string.
FibonacciWords(A,B)
for the
Fibonacci words on the letter A
and B
.
We now demonstrate how to do sample computations with predefined combinatorial algebras, starting with the algebra of symmetric functions. Note that we really consider those predefined algebras as mere examples of use of this package; important and useful examples of course, but just examples.
We define the ring of symmetric functions over the rational numbers:
>> S := examples::SymmetricFunctions(Dom::Rational);
examples::SymmetricFunctions(Dom::Rational)
This ring has several remarkable families like the symmetric power-sums pk: recall that the symmetric power-sum pk expands on any given specified alphabet (i.e. set of variables) as the sum of all the variables elevated to the power k; furthermore, given a partition λ:=(λ1,...,λk), the product of pλ1...pλk is denoted by pλ:
>> p1 := S::p([1]);
p1([x, y, z]);
p[1] x + y + z
>> p2 := S::p([2]);
p2([x, y, z])
p[2] 2 2 2 x + y + z
>> p421 := S::p([4, 2, 1]);
p421([x,y,z])
p[4, 2, 1] 2 2 2 4 4 4 (x + y + z ) (x + y + z ) (x + y + z)Note: the product being commutative, the order in which the terms appear in the expansion above depends on MuPAD internal order, and is mathematically irrelevant.
Actually, the ring of symmetric functions is the free commutative algebra on the symmetric power-sums:
>> p2 * p1 * p2 * p421
p[4, 2, 2, 2, 1, 1]
Note that any expression is immediately expanded by the system:
>> (p421 + 3*p2)*(1/4*p1 - p2)
1/4 p[4, 2, 1, 1] - 3 p[2, 2] - p[4, 2, 2, 1] + 3/4 p[2, 1]
This happens because the call S::p([1])
returns a typed
object, for which the standard arithmetic operators are overloaded:
>> domtype(p1);
examples::SymmetricFunctionsTools::powersum(Dom::Rational)
>> S::p
examples::SymmetricFunctionsTools::powersum(Dom::Rational)
That is, S::powersum
(or S::p
for short) really
represents the domain of symmetric functions expanded on the
power-sums basis. If at some time you do not want the expansion to
take place, the objects can always be converted to expressions:
>> f := (expr(p421) + 3*expr(p2))*(1/4*expr(p1) - expr(p2))
/ p[1] \ (3 p[2] + p[4, 2, 1]) | ---- - p[2] | \ 4 /
Here is how to convert efficiently the expression back into an
object of the domain S::p
:
>> f := eval(subs(f, p = S::p::domainWrapper))
1/4 p[4, 2, 1, 1] - 3 p[2, 2] - p[4, 2, 2, 1] + 3/4 p[2, 1]
This requires the use of a small trick because of the indexed
notation for the basis elements. S::p::domainWrapper is a special
MuPAD object which, when used as S::p::domainWrapper[3,2]
,
returns a call to S::p([3,2]).
Of course, examples::SymmetricFunctions
provides the other
classical bases of symmetric functions, like the elementary
symmetric functions S::e
, the monomial symmetric functions
S::m
, the homogeneous symmetric functions S::h
, the
Schur functions S::s
, etc.:
>> expand(S::e([2])([x,y,z]))
x y + x z + y z
>> expand(S::m([2, 1])([x,y,z]))
2 2 2 2 2 2 x y + x y + x z + x z + y z + y z
>> expand(S::h([2])([x,y,z]))
2 2 2 x y + x z + y z + x + y + z
>> expand(S::s([2])([x,y,z]))
2 2 2 x y + x z + y z + x + y + z
Here is how to convert from one basis to the other
:>> f := S::p([4]);
S::e(f);
S::h(f);
S::s(f);
S::m(f)
p[4] e[1, 1, 1, 1] - 4 e[2, 1, 1] + 4 e[3, 1] - 4 e[4] + 2 e[2, 2] - h[1, 1, 1, 1] + 4 h[2, 1, 1] - 4 h[3, 1] + 4 h[4] - 2 h[2, 2] - s[1, 1, 1, 1] + s[2, 1, 1] - s[3, 1] + s[4] m[4]
When multiplying two symmetric functions which are not expressed in the same basis, the system will make an implicit conversion, and return the result in one or the other of the two bases:
>> S::m([2]) * S::s([2])
- s[2, 1, 1] + s[4] + s[2, 2]
If you want to force the product to be done on a given basis, you can call the proper conversion explicitly:
>> S::m(S::m([2])) * S::s([2]);
m[2, 1, 1] + m[3, 1] + m[4] + 2 m[2, 2]
Now, we can combine everything, and do some complicated calculation:
>> S::p( S::m([1]) * ( S::e([3])*S::s([2]) + S::s([3]) ) )
1/12 p[1, 1, 1, 1, 1, 1] + 1/6 p[3, 2, 1] - 1/6 p[2, 1, 1, 1, 1] + 1/6 p[3, 1, 1, 1] - 1/4 p[2, 2, 1, 1] + 1/6 p[1, 1, 1, 1] + 1/2 p[2, 1, 1] + 1/3 p[3, 1]
Finally, there is some basic support for the Hall-Littlewood functions, in the P and Q' basis, which we demonstrate now. We need to take some ground field which contains the parameter t of those functions. The simplest (and actually most efficient with the current MuPAD version), is to take the full field of expressions as coefficient ring:
>> S := examples::SymmetricFunctions(Dom::ExpressionField()):
Here is the Hall-Littlewood function Q'(3,2,1,1):
>> el := S::QP([3, 2, 1, 1])
QP[3, 2, 1, 1]
The expansion of el
in terms of Schur functions reads as:
>> S::s(el)
2 t s[3, 2, 2] + (t + t ) s[3, 3, 1] + t s[4, 1, 1, 1] + 2 3 4 2 3 4 (t + t + t ) s[4, 3] + (t + t + t ) s[5, 1, 1] + 3 4 5 4 5 6 (2 t + t + t ) s[5, 2] + (t + t + t ) s[6, 1] + 7 2 3 t s[7] + s[3, 2, 1, 1] + (t + 2 t + t ) s[4, 2, 1]
The expansion of el
on the alphabet (q,qt) reads as:
>> expand(el([q, q*t]))
7 5 7 6 7 7 7 8 7 9 7 10 4 q t + 7 q t + 10 q t + 9 q t + 6 q t + 5 q t + 7 11 7 12 7 13 7 14 3 q t + 2 q t + q t + q t
We now turn to the central feature of the MuPAD-Combinat package: the ability to easily implement new combinatorial algebras. We start by the free associative algebra over the rational numbers generated by non commutative letters a,b,c,d,.... Its basis is indexed by words, and the product of two basis elements is obtained by concatenating the corresponding words:
>> domain FreeAlgebra // line 1
inherits Dom::FreeModule(Dom::Rational, combinat::words); // line 2
category Cat::AlgebraWithBasis(Dom::Rational); // line 3
one := dom::term([]); // line 5
mult2Basis := dom::term @ _concat; // line 6
end_domain: // line 7
We will explain the bits of this definition in a minute after a few examples of use. Let us define two elements of the free algebra:
>> x := FreeAlgebra([a, b, c]);
y := FreeAlgebra([d, e])
B([a, b, c]) B([d, e])The
B
just stands for the name of the basis. We can compute
linear combinations and products of x
and y
:
>> 3 * x;
x + y;
x * y
3 B([a, b, c]) B([d, e]) + B([a, b, c]) B([a, b, c, d, e])Here is a more complicated expression:
>> x * (2*x + y) + (3 + y/2)^2
1/4 B([d, e, d, e]) + 2 B([a, b, c, a, b, c]) + B([a, b, c, d, e]) + 3 B([d, e]) + 9 B([])
Note how the 3
in the expression is automatically converted
into an element of the domain; declaring that FreeAlgebra
was
an algebra (with a unit) automatically defined the natural embedding
of the coefficient ring into it.
We turn to the explanation of the implementation of FreeAlgebra
above. Line 1 states that we are declaring a new domain called
FreeAlgebra
(a new class in the usual object oriented
terminology). Line 2 lets FreeAlgebra
inherit its
implementation from the free module over the rationals
(Dom::Rational
) with basis indexed by words
(combinat::words
). Line 3 states that FreeAlgebra
is
actually an algebra with a distinguished basis; this allows, in
particular, to define the multiplication by linearity on the basis.
Line 5 defines that the unit of FreeAlgebra
is the empty word
(dom
refers to the domain being defined, and dom::term
is
a constructor that takes an element of the basis, and returns it as an
element of the domain). Finally, line 6 states that two elements of
the basis are multiplied by concatenating them, and making an element
of the domain with the result (@
denotes the composition of
functions). That's it.
Let us define the free commutative algebra on the letters a,b,c,...:
>> domain FreeCommutativeAlgebra
inherits Dom::FreeModule(Dom::Rational, combinat::words);
category Cat::AlgebraWithBasis(Dom::Rational);
one := dom::term([]);
straightenBasis := dom::term @ sort;
mult2Basis := dom::straightenBasis @ _concat;
end_domain:
Note that we cheated a little bit: we declared that the basis of
FreeCommutativeAlgebra
consisted of words, whereas it really
consists of words up to permutation of its letters: B([a,b])
and B([b,a])
represent the same element of the algebra. A
careful implementation should define the combinatorial class of
words up to permutation, and use it as the basis of
FreeCommutativeAlgebra
.
To enforce the uniqueness of the representation, we straighten the
words in the basis by sorting them. This is the job of the
straightenBasis
constructor.
>> x := FreeCommutativeAlgebra([a, b]);
y := FreeCommutativeAlgebra([c, b, a])
B([a, b]) B([a, b, c])The product of two words is then defined by concatenating them and straightening the result:
>> x * y;
y * x
B([a, a, b, b, c]) B([a, a, b, b, c])If efficiency was at a premium, we could have used the MuPAD function
listlib::merge
which merges sorted lists instead.
Note that two elements of FreeCommutativeAlgebra
and of
FreeAlgebra
may happen to be printed out the same way:
>> x := FreeAlgebra([a]);
y := FreeCommutativeAlgebra([a])
B([a]) B([a])and even share the exact same internal representation:
>> bool(extop(x) = extop(y))
TRUEHowever, they are not equal, because they are not in the same domain:
>> bool(x = y);
domtype(x), domtype(y)
FALSE FreeAlgebra, FreeCommutativeAlgebraSo, even if they share the same name of basis, there is no risk of confusion; for example we are not allowed to multiply them together:
>> x * y
Error: Don't know how to multiply a FreeAlgebra by a FreeCommu\ tativeAlgebraOf course, this is still confusing for the user. He or she may always customize the basis names (as many other things) at any time should he or she wish to do so:
>> FreeAlgebra::basisName := hold(T):
FreeCommutativeAlgebra::basisName := hold(S):
x, y
T([a]), S([a])Here,
T
stands for ``Tensor algebra'', while S
stands for
``Symmetric algebra''. The hold
is there for safety, to avoid
trouble if one of the identifiers T
or S
is assigned a
value.
We can define the natural evaluation morphism from FreeAlgebra
to FreeCommutativeAlgebra
by linearity on the words; a word
itself is simply sorted, and converted into an element of
FreeCommutativeAlgebra
:
>> evaluation := operators::makeLinear(FreeCommutativeAlgebra::term @ sort,
Source = FreeAlgebra,
ImageSet = FreeCommutativeAlgebra):
Let us apply this morphism to the sum of two words which only differ by a permutation:
>> x := FreeAlgebra([c, b, a]) + FreeAlgebra([c, a, b]);
T([c, a, b]) + T([c, b, a])
>> evaluation(x);
2 S([a, b, c])The evaluation morphism is actually quite canonical, so it can make sense to declare it as a conversion to the system. This can be achieved with the
operators::overloaded::declareConversion
function:
>> operators::overloaded::declareConversion(FreeAlgebra,
FreeCommutativeAlgebra,
evaluation):
FreeCommutativeAlgebra(x)
2 S([a, b, c])Here, the conversion has been declared as implicit. If an expression mixes elements of
FreeAlgebra
and FreeCommutativeAlgebra
,
the former are automatically converted into
FreeCommutativeAlgebra
:
>> FreeCommutativeAlgebra([a, b]) + FreeAlgebra([c,b,a])
S([a, b]) + S([a, b, c])Of course, such a feature is questionable. Depending on the context, it can prove very practical, or on the contrary dangerous. The user is the only judge, and she or he can restrict the scope of this conversion by using the
Explicit
option. In this case, the
conversion will only be applied if requested explicitly by
convert
or by new
:
>> operators::overloaded::declareConversion(FreeAlgebra, FreeCommutativeAlgebra,
evaluation, Explicit):
FreeCommutativeAlgebra(x);
2 S([a, b, c])
>> FreeCommutativeAlgebra([a, b]) + FreeAlgebra([c,b,a])
Error: Don't know how to add a FreeCommutativeAlgebra and a Fr\ eeAlgebraTypically, for symmetric functions, we only provided explicit conversions to construct symmetric functions from partitions because those conversions are not canonical at all: the Schur function
s[3,2,1]
obtained by converting the partition [3,2,1]
has
nothing to do with the elementary function e[3,2,1]
. We refer to
the design notes and to the
documentation of the operators::overloaded
library for details
on the mechanism we use for defining automatic conversions and
overloaded operators and functions. Note that it is not (yet)
completely possible to declare new conversions as above when the
target domain of the conversion is one of the predefined domains of
the MuPAD library.
To continue our exploration, we implement variations on the two
previous domains, where we assume that the algebra generators are
indexed by 1,2,...
. The basis elements of the free algebra
and of the free commutative algebra are now respectively indexed by
compositions and partitions.
>> domain FreeAlgebraInteger
inherits Dom::FreeModule(Dom::Rational, combinat::compositions);
category Cat::AlgebraWithBasis(Dom::Rational);
basisName := hold(E);
exprTerm := dom::exprTermIndex;
one := dom::term([]);
mult2Basis := dom::term @ _concat;
end_domain:
domain FreeCommutativeAlgebraInteger
inherits Dom::FreeModule(Dom::Rational, combinat::partitions);
category Cat::AlgebraWithBasis(Dom::Rational);
basisName := hold(e);
exprTerm := dom::exprTermIndex;
one := dom::term([]);
straightenBasis := dom::term @ revert @ sort;
mult2Basis := dom::straightenBasis @ _concat;
end_domain:
The reader may have recognized here respectively the commutative and non commutative symmetric functions, expanded on the elementary symmetric functions; hence the basis names. To shorten the notations, we define two aliases, and declare the same evaluation conversion as before:
>> alias(NCSF = FreeAlgebraInteger,
SF = FreeCommutativeAlgebraInteger):
operators::overloaded::declareConversion(NCSF, SF,
operators::makeLinear(SF::straightenBasis,
Source = NCSF,
ImageSet = SF)):
>> x := NCSF([1, 3, 2]);
E[1, 3, 2]
>> y := SF ([1, 3, 2])
e[3, 2, 1]
>> SF(x)
e[3, 2, 1]
>> bool(SF(x)=y)
TRUE
Let us analyze the differences with our previous implementation of
the free algebras. First, we chose an indexed notation for the
basis elements, as this notation is more compact and quite
conventional in other systems. This is the purpose of the line
exprTerm := dom::exprTermIndex
: the method exprTerm
of the domain is called to convert a term into an
expression, as well as to print a term if there is no
printTerm
method; dom::exprTermIndex
is a possible
implementation of exprTerm
, inherited from the category,
which gives indexed notations. The other difference is that,
following the usual convention, the integers in the partitions are
sorted decreasingly. Here, this is suboptimally achieved by
reverting the list after sorting it in the
SF::straightenBasis
method.
A disadvantage of this implementation of SF
is that elements
with many repetitions are not represented compactly:
>> SF([1])^10
e[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]One might prefer to use another basis for SF, where the element above would be represented as the first generator to the power 10. This can be done via the usual exponent notation for partitions. The basis of the algebra now consists of integer vectors. The product of two elements is simply obtained by adding up the vectors part by part, which can conveniently be implemented using the
zip
MuPAD
function.
>> domain SFExp
inherits Dom::FreeModule(Dom::Rational, combinat::integerVectors);
category Cat::AlgebraWithBasis(Dom::Rational);
basisName := hold(e);
exprTerm := dom::exprTermIndex;
one := dom::term([]);
mult2Basis := (v1,v2) -> dom::term(zip(v1,v2,_plus,0));
end_domain:
Let us do some computations:
>> SFExp([1]);
SFExp([2, 0, 1])*SFExp([1, 1]);
SFExp([1])^10
e[1] e[3, 1, 1] e[10]This notation could be confusing, so we override it:
>> SFExp::exprTerm := v -> _mult(dom::basisName.i^v[i] $ i=1..nops(v)):
SFExp([1]);
SFExp([2, 0, 1])*SFExp([1, 1]);
SFExp([1])^10
e1 3 e1 e2 e3 10 e1As is, the elements of this algebra are not uniquely represented. For example, the first generator of the algebra can be represented by any of
[1], [1,0], [1,0,0], ...
:
>> SFExp([1]), SFExp([1, 0]), SFExp([1, 0, 0]);
bool(SFExp([1]) = SFExp([1, 0]))
e1, e1, e1 FALSE
We leave it up as an exercise for the reader to fix this bug by
implementing a straightenBasis
method which strips out the
trailing zeroes of the basis elements.
Of course, SF
and SFExp
really represent the same
algebra; only the internal data representation changes. So, we
provide as conversions the reciprocal isomorphisms obtained by
extending by linearity the bijections
combinat::partitions::fromExp
and
combinat::partitions::toExp
:
>> operators::overloaded::declareConversion(SFExp, SF,
operators::makeLinear(SF::term @ combinat::partitions::fromExp,
Source = SFExp,
ImageSet = SF)):
operators::overloaded::declareConversion(SF, SFExp,
operators::makeLinear(SFExp::term @ combinat::partitions::toExp,
Source = SF,
ImageSet = SFExp)):
Here is a simple conversion:
>> SF([4, 3, 3, 1]), SFExp(SF([4, 3, 3, 1]))
2 e[4, 3, 3, 1], e1 e3 e4Let us check on an example that the conversions are indeed morphisms:
>> x := SF([3, 1]):
y := SF([4, 3, 2]):
x * y, SF( SFExp(x) * SFExp(y) )
e[4, 3, 3, 2, 1], e[4, 3, 3, 2, 1]We can write expressions that mix elements of both domains, and let the system find a way to convert them appropriately:
>> ( 1 + SF([3, 1])*x ) * SFExp([2, 0]) + SFExp([1])
4 2 2 e1 e3 + e1 + e1A priori, the representation of the result cannot be predicted; it depends on how the overloading mechanisms choose to resolve the conversions. If the user prefers one of the representations, she or he can take over the control at any level of the expression by forcing proper conversions:
>> (1 + SF([3, 1])*x) * SF( SFExp([2, 0]) ) + SF( SFExp([1]) ) ;
SF( (1 + SF([3, 1])*x) * SFExp([2, 0]) ) + SF( SFExp([1]) ) ;
SF( (1 + SF([3, 1])*x) * SFExp([2, 0]) + SFExp([1]) ) ;
e[3, 3, 1, 1, 1, 1] + e[1, 1] + e[1] e[3, 3, 1, 1, 1, 1] + e[1, 1] + e[1] e[3, 3, 1, 1, 1, 1] + e[1, 1] + e[1]
The implicit conversions are automatically applied transitively. As a
consequence, we have without further work a conversion from NCSF
to SF
:
>> SFExp(NCSF([1, 4, 2, 2]))
2 e1 e2 e4Another consequence is that, when there are n different representations for an algebra (say the symmetric functions expressed on any of the
p
, e
, m
, h
, or s
basis),
it is enough to implement 2n conversions to be able to get all the
n(n-1) possible conversions. Of course, it is still possible to
implement some extra direct conversions for improved efficiency; when
there are several ways to convert an element from one domain to
another, the system always uses one of the shortest ones.
We conclude by a typical computation which involves basic linear algebra: we find the minimal polynomial of an element of the group algebra of the symmetric group of order 4:
>> SymGroup4 := subClass(permutations, 4):
domain AlgSymGroup4
inherits Dom::FreeModule(Dom::ExpressionField(), SymGroup4);
category Cat::AlgebraWithBasis(Dom::ExpressionField());
basisName := hold(p);
exprTerm := dom::exprTermIndex;
one := dom::term([1,2,3,4]);
mult2Basis := dom::term @ permutations::mult2;
end_domain:
x := AlgSymGroup4([2,3,4,1]):
y := x^0 + a1*x^1 + a2*x^2 + a3*x^3 + a4*x^4;
a3 p[4, 1, 2, 3] + a2 p[3, 4, 1, 2] + a1 p[2, 3, 4, 1] + (a4 + 1) p[1, 2, 3, 4]
>> solve([coeff(y)], [a1,a2,a3,a4]);
{[a1 = 0, a2 = 0, a3 = 0, a4 = -1]}
>> subs(z^0 + a1*z^1 + a2*z^2 + a3*z^3 + a4*z^4, op(last(1),1))
4 1 - z
We conclude this guided tour by a summary of the current features. A first part of the package consists of predefined libraries to count, enumerate, and manipulate standard combinatorial objects (partitions, compositions, sets, words, permutations, tableaux, trees, ...), together with generic tools to help define new combinatorial classes:
A second part consists of tools to build combinatorial algebras. The typical usage is to take a vector space with basis indexed by some combinatorial objects, to define a product for two basis elements and to extend the product by bilinearity. The system automatically takes care of the parsing of algebraic elements, of extensions of functions by linearity, bi-linearity or associativity, of conversions between different bases, etc. Similarly, one may define coproducts, antipods, and similar operators, to implement stronger structures such as Hopf algebras. Some preliminary work has been done to manipulate Lie algebras as well. In short, the system takes care of the algebraic work, so that one can concentrate on the underlying combinatorics rather than on the programming.
As examples of usage and applications, we provide a library for the algebra of symmetric functions, and we have (currently undocumented) libraries for the algebra of non commutative symmetric functions, the algebra of (free) quasi-symmetric functions, the Loday-Ronco algebra, the Weyl algebra, the rational Steenrod algebra, the type A Hecke and Hecke-Clifford algebras, as well as invariant rings of permutation groups, and group algebras.
In the mid-term (or short term if you help us!)we plan to provide predefined libraries for the free symmetric algebra, the algebra of matrix quasi symmetric functions, the descents and peaks algebras, general Ore-algebras, the symmetric Weyl algebra, the algebra of multi-symmetric functions, the divided power algebra, free Lie algebras, etc.
Finally, we provide a library for manipulating weighted automatons, and related semi-rings.
MuPAD Combinat, an open source algebraic combinatorics package