4. Basic Concepts

This chapter will describe the different fundamental objects that can be handled by GiNaC. But before doing so, it is worthwhile introducing you to the more commonly used class of expressions, representing a flexible meta-class for storing all mathematical objects.

4.1 Expressions

The most common class of objects a user deals with is the expression ex, representing a mathematical object like a variable, number, function, sum, product, etc... Expressions may be put together to form new expressions, passed as arguments to functions, and so on. Here is a little collection of valid expressions:

ex MyEx1 = 5;                       // simple number
ex MyEx2 = x + 2*y;                 // polynomial in x and y
ex MyEx3 = (x + 1)/(x - 1);         // rational expression
ex MyEx4 = sin(x + 2*y) + 3*z + 41; // containing a function
ex MyEx5 = MyEx4 + 1;               // similar to above

Expressions are handles to other more fundamental objects, that often contain other expressions thus creating a tree of expressions (See section A. Internal Structures, for particular examples). Most methods on ex therefore run top-down through such an expression tree. For example, the method has() scans recursively for occurrences of something inside an expression. Thus, if you have declared MyEx4 as in the example above MyEx4.has(y) will find y inside the argument of sin and hence return true.

The next sections will outline the general picture of GiNaC's class hierarchy and describe the classes of objects that are handled by ex.

4.1.1 Note: Expressions and STL containers

GiNaC expressions (ex objects) have value semantics (they can be assigned, reassigned and copied like integral types) but the operator < doesn't provide a well-defined ordering on them. In STL-speak, expressions are `Assignable' but not `LessThanComparable'.

This implies that in order to use expressions in sorted containers such as std::map<> and std::set<> you have to supply a suitable comparison predicate. GiNaC provides such a predicate, called ex_is_less. For example, a set of expressions should be defined as std::set<ex, ex_is_less>.

Unsorted containers such as std::vector<> and std::list<> don't pose a problem. A std::vector<ex> works as expected.

See section 5.1 Getting information about expressions, for more about comparing and ordering expressions.

4.2 Automatic evaluation and canonicalization of expressions

GiNaC performs some automatic transformations on expressions, to simplify them and put them into a canonical form. Some examples:

ex MyEx1 = 2*x - 1 + x;  // 3*x-1
ex MyEx2 = x - x;        // 0
ex MyEx3 = cos(2*Pi);    // 1
ex MyEx4 = x*y/x;        // y

This behavior is usually referred to as automatic or anonymous evaluation. GiNaC only performs transformations that are

• at most of complexity
• algebraically correct, possibly except for a set of measure zero (e.g. is transformed to although this is incorrect for )

There are two types of automatic transformations in GiNaC that may not behave in an entirely obvious way at first glance:

• The terms of sums and products (and some other things like the arguments of symmetric functions, the indices of symmetric tensors etc.) are re-ordered into a canonical form that is deterministic, but not lexicographical or in any other way easy to guess (it almost always depends on the number and order of the symbols you define). However, constructing the same expression twice, either implicitly or explicitly, will always result in the same canonical form.
• Expressions of the form 'number times sum' are automatically expanded (this has to do with GiNaC's internal representation of sums and products). For example

  ex MyEx5 = 2*(x + y); // 2*x+2*y ex MyEx6 = z*(x + y); // z*(x+y)

The general rule is that when you construct expressions, GiNaC automatically creates them in canonical form, which might differ from the form you typed in your program. This may create some awkward looking output (`-y+x' instead of `x-y') but allows for more efficient operation and usually yields some immediate simplifications.

Internally, the anonymous evaluator in GiNaC is implemented by the methods

ex ex::eval(int level = 0) const;
ex basic::eval(int level = 0) const;

but unless you are extending GiNaC with your own classes or functions, there should never be any reason to call them explicitly. All GiNaC methods that transform expressions, like subs() or normal(), automatically re-evaluate their results.

4.3 Error handling

GiNaC reports run-time errors by throwing C++ exceptions. All exceptions generated by GiNaC are subclassed from the standard exception class defined in the `<stdexcept>' header. In addition to the predefined logic_error, domain_error, out_of_range, invalid_argument, runtime_error, range_error and overflow_error types, GiNaC also defines a pole_error exception that gets thrown when trying to evaluate a mathematical function at a singularity.

The pole_error class has a member function

int pole_error::degree() const;

that returns the order of the singularity (or 0 when the pole is logarithmic or the order is undefined).

When using GiNaC it is useful to arrange for exceptions to be caught in the main program even if you don't want to do any special error handling. Otherwise whenever an error occurs in GiNaC, it will be delegated to the default exception handler of your C++ compiler's run-time system which usually only aborts the program without giving any information what went wrong.

Here is an example for a main() function that catches and prints exceptions generated by GiNaC:

#include <iostream>
#include <stdexcept>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

int main()
{
    try {
        ...
        // code using GiNaC
        ...
    } catch (exception &p) {
        cerr << p.what() << endl;
        return 1;
    }
    return 0;
}

4.4 The Class Hierarchy

GiNaC's class hierarchy consists of several classes representing mathematical objects, all of which (except for ex and some helpers) are internally derived from one abstract base class called basic. You do not have to deal with objects of class basic, instead you'll be dealing with symbols, numbers, containers of expressions and so on.

To get an idea about what kinds of symbolic composites may be built we have a look at the most important classes in the class hierarchy and some of the relations among the classes:

The abstract classes shown here (the ones without drop-shadow) are of no interest for the user. They are used internally in order to avoid code duplication if two or more classes derived from them share certain features. An example is expairseq, a container for a sequence of pairs each consisting of one expression and a number (numeric). What is visible to the user are the derived classes add and mul, representing sums and products. See section A. Internal Structures, where these two classes are described in more detail. The following table shortly summarizes what kinds of mathematical objects are stored in the different classes:

4.5 Symbols

Symbolic indeterminates, or symbols for short, are for symbolic manipulation what atoms are for chemistry.

A typical symbol definition looks like this:

symbol x("x");

This definition actually contains three very different things:

• a C++ variable named x
• a symbol object stored in this C++ variable; this object represents the symbol in a GiNaC expression
• the string "x" which is the name of the symbol, used (almost) exclusively for printing expressions holding the symbol

Symbols have an explicit name, supplied as a string during construction, because in C++, variable names can't be used as values, and the C++ compiler throws them away during compilation.

It is possible to omit the symbol name in the definition:

symbol x;

In this case, GiNaC will assign the symbol an internal, unique name of the form symbolNNN. This won't affect the usability of the symbol but the output of your calculations will become more readable if you give your symbols sensible names (for intermediate expressions that are only used internally such anonymous symbols can be quite useful, however).

Now, here is one important property of GiNaC that differentiates it from other computer algebra programs you may have used: GiNaC does not use the names of symbols to tell them apart, but a (hidden) serial number that is unique for each newly created symbol object. In you want to use one and the same symbol in different places in your program, you must only create one symbol object and pass that around. If you create another symbol, even if it has the same name, GiNaC will treat it as a different indeterminate.

Observe:

ex f(int n)
{
    symbol x("x");
    return pow(x, n);
}

int main()
{
    symbol x("x");
    ex e = f(6);

    cout << e << endl;
     // prints "x^6" which looks right, but...

    cout << e.degree(x) << endl;
     // ...this doesn't work. The symbol "x" here is different from the one
     // in f() and in the expression returned by f(). Consequently, it
     // prints "0".
}

One possibility to ensure that f() and main() use the same symbol is to pass the symbol as an argument to f():

ex f(int n, const ex & x)
{
    return pow(x, n);
}

int main()
{
    symbol x("x");

    // Now, f() uses the same symbol.
    ex e = f(6, x);

    cout << e.degree(x) << endl;
     // prints "6", as expected
}

Another possibility would be to define a global symbol x that is used by both f() and main(). If you are using global symbols and multiple compilation units you must take special care, however. Suppose that you have a header file `globals.h' in your program that defines a symbol x("x");. In this case, every unit that includes `globals.h' would also get its own definition of x (because header files are just inlined into the source code by the C++ preprocessor), and hence you would again end up with multiple equally-named, but different, symbols. Instead, the `globals.h' header should only contain a declaration like extern symbol x;, with the definition of x moved into a C++ source file such as `globals.cpp'.

A different approach to ensuring that symbols used in different parts of your program are identical is to create them with a factory function like this one:

const symbol & get_symbol(const string & s)
{
    static map<string, symbol> directory;
    map<string, symbol>::iterator i = directory.find(s);
    if (i != directory.end())
        return i->second;
    else
        return directory.insert(make_pair(s, symbol(s))).first->second;
}

This function returns one newly constructed symbol for each name that is passed in, and it returns the same symbol when called multiple times with the same name. Using this symbol factory, we can rewrite our example like this:

ex f(int n)
{
    return pow(get_symbol("x"), n);
}

int main()
{
    ex e = f(6);

    // Both calls of get_symbol("x") yield the same symbol.
    cout << e.degree(get_symbol("x")) << endl;
     // prints "6"
}

Instead of creating symbols from strings we could also have get_symbol() take, for example, an integer number as its argument. In this case, we would probably want to give the generated symbols names that include this number, which can be accomplished with the help of an ostringstream.

In general, if you're getting weird results from GiNaC such as an expression `x-x' that is not simplified to zero, you should check your symbol definitions.

As we said, the names of symbols primarily serve for purposes of expression output. But there are actually two instances where GiNaC uses the names for identifying symbols: When constructing an expression from a string, and when recreating an expression from an archive (see section 5.15 Input and output of expressions).

In addition to its name, a symbol may contain a special string that is used in LaTeX output:

symbol x("x", "\\Box");

This creates a symbol that is printed as "x" in normal output, but as "\Box" in LaTeX code (See section 5.15 Input and output of expressions, for more information about the different output formats of expressions in GiNaC). GiNaC automatically creates proper LaTeX code for symbols having names of greek letters (`alpha', `mu', etc.).

Symbols in GiNaC can't be assigned values. If you need to store results of calculations and give them a name, use C++ variables of type ex. If you want to replace a symbol in an expression with something else, you can invoke the expression's .subs() method (see section 5.3 Substituting expressions).

By default, symbols are expected to stand in for complex values, i.e. they live in the complex domain. As a consequence, operations like complex conjugation, for example (see section 5.13 Complex Conjugation), do not evaluate if applied to such symbols. Likewise log(exp(x)) does not evaluate to x, because of the unknown imaginary part of x. On the other hand, if you are sure that your symbols will hold only real values, you would like to have such functions evaluated. Therefore GiNaC allows you to specify the domain of the symbol. Instead of symbol x("x"); you can write realsymbol x("x"); to tell GiNaC that x stands in for real values.

4.6 Numbers

For storing numerical things, GiNaC uses Bruno Haible's library CLN. The classes therein serve as foundation classes for GiNaC. CLN stands for Class Library for Numbers or alternatively for Common Lisp Numbers. In order to find out more about CLN's internals, the reader is referred to the documentation of that library. See Info file `cln', node `Introduction', for more information. Suffice to say that it is by itself build on top of another library, the GNU Multiple Precision library GMP, which is an extremely fast library for arbitrary long integers and rationals as well as arbitrary precision floating point numbers. It is very commonly used by several popular cryptographic applications. CLN extends GMP by several useful things: First, it introduces the complex number field over either reals (i.e. floating point numbers with arbitrary precision) or rationals. Second, it automatically converts rationals to integers if the denominator is unity and complex numbers to real numbers if the imaginary part vanishes and also correctly treats algebraic functions. Third it provides good implementations of state-of-the-art algorithms for all trigonometric and hyperbolic functions as well as for calculation of some useful constants.

The user can construct an object of class numeric in several ways. The following example shows the four most important constructors. It uses construction from C-integer, construction of fractions from two integers, construction from C-float and construction from a string:

#include <iostream>
#include <ginac/ginac.h>
using namespace GiNaC;

int main()
{
    numeric two = 2;                      // exact integer 2
    numeric r(2,3);                       // exact fraction 2/3
    numeric e(2.71828);                   // floating point number
    numeric p = "3.14159265358979323846"; // constructor from string
    // Trott's constant in scientific notation:
    numeric trott("1.0841015122311136151E-2");
    
    std::cout << two*p << std::endl;  // floating point 6.283...
    ...

The imaginary unit in GiNaC is a predefined numeric object with the name I:

    ...
    numeric z1 = 2-3*I;                    // exact complex number 2-3i
    numeric z2 = 5.9+1.6*I;                // complex floating point number
}

It may be tempting to construct fractions by writing numeric r(3/2). This would, however, call C's built-in operator / for integers first and result in a numeric holding a plain integer 1. Never use the operator / on integers unless you know exactly what you are doing! Use the constructor from two integers instead, as shown in the example above. Writing numeric(1)/2 may look funny but works also.

We have seen now the distinction between exact numbers and floating point numbers. Clearly, the user should never have to worry about dynamically created exact numbers, since their `exactness' always determines how they ought to be handled, i.e. how `long' they are. The situation is different for floating point numbers. Their accuracy is controlled by one global variable, called Digits. (For those readers who know about Maple: it behaves very much like Maple's Digits). All objects of class numeric that are constructed from then on will be stored with a precision matching that number of decimal digits:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

void foo()
{
    numeric three(3.0), one(1.0);
    numeric x = one/three;

    cout << "in " << Digits << " digits:" << endl;
    cout << x << endl;
    cout << Pi.evalf() << endl;
}

int main()
{
    foo();
    Digits = 60;
    foo();
    return 0;
}

The above example prints the following output to screen:

in 17 digits:
0.33333333333333333334
3.1415926535897932385
in 60 digits:
0.33333333333333333333333333333333333333333333333333333333333333333334
3.1415926535897932384626433832795028841971693993751058209749445923078

Note that the last number is not necessarily rounded as you would naively expect it to be rounded in the decimal system. But note also, that in both cases you got a couple of extra digits. This is because numbers are internally stored by CLN as chunks of binary digits in order to match your machine's word size and to not waste precision. Thus, on architectures with different word size, the above output might even differ with regard to actually computed digits.

It should be clear that objects of class numeric should be used for constructing numbers or for doing arithmetic with them. The objects one deals with most of the time are the polymorphic expressions ex.

4.6.1 Tests on numbers

Once you have declared some numbers, assigned them to expressions and done some arithmetic with them it is frequently desired to retrieve some kind of information from them like asking whether that number is integer, rational, real or complex. For those cases GiNaC provides several useful methods. (Internally, they fall back to invocations of certain CLN functions.)

As an example, let's construct some rational number, multiply it with some multiple of its denominator and test what comes out:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

// some very important constants:
const numeric twentyone(21);
const numeric ten(10);
const numeric five(5);

int main()
{
    numeric answer = twentyone;

    answer /= five;
    cout << answer.is_integer() << endl;  // false, it's 21/5
    answer *= ten;
    cout << answer.is_integer() << endl;  // true, it's 42 now!
}

Note that the variable answer is constructed here as an integer by numeric's copy constructor but in an intermediate step it holds a rational number represented as integer numerator and integer denominator. When multiplied by 10, the denominator becomes unity and the result is automatically converted to a pure integer again. Internally, the underlying CLN is responsible for this behavior and we refer the reader to CLN's documentation. Suffice to say that the same behavior applies to complex numbers as well as return values of certain functions. Complex numbers are automatically converted to real numbers if the imaginary part becomes zero. The full set of tests that can be applied is listed in the following table.

4.6.2 Numeric functions

The following functions can be applied to numeric objects and will be evaluated immediately:

Most of these functions are also available as symbolic functions that can be used in expressions (see section 4.10 Mathematical functions) or, like gcd(), as polynomial algorithms.

4.6.3 Converting numbers

Sometimes it is desirable to convert a numeric object back to a built-in arithmetic type (int, double, etc.). The numeric class provides a couple of methods for this purpose:

int numeric::to_int() const;
long numeric::to_long() const;
double numeric::to_double() const;
cln::cl_N numeric::to_cl_N() const;

to_int() and to_long() only work when the number they are applied on is an exact integer. Otherwise the program will halt with a message like `Not a 32-bit integer'. to_double() applied on a rational number will return a floating-point approximation. Both to_int()/to_long() and to_double() discard the imaginary part of complex numbers.

4.7 Constants

Constants behave pretty much like symbols except that they return some specific number when the method .evalf() is called.

The predefined known constants are:

4.8 Sums, products and powers

Simple rational expressions are written down in GiNaC pretty much like in other CAS or like expressions involving numerical variables in C. The necessary operators +, -, * and / have been overloaded to achieve this goal. When you run the following code snippet, the constructor for an object of type mul is automatically called to hold the product of a and b and then the constructor for an object of type add is called to hold the sum of that mul object and the number one:

    ...
    symbol a("a"), b("b");
    ex MyTerm = 1+a*b;
    ...

For exponentiation, you have already seen the somewhat clumsy (though C-ish) statement pow(x,2); to represent x squared. This direct construction is necessary since we cannot safely overload the constructor ^ in C++ to construct a power object. If we did, it would have several counterintuitive and undesired effects:

• Due to C's operator precedence, 2*x^2 would be parsed as (2*x)^2.
• Due to the binding of the operator ^, x^a^b would result in (x^a)^b. This would be confusing since most (though not all) other CAS interpret this as x^(a^b).
• Also, expressions involving integer exponents are very frequently used, which makes it even more dangerous to overload ^ since it is then hard to distinguish between the semantics as exponentiation and the one for exclusive or. (It would be embarrassing to return 1 where one has requested 2^3.)

All effects are contrary to mathematical notation and differ from the way most other CAS handle exponentiation, therefore overloading ^ is ruled out for GiNaC's C++ part. The situation is different in ginsh, there the exponentiation-^ exists. (Also note that the other frequently used exponentiation operator ** does not exist at all in C++).

To be somewhat more precise, objects of the three classes described here, are all containers for other expressions. An object of class power is best viewed as a container with two slots, one for the basis, one for the exponent. All valid GiNaC expressions can be inserted. However, basic transformations like simplifying pow(pow(x,2),3) to x^6 automatically are only performed when this is mathematically possible. If we replace the outer exponent three in the example by some symbols a, the simplification is not safe and will not be performed, since a might be 1/2 and x negative.

Objects of type add and mul are containers with an arbitrary number of slots for expressions to be inserted. Again, simple and safe simplifications are carried out like transforming 3*x+4-x to 2*x+4.

4.9 Lists of expressions

The GiNaC class lst serves for holding a list of arbitrary expressions. They are not as ubiquitous as in many other computer algebra packages, but are sometimes used to supply a variable number of arguments of the same type to GiNaC methods such as subs() and some matrix constructors, so you should have a basic understanding of them.

Lists can be constructed by assigning a comma-separated sequence of expressions:

{
    symbol x("x"), y("y");
    lst l;
    l = x, 2, y, x+y;
    // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y',
    // in that order
    ...

There are also constructors that allow direct creation of lists of up to 16 expressions, which is often more convenient but slightly less efficient:

    ...
    // This produces the same list 'l' as above:
    // lst l(x, 2, y, x+y);
    // lst l = lst(x, 2, y, x+y);
    ...

Use the nops() method to determine the size (number of expressions) of a list and the op() method or the [] operator to access individual elements:

    ...
    cout << l.nops() << endl;                // prints '4'
    cout << l.op(2) << " " << l[0] << endl;  // prints 'y x'
    ...

As with the standard list<T> container, accessing random elements of a lst is generally an operation of order . Faster read-only sequential access to the elements of a list is possible with the iterator types provided by the lst class:

typedef ... lst::const_iterator;
typedef ... lst::const_reverse_iterator;
lst::const_iterator lst::begin() const;
lst::const_iterator lst::end() const;
lst::const_reverse_iterator lst::rbegin() const;
lst::const_reverse_iterator lst::rend() const;

For example, to print the elements of a list individually you can use:

    ...
    // O(N)
    for (lst::const_iterator i = l.begin(); i != l.end(); ++i)
        cout << *i << endl;
    ...

which is one order faster than

    ...
    // O(N^2)
    for (size_t i = 0; i < l.nops(); ++i)
        cout << l.op(i) << endl;
    ...

These iterators also allow you to use some of the algorithms provided by the C++ standard library:

    ...
    // print the elements of the list (requires #include <iterator>)
    std::copy(l.begin(), l.end(), ostream_iterator<ex>(cout, "\n"));

    // sum up the elements of the list (requires #include <numeric>)
    ex sum = std::accumulate(l.begin(), l.end(), ex(0));
    cout << sum << endl;  // prints '2+2*x+2*y'
    ...

lst is one of the few GiNaC classes that allow in-place modifications (the only other one is matrix). You can modify single elements:

    ...
    l[1] = 42;       // l is now {x, 42, y, x+y}
    l.let_op(1) = 7; // l is now {x, 7, y, x+y}
    ...

You can append or prepend an expression to a list with the append() and prepend() methods:

    ...
    l.append(4*x);   // l is now {x, 7, y, x+y, 4*x}
    l.prepend(0);    // l is now {0, x, 7, y, x+y, 4*x}
    ...

You can remove the first or last element of a list with remove_first() and remove_last():

    ...
    l.remove_first();   // l is now {x, 7, y, x+y, 4*x}
    l.remove_last();    // l is now {x, 7, y, x+y}
    ...

You can remove all the elements of a list with remove_all():

    ...
    l.remove_all();     // l is now empty
    ...

You can bring the elements of a list into a canonical order with sort():

    ...
    lst l1, l2;
    l1 = x, 2, y, x+y;
    l2 = 2, x+y, x, y;
    l1.sort();
    l2.sort();
    // l1 and l2 are now equal
    ...

Finally, you can remove all but the first element of consecutive groups of elements with unique():

    ...
    lst l3;
    l3 = x, 2, 2, 2, y, x+y, y+x;
    l3.unique();        // l3 is now {x, 2, y, x+y}
}

4.10 Mathematical functions

There are quite a number of useful functions hard-wired into GiNaC. For instance, all trigonometric and hyperbolic functions are implemented (See section 5.12 Predefined mathematical functions, for a complete list).

These functions (better called pseudofunctions) are all objects of class function. They accept one or more expressions as arguments and return one expression. If the arguments are not numerical, the evaluation of the function may be halted, as it does in the next example, showing how a function returns itself twice and finally an expression that may be really useful:

    ...
    symbol x("x"), y("y");    
    ex foo = x+y/2;
    cout << tgamma(foo) << endl;
     // -> tgamma(x+(1/2)*y)
    ex bar = foo.subs(y==1);
    cout << tgamma(bar) << endl;
     // -> tgamma(x+1/2)
    ex foobar = bar.subs(x==7);
    cout << tgamma(foobar) << endl;
     // -> (135135/128)*Pi^(1/2)
    ...

Besides evaluation most of these functions allow differentiation, series expansion and so on. Read the next chapter in order to learn more about this.

It must be noted that these pseudofunctions are created by inline functions, where the argument list is templated. This means that whenever you call GiNaC::sin(1) it is equivalent to sin(ex(1)) and will therefore not result in a floating point number. Unless of course the function prototype is explicitly overridden -- which is the case for arguments of type numeric (not wrapped inside an ex). Hence, in order to obtain a floating point number of class numeric you should call sin(numeric(1)). This is almost the same as calling sin(1).evalf() except that the latter will return a numeric wrapped inside an ex.

4.11 Relations

Sometimes, a relation holding between two expressions must be stored somehow. The class relational is a convenient container for such purposes. A relation is by definition a container for two ex and a relation between them that signals equality, inequality and so on. They are created by simply using the C++ operators ==, !=, <, <=, > and >= between two expressions.

See section 4.10 Mathematical functions, for examples where various applications of the .subs() method show how objects of class relational are used as arguments. There they provide an intuitive syntax for substitutions. They are also used as arguments to the ex::series method, where the left hand side of the relation specifies the variable to expand in and the right hand side the expansion point. They can also be used for creating systems of equations that are to be solved for unknown variables. But the most common usage of objects of this class is rather inconspicuous in statements of the form if (expand(pow(a+b,2))==a*a+2*a*b+b*b) {...}. Here, an implicit conversion from relational to bool takes place. Note, however, that == here does not perform any simplifications, hence expand() must be called explicitly.

4.12 Integrals

An object of class integral can be used to hold a symbolic integral. If you want to symbolically represent the integral of x*x from 0 to 1, you would write this as

integral(x, 0, 1, x*x)

.eval_integ()

.evalf()

ex integral::relative_integration_error

int integral::max_integration_level

ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f,
const ex & error)

If you know that an expression holds an integral, you can get the integration variable, the left boundary, right boundary and integrant by respectively calling .op(0), .op(1), .op(2), and .op(3). Differentiating integrals with respect to variables works as expected. Note that it makes no sense to differentiate an integral with respect to the integration variable.

4.13 Matrices

A matrix is a two-dimensional array of expressions. The elements of a matrix with rows and columns are accessed with two unsigned indices, the first one in the range 0..., the second one in the range 0....

There are a couple of ways to construct matrices, with or without preset elements. The constructor

matrix::matrix(unsigned r, unsigned c);

creates a matrix with `r' rows and `c' columns with all elements set to zero.

The fastest way to create a matrix with preinitialized elements is to assign a list of comma-separated expressions to an empty matrix (see below for an example). But you can also specify the elements as a (flat) list with

matrix::matrix(unsigned r, unsigned c, const lst & l);

The function

ex lst_to_matrix(const lst & l);

constructs a matrix from a list of lists, each list representing a matrix row.

There is also a set of functions for creating some special types of matrices:

ex diag_matrix(const lst & l);
ex unit_matrix(unsigned x);
ex unit_matrix(unsigned r, unsigned c);
ex symbolic_matrix(unsigned r, unsigned c, const string & base_name);
ex symbolic_matrix(unsigned r, unsigned c, const string & base_name,
                   const string & tex_base_name);

diag_matrix() constructs a diagonal matrix given the list of diagonal elements. unit_matrix() creates an `x' by `x' (or `r' by `c') unit matrix. And finally, symbolic_matrix constructs a matrix filled with newly generated symbols made of the specified base name and the position of each element in the matrix.

Matrix elements can be accessed and set using the parenthesis (function call) operator:

const ex & matrix::operator()(unsigned r, unsigned c) const;
ex & matrix::operator()(unsigned r, unsigned c);

It is also possible to access the matrix elements in a linear fashion with the op() method. But C++-style subscripting with square brackets `[]' is not available.

Here are a couple of examples for constructing matrices:

{
    symbol a("a"), b("b");

    matrix M(2, 2);
    M = a, 0,
        0, b;
    cout << M << endl;
     // -> [[a,0],[0,b]]

    matrix M2(2, 2);
    M2(0, 0) = a;
    M2(1, 1) = b;
    cout << M2 << endl;
     // -> [[a,0],[0,b]]

    cout << matrix(2, 2, lst(a, 0, 0, b)) << endl;
     // -> [[a,0],[0,b]]

    cout << lst_to_matrix(lst(lst(a, 0), lst(0, b))) << endl;
     // -> [[a,0],[0,b]]

    cout << diag_matrix(lst(a, b)) << endl;
     // -> [[a,0],[0,b]]

    cout << unit_matrix(3) << endl;
     // -> [[1,0,0],[0,1,0],[0,0,1]]

    cout << symbolic_matrix(2, 3, "x") << endl;
     // -> [[x00,x01,x02],[x10,x11,x12]]
}

There are three ways to do arithmetic with matrices. The first (and most direct one) is to use the methods provided by the matrix class:

matrix matrix::add(const matrix & other) const;
matrix matrix::sub(const matrix & other) const;
matrix matrix::mul(const matrix & other) const;
matrix matrix::mul_scalar(const ex & other) const;
matrix matrix::pow(const ex & expn) const;
matrix matrix::transpose() const;

All of these methods return the result as a new matrix object. Here is an example that calculates for three matrices , and :

{
    matrix A(2, 2), B(2, 2), C(2, 2);
    A =  1, 2,
         3, 4;
    B = -1, 0,
         2, 1;
    C =  8, 4,
         2, 1;

    matrix result = A.mul(B).sub(C.mul_scalar(2));
    cout << result << endl;
     // -> [[-13,-6],[1,2]]
    ...
}

The second (and probably the most natural) way is to construct an expression containing matrices with the usual arithmetic operators and pow(). For efficiency reasons, expressions with sums, products and powers of matrices are not automatically evaluated in GiNaC. You have to call the method

ex ex::evalm() const;

to obtain the result:

{
    ...
    ex e = A*B - 2*C;
    cout << e << endl;
     // -> [[1,2],[3,4]]*[[-1,0],[2,1]]-2*[[8,4],[2,1]]
    cout << e.evalm() << endl;
     // -> [[-13,-6],[1,2]]
    ...
}

The non-commutativity of the product A*B in this example is automatically recognized by GiNaC. There is no need to use a special operator here. See section 4.15 Non-commutative objects, for more information about dealing with non-commutative expressions.

Finally, you can work with indexed matrices and call simplify_indexed() to perform the arithmetic:

{
    ...
    idx i(symbol("i"), 2), j(symbol("j"), 2), k(symbol("k"), 2);
    e = indexed(A, i, k) * indexed(B, k, j) - 2 * indexed(C, i, j);
    cout << e << endl;
     // -> -2*[[8,4],[2,1]].i.j+[[-1,0],[2,1]].k.j*[[1,2],[3,4]].i.k
    cout << e.simplify_indexed() << endl;
     // -> [[-13,-6],[1,2]].i.j
}

Using indices is most useful when working with rectangular matrices and one-dimensional vectors because you don't have to worry about having to transpose matrices before multiplying them. See section 4.14 Indexed objects, for more information about using matrices with indices, and about indices in general.

The matrix class provides a couple of additional methods for computing determinants, traces, characteristic polynomials and ranks:

ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
ex matrix::trace() const;
ex matrix::charpoly(const ex & lambda) const;
unsigned matrix::rank() const;

The `algo' argument of determinant() allows to select between different algorithms for calculating the determinant. The asymptotic speed (as parametrized by the matrix size) can greatly differ between those algorithms, depending on the nature of the matrix' entries. The possible values are defined in the `flags.h' header file. By default, GiNaC uses a heuristic to automatically select an algorithm that is likely (but not guaranteed) to give the result most quickly.

Matrices may also be inverted using the ex matrix::inverse() method and linear systems may be solved with:

matrix matrix::solve(const matrix & vars, const matrix & rhs,
                     unsigned algo=solve_algo::automatic) const;

Assuming the matrix object this method is applied on is an m times n matrix, then vars must be a n times p matrix of symbolic indeterminates and rhs a m times p matrix. The returned matrix then has dimension n times p and in the case of an underdetermined system will still contain some of the indeterminates from vars. If the system is overdetermined, an exception is thrown.

4.14 Indexed objects

GiNaC allows you to handle expressions containing general indexed objects in arbitrary spaces. It is also able to canonicalize and simplify such expressions and perform symbolic dummy index summations. There are a number of predefined indexed objects provided, like delta and metric tensors.

There are few restrictions placed on indexed objects and their indices and it is easy to construct nonsense expressions, but our intention is to provide a general framework that allows you to implement algorithms with indexed quantities, getting in the way as little as possible.

4.14.1 Indexed quantities and their indices

Indexed expressions in GiNaC are constructed of two special types of objects, index objects and indexed objects.

• Index objects are of class idx or a subclass. Every index has a value and a dimension (which is the dimension of the space the index lives in) which can both be arbitrary expressions but are usually a number or a simple symbol. In addition, indices of class varidx have a variance (they can be co- or contravariant), and indices of class spinidx have a variance and can be dotted or undotted.

• Indexed objects are of class indexed or a subclass. They contain a base expression (which is the expression being indexed), and one or more indices.

Please notice: when printing expressions, covariant indices and indices without variance are denoted `.i' while contravariant indices are denoted `~i'. Dotted indices have a `*' in front of the index value. In the following, we are going to use that notation in the text so instead of we will write `A~i.j.k'. Index dimensions are not visible in the output.

A simple example shall illustrate the concepts:

#include <iostream>
#include <ginac/ginac.h>
using namespace std;
using namespace GiNaC;

int main()
{
    symbol i_sym("i"), j_sym("j");
    idx i(i_sym, 3), j(j_sym, 3);

    symbol A("A");
    cout << indexed(A, i, j) << endl;
     // -> A.i.j
    cout << index_dimensions << indexed(A, i, j) << endl;
     // -> A.i[3].j[3]
    cout << dflt; // reset cout to default output format (dimensions hidden)
    ...

The idx constructor takes two arguments, the index value and the index dimension. First we define two index objects, i and j, both with the numeric dimension 3. The value of the index i is the symbol i_sym (which prints as `i') and the value of the index j is the symbol j_sym (which prints as `j'). Next we construct an expression containing one indexed object, `A.i.j'. It has the symbol A as its base expression and the two indices i and j.

The dimensions of indices are normally not visible in the output, but one can request them to be printed with the index_dimensions manipulator, as shown above.

Note the difference between the indices i and j which are of class idx, and the index values which are the symbols i_sym and j_sym. The indices of indexed objects cannot directly be symbols or numbers but must be index objects. For example, the following is not correct and will raise an exception:

symbol i("i"), j("j");
e = indexed(A, i, j); // ERROR: indices must be of type idx

You can have multiple indexed objects in an expression, index values can be numeric, and index dimensions symbolic:

    ...
    symbol B("B"), dim("dim");
    cout << 4 * indexed(A, i)
          + indexed(B, idx(j_sym, 4), idx(2, 3), idx(i_sym, dim)) << endl;
     // -> B.j.2.i+4*A.i
    ...

B has a 4-dimensional symbolic index `k', a 3-dimensional numeric index of value 2, and a symbolic index `i' with the symbolic dimension `dim'. Note that GiNaC doesn't automatically notify you that the free indices of `A' and `B' in the sum don't match (you have to call simplify_indexed() for that, see below).

In fact, base expressions, index values and index dimensions can be arbitrary expressions:

    ...
    cout << indexed(A+B, idx(2*i_sym+1, dim/2)) << endl;
     // -> (B+A).(1+2*i)
    ...

It's also possible to construct nonsense like `Pi.sin(x)'. You will not get an error message from this but you will probably not be able to do anything useful with it.

The methods

ex idx::get_value();
ex idx::get_dimension();

return the value and dimension of an idx object. If you have an index in an expression, such as returned by calling .op() on an indexed object, you can get a reference to the idx object with the function ex_to<idx>() on the expression.

There are also the methods

bool idx::is_numeric();
bool idx::is_symbolic();
bool idx::is_dim_numeric();
bool idx::is_dim_symbolic();

for checking whether the value and dimension are numeric or symbolic (non-numeric). Using the info() method of an index (see 5.1 Getting information about expressions) returns information about the index value.

If you need co- and contravariant indices, use the varidx class:

    ...
    symbol mu_sym("mu"), nu_sym("nu");
    varidx mu(mu_sym, 4), nu(nu_sym, 4); // default is contravariant ~mu, ~nu
    varidx mu_co(mu_sym, 4, true);       // covariant index .mu

    cout << indexed(A, mu, nu) << endl;
     // -> A~mu~nu
    cout << indexed(A, mu_co, nu) << endl;
     // -> A.mu~nu
    cout << indexed(A, mu.toggle_variance(), nu) << endl;
     // -> A.mu~nu
    ...

A varidx is an idx with an additional flag that marks it as co- or contravariant. The default is a contravariant (upper) index, but this can be overridden by supplying a third argument to the varidx constructor. The two methods

bool varidx::is_covariant();
bool varidx::is_contravariant();

allow you to check the variance of a varidx object (use ex_to<varidx>() to get the object reference from an expression). There's also the very useful method

ex varidx::toggle_variance();

which makes a new index with the same value and dimension but the opposite variance. By using it you only have to define the index once.

The spinidx class provides dotted and undotted variant indices, as used in the Weyl-van-der-Waerden spinor formalism:

    ...
    symbol K("K"), C_sym("C"), D_sym("D");
    spinidx C(C_sym, 2), D(D_sym);          // default is 2-dimensional,
                                            // contravariant, undotted
    spinidx C_co(C_sym, 2, true);           // covariant index
    spinidx D_dot(D_sym, 2, false, true);   // contravariant, dotted
    spinidx D_co_dot(D_sym, 2, true, true); // covariant, dotted

    cout << indexed(K, C, D) << endl;
     // -> K~C~D
    cout << indexed(K, C_co, D_dot) << endl;
     // -> K.C~*D
    cout << indexed(K, D_co_dot, D) << endl;
     // -> K.*D~D
    ...

A spinidx is a varidx with an additional flag that marks it as dotted or undotted. The default is undotted but this can be overridden by supplying a fourth argument to the spinidx constructor. The two methods

bool spinidx::is_dotted();
bool spinidx::is_undotted();

allow you to check whether or not a spinidx object is dotted (use ex_to<spinidx>() to get the object reference from an expression). Finally, the two methods

ex spinidx::toggle_dot();
ex spinidx::toggle_variance_dot();

create a new index with the same value and dimension but opposite dottedness and the same or opposite variance.

4.14.2 Substituting indices

Sometimes you will want to substitute one symbolic index with another symbolic or numeric index, for example when calculating one specific element of a tensor expression. This is done with the .subs() method, as it is done for symbols (see 5.3 Substituting expressions).

You have two possibilities here. You can either substitute the whole index by another index or expression:

    ...
    ex e = indexed(A, mu_co);
    cout << e << " becomes " << e.subs(mu_co == nu) << endl;
     // -> A.mu becomes A~nu
    cout << e << " becomes " << e.subs(mu_co == varidx(0, 4)) << endl;
     // -> A.mu becomes A~0
    cout << e << " becomes " << e.subs(mu_co == 0) << endl;
     // -> A.mu becomes A.0
    ...

The third example shows that trying to replace an index with something that is not an index will substitute the index value instead.

Alternatively, you can substitute the symbol of a symbolic index by another expression:

    ...
    ex e = indexed(A, mu_co);
    cout << e << " becomes " << e.subs(mu_sym == nu_sym) << endl;
     // -> A.mu becomes A.nu
    cout << e << " becomes " << e.subs(mu_sym == 0) << endl;
     // -> A.mu becomes A.0
    ...

As you see, with the second method only the value of the index will get substituted. Its other properties, including its dimension, remain unchanged. If you want to change the dimension of an index you have to substitute the whole index by another one with the new dimension.

Finally, substituting the base expression of an indexed object works as expected:

    ...
    ex e = indexed(A, mu_co);
    cout << e << " becomes " << e.subs(A == A+B) << endl;
     // -> A.mu becomes (B+A).mu
    ...

4.14.3 Symmetries

Indexed objects can have certain symmetry properties with respect to their indices. Symmetries are specified as a tree of objects of class symmetry that is constructed with the helper functions

symmetry sy_none(...);
symmetry sy_symm(...);
symmetry sy_anti(...);
symmetry sy_cycl(...);

sy_none() stands for no symmetry, sy_symm() and sy_anti() specify fully symmetric or antisymmetric, respectively, and sy_cycl() represents a cyclic symmetry. Each of these functions accepts up to four arguments which can be either symmetry objects themselves or unsigned integer numbers that represent an index position (counting from 0). A symmetry specification that consists of only a single sy_symm(), sy_anti() or sy_cycl() with no arguments specifies the respective symmetry for all indices.

Here are some examples of symmetry definitions:

    ...
    // No symmetry:
    e = indexed(A, i, j);
    e = indexed(A, sy_none(), i, j);     // equivalent
    e = indexed(A, sy_none(0, 1), i, j); // equivalent

    // Symmetric in all three indices:
    e = indexed(A, sy_symm(), i, j, k);
    e = indexed(A, sy_symm(0, 1, 2), i, j, k); // equivalent
    e = indexed(A, sy_symm(2, 0, 1), i, j, k); // same symmetry, but yields a
                                               // different canonical order

    // Symmetric in the first two indices only:
    e = indexed(A, sy_symm(0, 1), i, j, k);
    e = indexed(A, sy_none(sy_symm(0, 1), 2), i, j, k); // equivalent

    // Antisymmetric in the first and last index only (index ranges need not
    // be contiguous):
    e = indexed(A, sy_anti(0, 2), i, j, k);
    e = indexed(A, sy_none(sy_anti(0, 2), 1), i, j, k); // equivalent

    // An example of a mixed symmetry: antisymmetric in the first two and
    // last two indices, symmetric when swapping the first and last index
    // pairs (like the Riemann curvature tensor):
    e = indexed(A, sy_symm(sy_anti(0, 1), sy_anti(2, 3)), i, j, k, l);

    // Cyclic symmetry in all three indices:
    e = indexed(A, sy_cycl(), i, j, k);
    e = indexed(A, sy_cycl(0, 1, 2), i, j, k); // equivalent

    // The following examples are invalid constructions that will throw
    // an exception at run time.

    // An index may not appear multiple times:
    e = indexed(A, sy_symm(0, 0, 1), i, j, k); // ERROR
    e = indexed(A, sy_none(sy_symm(0, 1), sy_anti(0, 2)), i, j, k); // ERROR

    // Every child of sy_symm(), sy_anti() and sy_cycl() must refer to the
    // same number of indices:
    e = indexed(A, sy_symm(sy_anti(0, 1), 2), i, j, k); // ERROR

    // And of course, you cannot specify indices which are not there:
    e = indexed(A, sy_symm(0, 1, 2, 3), i, j, k); // ERROR
    ...

If you need to specify more than four indices, you have to use the .add() method of the symmetry class. For example, to specify full symmetry in the first six indices you would write sy_symm(0, 1, 2, 3).add(4).add(5).

If an indexed object has a symmetry, GiNaC will automatically bring the indices into a canonical order which allows for some immediate simplifications:

    ...
    cout << indexed(A, sy_symm(), i, j)
          + indexed(A, sy_symm(), j, i) << endl;
     // -> 2*A.j.i
    cout << indexed(B, sy_anti(), i, j)
          + indexed(B, sy_anti(), j, i) << endl;
     // -> 0
    cout << indexed(B, sy_anti(), i, j, k)
          - indexed(B, sy_anti(), j, k, i) << endl;
     // -> 0
    ...

4.14.4 Dummy indices

GiNaC treats certain symbolic index pairs as dummy indices meaning that a summation over the index range is implied. Symbolic indices which are not dummy indices are called free indices. Numeric indices are neither dummy nor free indices.

To be recognized as a dummy index pair, the two indices must be of the same class and their value must be the same single symbol (an index like `2*n+1' is never a dummy index). If the indices are of class varidx they must also be of opposite variance; if they are of class spinidx they must be both dotted or both undotted.

The method .get_free_indices() returns a vector containing the free indices of an expression. It also checks that the free indices of the terms of a sum are consistent:

{
    symbol A("A"), B("B"), C("C");

    symbol i_sym("i"), j_sym("j"), k_sym("k"), l_sym("l");
    idx i(i_sym, 3), j(j_sym, 3), k(k_sym, 3), l(l_sym, 3);

    ex e = indexed(A, i, j) * indexed(B, j, k) + indexed(C, k, l, i, l);
    cout << exprseq(e.get_free_indices()) << endl;
     // -> (.i,.k)
     // 'j' and 'l' are dummy indices

    symbol mu_sym("mu"), nu_sym("nu"), rho_sym("rho"), sigma_sym("sigma");
    varidx mu(mu_sym, 4), nu(nu_sym, 4), rho(rho_sym, 4), sigma(sigma_sym, 4);

    e = indexed(A, mu, nu) * indexed(B, nu.toggle_variance(), rho)
      + indexed(C, mu, sigma, rho, sigma.toggle_variance());
    cout << exprseq(e.get_free_indices()) << endl;
     // -> (~mu,~rho)
     // 'nu' is a dummy index, but 'sigma' is not

    e = indexed(A, mu, mu);
    cout << exprseq(e.get_free_indices()) << endl;
     // -> (~mu)
     // 'mu' is not a dummy index because it appears twice with the same
     // variance

    e = indexed(A, mu, nu) + 42;
    cout << exprseq(e.get_free_indices()) << endl; // ERROR
     // this will throw an exception:
     // "add::get_free_indices: inconsistent indices in sum"
}

4.14.5 Simplifying indexed expressions

In addition to the few automatic simplifications that GiNaC performs on indexed expressions (such as re-ordering the indices of symmetric tensors and calculating traces and convolutions of matrices and predefined tensors) there is the method

ex ex::simplify_indexed();
ex ex::simplify_indexed(const scalar_products & sp);

that performs some more expensive operations:

• it checks the consistency of free indices in sums in the same way get_free_indices() does
• it tries to give dummy indices that appear in different terms of a sum the same name to allow simplifications like
• it (symbolically) calculates all possible dummy index summations/contractions with the predefined tensors (this will be explained in more detail in the next section)
• it detects contractions that vanish for symmetry reasons, for example the contraction of a symmetric and a totally antisymmetric tensor
• as a special case of dummy index summation, it can replace scalar products of two tensors with a user-defined value

The last point is done with the help of the scalar_products class which is used to store scalar products with known values (this is not an arithmetic class, you just pass it to simplify_indexed()):

{
    symbol A("A"), B("B"), C("C"), i_sym("i");
    idx i(i_sym, 3);

    scalar_products sp;
    sp.add(A, B, 0); // A and B are orthogonal
    sp.add(A, C, 0); // A and C are orthogonal
    sp.add(A, A, 4); // A^2 = 4 (A has length 2)

    e = indexed(A + B, i) * indexed(A + C, i);
    cout << e << endl;
     // -> (B+A).i*(A+C).i

    cout << e.expand(expand_options::expand_indexed).simplify_indexed(sp)
         << endl;
     // -> 4+C.i*B.i
}

The scalar_products object sp acts as a storage for the scalar products added to it with the .add() method. This method takes three arguments: the two expressions of which the scalar product is taken, and the expression to replace it with. After sp.add(A, B, 0), simplify_indexed() will replace all scalar products of indexed objects that have the symbols A and B as base expressions with the single value 0. The number, type and dimension of the indices don't matter; `A~mu~nu*B.mu.nu' would also be replaced by 0.

The example above also illustrates a feature of the expand() method: if passed the expand_indexed option it will distribute indices over sums, so `(A+B).i' becomes `A.i+B.i'.

4.14.6 Predefined tensors

Some frequently used special tensors such as the delta, epsilon and metric tensors are predefined in GiNaC. They have special properties when contracted with other tensor expressions and some of them have constant matrix representations (they will evaluate to a number when numeric indices are specified).

4.14.6.1 Delta tensor

The delta tensor takes two indices, is symmetric and has the matrix representation diag(1, 1, 1, ...). It is constructed by the function delta_tensor():

{
    symbol A("A"), B("B");

    idx i(symbol("i"), 3), j(symbol("j"), 3),
        k(symbol("k"), 3), l(symbol("l"), 3);

    ex e = indexed(A, i, j) * indexed(B, k, l)
         * delta_tensor(i, k) * delta_tensor(j, l);
    cout << e.simplify_indexed() << endl;
     // -> B.i.j*A.i.j

    cout << delta_tensor(i, i) << endl;
     // -> 3
}

4.14.6.2 General metric tensor

The function metric_tensor() creates a general symmetric metric tensor with two indices that can be used to raise/lower tensor indices. The metric tensor is denoted as `g' in the output and if its indices are of mixed variance it is automatically replaced by a delta tensor:

{
    symbol A("A");

    varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);

    ex e = metric_tensor(mu, nu) * indexed(A, nu.toggle_variance(), rho);
    cout << e.simplify_indexed() << endl;
     // -> A~mu~rho

    e = delta_tensor(mu, nu.toggle_variance()) * metric_tensor(nu, rho);
    cout << e.simplify_indexed() << endl;
     // -> g~mu~rho

    e = metric_tensor(mu.toggle_variance(), nu.toggle_variance())
      * metric_tensor(nu, rho);
    cout << e.simplify_indexed() << endl;
     // -> delta.mu~rho

    e = metric_tensor(nu.toggle_variance(), rho.toggle_variance())
      * metric_tensor(mu, nu) * (delta_tensor(mu.toggle_variance(), rho)
        + indexed(A, mu.toggle_variance(), rho));
    cout << e.simplify_indexed() << endl;
     // -> 4+A.rho~rho
}

4.14.6.3 Minkowski metric tensor

The Minkowski metric tensor is a special metric tensor with a constant matrix representation which is either diag(1, -1, -1, ...) (negative signature, the default) or diag(-1, 1, 1, ...) (positive signature). It is created with the function lorentz_g() (although it is output as `eta'):

{
    varidx mu(symbol("mu"), 4);

    e = delta_tensor(varidx(0, 4), mu.toggle_variance())
      * lorentz_g(mu, varidx(0, 4));       // negative signature
    cout << e.simplify_indexed() << endl;
     // -> 1

    e = delta_tensor(varidx(0, 4), mu.toggle_variance())
      * lorentz_g(mu, varidx(0, 4), true); // positive signature
    cout << e.simplify_indexed() << endl;
     // -> -1
}

4.14.6.4 Spinor metric tensor

The function spinor_metric() creates an antisymmetric tensor with two indices that is used to raise/lower indices of 2-component spinors. It is output as `eps':

{
    symbol psi("psi");

    spinidx A(symbol("A")), B(symbol("B")), C(symbol("C"));
    ex A_co = A.toggle_variance(), B_co = B.toggle_variance();

    e = spinor_metric(A, B) * indexed(psi, B_co);
    cout << e.simplify_indexed() << endl;
     // -> psi~A

    e = spinor_metric(A, B) * indexed(psi, A_co);
    cout << e.simplify_indexed() << endl;
     // -> -psi~B

    e = spinor_metric(A_co, B_co) * indexed(psi, B);
    cout << e.simplify_indexed() << endl;
     // -> -psi.A

    e = spinor_metric(A_co, B_co) * indexed(psi, A);
    cout << e.simplify_indexed() << endl;
     // -> psi.B

    e = spinor_metric(A_co, B_co) * spinor_metric(A, B);
    cout << e.simplify_indexed() << endl;
     // -> 2

    e = spinor_metric(A_co, B_co) * spinor_metric(B, C);
    cout << e.simplify_indexed() << endl;
     // -> -delta.A~C
}

The matrix representation of the spinor metric is [[0, 1], [-1, 0]].

4.14.6.5 Epsilon tensor

The epsilon tensor is totally antisymmetric, its number of indices is equal to the dimension of the index space (the indices must all be of the same numeric dimension), and `eps.1.2.3...' (resp. `eps~0~1~2...') is defined to be 1. Its behavior with indices that have a variance also depends on the signature of the metric. Epsilon tensors are output as `eps'.

There are three functions defined to create epsilon tensors in 2, 3 and 4 dimensions:

ex epsilon_tensor(const ex & i1, const ex & i2);
ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3);
ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4,
               bool pos_sig = false);

The first two functions create an epsilon tensor in 2 or 3 Euclidean dimensions, the last function creates an epsilon tensor in a 4-dimensional Minkowski space (the last bool argument specifies whether the metric has negative or positive signature, as in the case of the Minkowski metric tensor):

{
    varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4),
           sig(symbol("sig"), 4), lam(symbol("lam"), 4), bet(symbol("bet"), 4);
    e = lorentz_eps(mu, nu, rho, sig) *
        lorentz_eps(mu.toggle_variance(), nu.toggle_variance(), lam, bet);
    cout << simplify_indexed(e) << endl;
     // -> 2*eta~bet~rho*eta~sig~lam-2*eta~sig~bet*eta~rho~lam

    idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
    symbol A("A"), B("B");
    e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(B, k);
    cout << simplify_indexed(e) << endl;
     // -> -B.k*A.j*eps.i.k.j
    e = epsilon_tensor(i, j, k) * indexed(A, j) * indexed(A, k);
    cout << simplify_indexed(e) << endl;
     // -> 0
}

4.14.7 Linear algebra

The matrix class can be used with indices to do some simple linear algebra (linear combinations and products of vectors and matrices, traces and scalar products):

{
    idx i(symbol("i"), 2), j(symbol("j"), 2);
    symbol x("x"), y("y");

    // A is a 2x2 matrix, X is a 2x1 vector
    matrix A(2, 2), X(2, 1);
    A = 1, 2,
        3, 4;
    X = x, y;

    cout << indexed(A, i, i) << endl;
     // -> 5

    ex e = indexed(A, i, j) * indexed(X, j);
    cout << e.simplify_indexed() << endl;
     // -> [[2*y+x],[4*y+3*x]].i

    e = indexed(A, i, j) * indexed(X, i) + indexed(X, j) * 2;
    cout << e.simplify_indexed() << endl;
     // -> [[3*y+3*x,6*y+2*x]].j
}

You can of course obtain the same results with the matrix::add(), matrix::mul() and matrix::trace() methods (see section 4.13 Matrices) but with indices you don't have to worry about transposing matrices.

Matrix indices always start at 0 and their dimension must match the number of rows/columns of the matrix. Matrices with one row or one column are vectors and can have one or two indices (it doesn't matter whether it's a row or a column vector). Other matrices must have two indices.

You should be careful when using indices with variance on matrices. GiNaC doesn't look at the variance and doesn't know that `F~mu~nu' and `F.mu.nu' are different matrices. In this case you should use only one form for `F' and explicitly multiply it with a matrix representation of the metric tensor.

4.15 Non-commutative objects

GiNaC is equipped to handle certain non-commutative algebras. Three classes of non-commutative objects are built-in which are mostly of use in high energy physics:

• Clifford (Dirac) algebra (class clifford)
• su(3) Lie algebra (class color)
• Matrices (unindexed) (class matrix)

The clifford and color classes are subclasses of indexed because the elements of these algebras usually carry indices. The matrix class is described in more detail in 4.13 Matrices.

Unlike most computer algebra systems, GiNaC does not primarily provide an operator (often denoted `&*') for representing inert products of arbitrary objects. Rather, non-commutativity in GiNaC is a property of the classes of objects involved, and non-commutative products are formed with the usual `*' operator, as are ordinary products. GiNaC is capable of figuring out by itself which objects commutate and will group the factors by their class. Consider this example:

    ...
    varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
    idx a(symbol("a"), 8), b(symbol("b"), 8);
    ex e = -dirac_gamma(mu) * (2*color_T(a)) * 8 * color_T(b) * dirac_gamma(nu);
    cout << e << endl;
     // -> -16*(gamma~mu*gamma~nu)*(T.a*T.b)
    ...

As can be seen, GiNaC pulls out the overall commutative factor `-16' and groups the non-commutative factors (the gammas and the su(3) generators) together while preserving the order of factors within each class (because Clifford objects commutate with color objects). The resulting expression is a commutative product with two factors that are themselves non-commutative products (`gamma~mu*gamma~nu' and `T.a*T.b'). For clarification, parentheses are placed around the non-commutative products in the output.

Non-commutative products are internally represented by objects of the class ncmul, as opposed to commutative products which are handled by the mul class. You will normally not have to worry about this distinction, though.

The advantage of this approach is that you never have to worry about using (or forgetting to use) a special operator when constructing non-commutative expressions. Also, non-commutative products in GiNaC are more intelligent than in other computer algebra systems; they can, for example, automatically canonicalize themselves according to rules specified in the implementation of the non-commutative classes. The drawback is that to work with other than the built-in algebras you have to implement new classes yourself. Symbols always commutate and it's not possible to construct non-commutative products using symbols to represent the algebra elements or generators. User-defined functions can, however, be specified as being non-commutative.

Information about the commutativity of an object or expression can be obtained with the two member functions

unsigned ex::return_type() const;
unsigned ex::return_type_tinfo() const;

The return_type() function returns one of three values (defined in the header file `flags.h'), corresponding to three categories of expressions in GiNaC:

• return_types::commutative: Commutates with everything. Most GiNaC classes are of this kind.
• return_types::noncommutative: Non-commutative, belonging to a certain class of non-commutative objects which can be determined with the return_type_tinfo() method. Expressions of this category commutate with everything except noncommutative expressions of the same class.
• return_types::noncommutative_composite: Non-commutative, composed of non-commutative objects of different classes. Expressions of this category don't commutate with any other noncommutative or noncommutative_composite expressions.

The value returned by the return_type_tinfo() method is valid only when the return type of the expression is noncommutative. It is a value that is unique to the class of the object and usually one of the constants in `tinfos.h', or derived therefrom.

Here are a couple of examples:

Note: the return_type_tinfo() of Clifford objects is only equal to TINFO_clifford for objects with a representation label of zero. Other representation labels yield a different return_type_tinfo(), but it's the same for any two objects with the same label. This is also true for color objects.

A last note: With the exception of matrices, positive integer powers of non-commutative objects are automatically expanded in GiNaC. For example, pow(a*b, 2) becomes `a*b*a*b' if `a' and `b' are non-commutative expressions).

4.15.1 Clifford algebra

Clifford algebras are supported in two flavours: Dirac gamma matrices (more physical) and generic Clifford algebras (more mathematical).

4.15.1.1 Dirac gamma matrices

ex dirac_gamma(const ex & mu, unsigned char rl = 0);

which takes two arguments: the index and a representation label in the range 0 to 255 which is used to distinguish elements of different Clifford algebras (this is also called a spin line index). Gammas with different labels commutate with each other. The dimension of the index can be 4 or (in the framework of dimensional regularization) any symbolic value. Spinor indices on Dirac gammas are not supported in GiNaC.

The unity element of a Clifford algebra is constructed by

ex dirac_ONE(unsigned char rl = 0);

Please notice: You must always use dirac_ONE() when referring to multiples of the unity element, even though it's customary to omit it. E.g. instead of dirac_gamma(mu)*(dirac_slash(q,4)+m) you have to write dirac_gamma(mu)*(dirac_slash(q,4)+m*dirac_ONE()). Otherwise, GiNaC will complain and/or produce incorrect results.

There is a special element `gamma5' that commutates with all other gammas, has a unit square, and in 4 dimensions equals `gamma~0 gamma~1 gamma~2 gamma~3', provided by

ex dirac_gamma5(unsigned char rl = 0);

The chiral projectors `(1+/-gamma5)/2' are also available as proper objects, constructed by

ex dirac_gammaL(unsigned char rl = 0);
ex dirac_gammaR(unsigned char rl = 0);

They observe the relations `gammaL^2 = gammaL', `gammaR^2 = gammaR', and `gammaL gammaR = gammaR gammaL = 0'.

Finally, the function

ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);

creates a term that represents a contraction of `e' with the Dirac Lorentz vector (it behaves like a term of the form `e.mu gamma~mu' with a unique index whose dimension is given by the dim argument). Such slashed expressions are printed with a trailing backslash, e.g. `e\'.

In products of dirac gammas, superfluous unity elements are automatically removed, squares are replaced by their values, and `gamma5', `gammaL' and `gammaR' are moved to the front.

The simplify_indexed() function performs contractions in gamma strings, for example

{
    ...
    symbol a("a"), b("b"), D("D");
    varidx mu(symbol("mu"), D);
    ex e = dirac_gamma(mu) * dirac_slash(a, D)
         * dirac_gamma(mu.toggle_variance());
    cout << e << endl;
     // -> gamma~mu*a\*gamma.mu
    e = e.simplify_indexed();
    cout << e << endl;
     // -> -D*a\+2*a\
    cout << e.subs(D == 4) << endl;
     // -> -2*a\
    ...
}

To calculate the trace of an expression containing strings of Dirac gammas you use one of the functions

ex dirac_trace(const ex & e, const std::set<unsigned char> & rls,
               const ex & trONE = 4);
ex dirac_trace(const ex & e, const lst & rll, const ex & trONE = 4);
ex dirac_trace(const ex & e, unsigned char rl = 0, const ex & trONE = 4);

These functions take the trace over all gammas in the specified set rls or list rll of representation labels, or the single label rl; gammas with other labels are left standing. The last argument to dirac_trace() is the value to be returned for the trace of the unity element, which defaults to 4.

The dirac_trace() function is a linear functional that is equal to the ordinary matrix trace only in dimensions. In particular, the functional is not cyclic in dimensions when acting on expressions containing `gamma5', so it's not a proper trace. This `gamma5' scheme is described in greater detail in The Role of gamma5 in Dimensional Regularization.

The value of the trace itself is also usually different in 4 and in dimensions:

{
    // 4 dimensions
    varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
    ex e = dirac_gamma(mu) * dirac_gamma(nu) *
           dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
    cout << dirac_trace(e).simplify_indexed() << endl;
     // -> -8*eta~rho~nu
}
...
{
    // D dimensions
    symbol D("D");
    varidx mu(symbol("mu"), D), nu(symbol("nu"), D), rho(symbol("rho"), D);
    ex e = dirac_gamma(mu) * dirac_gamma(nu) *
           dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
    cout << dirac_trace(e).simplify_indexed() << endl;
     // -> 8*eta~rho~nu-4*eta~rho~nu*D
}

Here is an example for using dirac_trace() to compute a value that appears in the calculation of the one-loop vacuum polarization amplitude in QED:

{
    symbol q("q"), l("l"), m("m"), ldotq("ldotq"), D("D");
    varidx mu(symbol("mu"), D), nu(symbol("nu"), D);

    scalar_products sp;
    sp.add(l, l, pow(l, 2));
    sp.add(l, q, ldotq);

    ex e = dirac_gamma(mu) *
           (dirac_slash(l, D) + dirac_slash(q, D) + m * dirac_ONE()) *    
           dirac_gamma(mu.toggle_variance()) *
           (dirac_slash(l, D) + m * dirac_ONE());   
    e = dirac_trace(e).simplify_indexed(sp);
    e = e.collect(lst(l, ldotq, m));
    cout << e << endl;
     // -> (8-4*D)*l^2+(8-4*D)*ldotq+4*D*m^2
}

The canonicalize_clifford() function reorders all gamma products that appear in an expression to a canonical (but not necessarily simple) form. You can use this to compare two expressions or for further simplifications:

{
    varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
    ex e = dirac_gamma(mu) * dirac_gamma(nu) + dirac_gamma(nu) * dirac_gamma(mu);
    cout << e << endl;
     // -> gamma~mu*gamma~nu+gamma~nu*gamma~mu

    e = canonicalize_clifford(e);
    cout << e << endl;
     // -> 2*ONE*eta~mu~nu
}

4.15.1.2 A generic Clifford algebra

A generic Clifford algebra, i.e. a dimensional algebra with generators satisfying the identities e~i e~j + e~j e~i = M(i, j) for some matrix (metric) , which may be non-symmetric and containing symbolic entries. Such generators are created by the function

    ex clifford_unit(const ex & mu, const ex & metr, unsigned char rl = 0);

where mu should be a varidx class object indexing the generators, metr defines the metric and can be represented by a square matrix, tensormetric or indexed class object, optional parameter rl allows to distinguish different Clifford algebras (which will commute with each other). Note that the call clifford_unit(mu, minkmetric()) creates something very close to dirac_gamma(mu). The method clifford::get_metric() returns a metric defining this Clifford number.

If the matrix is in fact symmetric you may prefer to create the Clifford algebra units with a call like that

    ex e = clifford_unit(mu, indexed(M, sy_symm(), i, j));

since this may yield some further automatic simplifications.

Individual generators of a Clifford algebra can be accessed in several ways. For example

{
    ... 
    varidx nu(symbol("nu"), 4);
    realsymbol s("s");
    ex M = diag_matrix(lst(1, -1, 0, s));
    ex e = clifford_unit(nu, M);
    ex e0 = e.subs(nu == 0);
    ex e1 = e.subs(nu == 1);
    ex e2 = e.subs(nu == 2);
    ex e3 = e.subs(nu == 3);
    ...
}

will produce four anti-commuting generators of a Clifford algebra with properties pow(e0, 2) = 1, pow(e1, 2) = -1, pow(e2, 2) = 0 and pow(e3, 2) = s.

A similar effect can be achieved from the function

    ex lst_to_clifford(const ex & v, const ex & mu,  const ex & metr,
                       unsigned char rl = 0);
    ex lst_to_clifford(const ex & v, const ex & e);

which converts a list or vector `v = (v~0, v~1, ..., v~n)' into the Clifford number `v~0 e.0 + v~1 e.1 + ... + v~n e.n' with `e.k' directly supplied in the second form of the procedure. In the first form the Clifford unit `e.k' is generated by the call of clifford_unit(mu, metr, rl). The previous code may be rewritten with the help of lst_to_clifford() as follows

{
    ...
    varidx nu(symbol("nu"), 4);
    realsymbol s("s");
    ex M = diag_matrix(lst(1, -1, 0, s));
    ex e0 = lst_to_clifford(lst(1, 0, 0, 0), nu, M);
    ex e1 = lst_to_clifford(lst(0, 1, 0, 0), nu, M);
    ex e2 = lst_to_clifford(lst(0, 0, 1, 0), nu, M);
    ex e3 = lst_to_clifford(lst(0, 0, 0, 1), nu, M);
  ...
}

There is the inverse function

    lst clifford_to_lst(const ex & e, const ex & c, bool algebraic = true);

which takes an expression e and tries to find a list `v = (v~0, v~1, ..., v~n)' such that `e = v~0 c.0 + v~1 c.1 + ... + v~n c.n' with respect to the given Clifford units c and with none of the `v~k' containing Clifford units c (of course, this may be impossible). This function can use an algebraic method (default) or a symbolic one. With the algebraic method the `v~k' are calculated as `(e c.k + c.k e)/pow(c.k, 2)'. If `pow(c.k, 2)' is zero or is not a numeric for some `k' then the method will be automatically changed to symbolic. The same effect is obtained by the assignment (algebraic = false) in the procedure call.

There are several functions for (anti-)automorphisms of Clifford algebras:

    ex clifford_prime(const ex & e)
    inline ex clifford_star(const ex & e) { return e.conjugate(); }
    inline ex clifford_bar(const ex & e) { return clifford_prime(e.conjugate()); }

The automorphism of a Clifford algebra clifford_prime() simply changes signs of all Clifford units in the expression. The reversion of a Clifford algebra clifford_star() coincides with the conjugate() method and effectively reverses the order of Clifford units in any product. Finally the main anti-automorphism of a Clifford algebra clifford_bar() is the composition of the previous two, i.e. it makes the reversion and changes signs of all Clifford units in a product. These functions correspond to the notations , e* and \bar{e} used in Clifford algebra textbooks.

The function

    ex clifford_norm(const ex & e);

calculates the norm of a Clifford number from the expression ||e||^2 = e \bar{e} The inverse of a Clifford expression is returned by the function

    ex clifford_inverse(const ex & e);

which calculates it as If then an exception is raised.

If a Clifford number happens to be a factor of dirac_ONE() then we can convert it to a "real" (non-Clifford) expression by the function

    ex remove_dirac_ONE(const ex & e);

The function canonicalize_clifford() works for a generic Clifford algebra in a similar way as for Dirac gammas.

The last provided function is

    ex clifford_moebius_map(const ex & a, const ex & b, const ex & c,
                            const ex & d, const ex & v, const ex & G,
                            unsigned char rl = 0);
    ex clifford_moebius_map(const ex & M, const ex & v, const ex & G,
                            unsigned char rl = 0);

It takes a list or vector v and makes the Moebius (conformal or linear-fractional) transformation `v -> (av+b)/(cv+d)' defined by the matrix `M = [[a, b], [c, d]]'. The parameter G defines the metric of the surrounding (pseudo-)Euclidean space. This can be a matrix or a Clifford unit, in the later case the parameter rl is ignored even if supplied. The returned value of this function is a list of components of the resulting vector.

LaTeX output for Clifford units looks like \clifford[1]{e}^{{\nu}}, where 1 is the representation_label and \nu is the index of the corresponding unit. This provides a flexible typesetting with a suitable defintion of the \clifford command. For example, the definition

    \newcommand{\clifford}[1][]{}

    \newcommand{\clifford}[2][]{\ifcase #1 #2\or \tilde{#2} \or \breve{#2} \fi}

4.15.2 Color algebra

For computations in quantum chromodynamics, GiNaC implements the base elements and structure constants of the su(3) Lie algebra (color algebra). The base elements are constructed by the function

ex color_T(const ex & a, unsigned char rl = 0);

which takes two arguments: the index and a representation label in the range 0 to 255 which is used to distinguish elements of different color algebras. Objects with different labels commutate with each other. The dimension of the index must be exactly 8 and it should be of class idx, not varidx.

The unity element of a color algebra is constructed by

ex color_ONE(unsigned char rl = 0);

Please notice: You must always use color_ONE() when referring to multiples of the unity element, even though it's customary to omit it. E.g. instead of color_T(a)*(color_T(b)*indexed(X,b)+1) you have to write color_T(a)*(color_T(b)*indexed(X,b)+color_ONE()). Otherwise, GiNaC may produce incorrect results.

The functions

ex color_d(const ex & a, const ex & b, const ex & c);
ex color_f(const ex & a, const ex & b, const ex & c);

create the symmetric and antisymmetric structure constants and which satisfy and .

These functions evaluate to their numerical values, if you supply numeric indices to them. The index values should be in the range from 1 to 8, not from 0 to 7. This departure from usual conventions goes along better with the notations used in physical literature.

There's an additional function

ex color_h(const ex & a, const ex & b, const ex & c);

which returns the linear combination `color_d(a, b, c)+I*color_f(a, b, c)'.

The function simplify_indexed() performs some simplifications on expressions containing color objects:

{
    ...
    idx a(symbol("a"), 8), b(symbol("b"), 8), c(symbol("c"), 8),
        k(symbol("k"), 8), l(symbol("l"), 8);

    e = color_d(a, b, l) * color_f(a, b, k);
    cout << e.simplify_indexed() << endl;
     // -> 0

    e = color_d(a, b, l) * color_d(a, b, k);
    cout << e.simplify_indexed() << endl;
     // -> 5/3*delta.k.l

    e = color_f(l, a, b) * color_f(a, b, k);
    cout << e.simplify_indexed() << endl;
     // -> 3*delta.k.l

    e = color_h(a, b, c) * color_h(a, b, c);
    cout << e.simplify_indexed() << endl;
     // -> -32/3

    e = color_h(a, b, c) * color_T(b) * color_T(c);
    cout << e.simplify_indexed() << endl;
     // -> -2/3*T.a

    e = color_h(a, b, c) * color_T(a) * color_T(b) * color_T(c);
    cout << e.simplify_indexed() << endl;
     // -> -8/9*ONE

    e = color_T(k) * color_T(a) * color_T(b) * color_T(k);
    cout << e.simplify_indexed() << endl;
     // -> 1/4*delta.b.a*ONE-1/6*T.a*T.b
    ...

To calculate the trace of an expression containing color objects you use one of the functions

ex color_trace(const ex & e, const std::set<unsigned char> & rls);
ex color_trace(const ex & e, const lst & rll);
ex color_trace(const ex & e, unsigned char rl = 0);

These functions take the trace over all color `T' objects in the specified set rls or list rll of representation labels, or the single label rl; `T's with other labels are left standing. For example:

    ...
    e = color_trace(4 * color_T(a) * color_T(b) * color_T(c));
    cout << e << endl;
     // -> -I*f.a.c.b+d.a.c.b
}

4.16 Hash Maps

For your convenience, GiNaC offers the container template exhashmap<T> that can be used as a drop-in replacement for the STL std::map<ex, T, ex_is_less>, using hash tables to provide faster, typically constant-time, element look-up than map<>.

exhashmap<> supports all map<> members and operations, with the following differences:

• no lower_bound() and upper_bound() methods
• no reverse iterators, no rbegin()/rend()
• no operator<(exhashmap, exhashmap)
• the comparison function object key_compare is hardcoded to ex_is_less
• the constructor exhashmap(size_t n) allows specifying the minimum initial hash table size (the actual table size after construction may be larger than the specified value)
• the method size_t bucket_count() returns the current size of the hash table
• insert() and erase() operations invalidate all iterators