[Previous] [Next] [Contents]

A guided tour through MuPAD-Combinat

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.

Two examples of combinatorial algebras

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    /
            /  \
           / \
        

MuPAD-Combinat, step by step

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):
      

Using predefined combinatorial functions and classes

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

Defining new combinatorial classes

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:

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, TRUE
And 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:

Using predefined combinatorial algebras

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

Defining new combinatorial algebras

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))
                                   TRUE
However, they are not equal, because they are not in the same domain:
>> bool(x = y);
   domtype(x), domtype(y)
                                   FALSE
      
                    FreeAlgebra, FreeCommutativeAlgebra
So, 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\
      tativeAlgebra
Of 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\
      eeAlgebra
Typically, 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
                                   e1
As 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  e4
Let 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  + e1
A 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  e4
Another 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.

A typical computation

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

Current features

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:

Most predefined libraries actually make use of these engines.

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.

[Previous] [Next] [Contents]


MuPAD Combinat, an open source algebraic combinatorics package