+@example
+ ...
+ l.remove_first(); // l is now @{x, 7, y, x+y, 4*x@}
+ l.remove_last(); // l is now @{x, 7, y, x+y@}
+ ...
+@end example
+
+You can remove all the elements of a list with @code{remove_all()}:
+
+@example
+ ...
+ l.remove_all(); // l is now empty
+ ...
+@end example
+
+You can bring the elements of a list into a canonical order with @code{sort()}:
+
+@example
+ ...
+ 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
+ ...
+@end example
+
+Finally, you can remove all but the first element of consecutive groups of
+elements with @code{unique()}:
+
+@example
+ ...
+ lst l3;
+ l3 = x, 2, 2, 2, y, x+y, y+x;
+ l3.unique(); // l3 is now @{x, 2, y, x+y@}
+@}
+@end example
+
+
+@node Mathematical functions, Relations, Lists, Basic concepts
+@c node-name, next, previous, up
+@section Mathematical functions
+@cindex @code{function} (class)
+@cindex trigonometric function
+@cindex hyperbolic function
+
+There are quite a number of useful functions hard-wired into GiNaC. For
+instance, all trigonometric and hyperbolic functions are implemented
+(@xref{Built-in functions}, for a complete list).
+
+These functions (better called @emph{pseudofunctions}) are all objects
+of class @code{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:
+
+@cindex Gamma function
+@cindex @code{subs()}
+@example
+ ...
+ 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)
+ ...
+@end example
+
+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 @code{GiNaC::sin(1)} it is equivalent to
+@code{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 @code{numeric}
+(not wrapped inside an @code{ex}). Hence, in order to obtain a floating
+point number of class @code{numeric} you should call
+@code{sin(numeric(1))}. This is almost the same as calling
+@code{sin(1).evalf()} except that the latter will return a numeric
+wrapped inside an @code{ex}.
+
+
+@node Relations, Integrals, Mathematical functions, Basic concepts
+@c node-name, next, previous, up
+@section Relations
+@cindex @code{relational} (class)
+
+Sometimes, a relation holding between two expressions must be stored
+somehow. The class @code{relational} is a convenient container for such
+purposes. A relation is by definition a container for two @code{ex} and
+a relation between them that signals equality, inequality and so on.
+They are created by simply using the C++ operators @code{==}, @code{!=},
+@code{<}, @code{<=}, @code{>} and @code{>=} between two expressions.
+
+@xref{Mathematical functions}, for examples where various applications
+of the @code{.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 @code{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 @code{if
+(expand(pow(a+b,2))==a*a+2*a*b+b*b) @{...@}}. Here, an implicit
+conversion from @code{relational} to @code{bool} takes place. Note,
+however, that @code{==} here does not perform any simplifications, hence
+@code{expand()} must be called explicitly.
+
+@node Integrals, Matrices, Relations, Basic concepts
+@c node-name, next, previous, up
+@section Integrals
+@cindex @code{integral} (class)
+
+An object of class @dfn{integral} can be used to hold a symbolic integral.
+If you want to symbolically represent the integral of @code{x*x} from 0 to
+1, you would write this as
+@example
+integral(x, 0, 1, x*x)
+@end example
+The first argument is the integration variable. It should be noted that
+GiNaC is not very good (yet?) at symbolically evaluating integrals. In
+fact, it can only integrate polynomials. An expression containing integrals
+can be evaluated symbolically by calling the
+@example
+.eval_integ()
+@end example
+method on it. Numerical evaluation is available by calling the
+@example
+.evalf()
+@end example
+method on an expression containing the integral. This will only evaluate
+integrals into a number if @code{subs}ing the integration variable by a
+number in the fourth argument of an integral and then @code{evalf}ing the
+result always results in a number. Of course, also the boundaries of the
+integration domain must @code{evalf} into numbers. It should be noted that
+trying to @code{evalf} a function with discontinuities in the integration
+domain is not recommended. The accuracy of the numeric evaluation of
+integrals is determined by the static member variable
+@example
+ex integral::relative_integration_error
+@end example
+of the class @code{integral}. The default value of this is 10^-8.
+The integration works by halving the interval of integration, until numeric
+stability of the answer indicates that the requested accuracy has been
+reached. The maximum depth of the halving can be set via the static member
+variable
+@example
+int integral::max_integration_level
+@end example
+The default value is 15. If this depth is exceeded, @code{evalf} will simply
+return the integral unevaluated. The function that performs the numerical
+evaluation, is also available as
+@example
+ex adaptivesimpson(const ex & x, const ex & a, const ex & b, const ex & f,
+ const ex & error)
+@end example
+This function will throw an exception if the maximum depth is exceeded. The
+last parameter of the function is optional and defaults to the
+@code{relative_integration_error}. To make sure that we do not do too
+much work if an expression contains the same integral multiple times,
+a lookup table is used.
+
+If you know that an expression holds an integral, you can get the
+integration variable, the left boundary, right boundary and integrand by
+respectively calling @code{.op(0)}, @code{.op(1)}, @code{.op(2)}, and
+@code{.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.
+
+@node Matrices, Indexed objects, Integrals, Basic concepts
+@c node-name, next, previous, up
+@section Matrices
+@cindex @code{matrix} (class)
+
+A @dfn{matrix} is a two-dimensional array of expressions. The elements of a
+matrix with @math{m} rows and @math{n} columns are accessed with two
+@code{unsigned} indices, the first one in the range 0@dots{}@math{m-1}, the
+second one in the range 0@dots{}@math{n-1}.
+
+There are a couple of ways to construct matrices, with or without preset
+elements. The constructor
+
+@example
+matrix::matrix(unsigned r, unsigned c);
+@end example
+
+creates a matrix with @samp{r} rows and @samp{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
+
+@example
+matrix::matrix(unsigned r, unsigned c, const lst & l);
+@end example
+
+The function
+
+@cindex @code{lst_to_matrix()}
+@example
+ex lst_to_matrix(const lst & l);
+@end example
+
+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:
+
+@cindex @code{diag_matrix()}
+@cindex @code{unit_matrix()}
+@cindex @code{symbolic_matrix()}
+@example
+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);
+@end example
+
+@code{diag_matrix()} constructs a diagonal matrix given the list of diagonal
+elements. @code{unit_matrix()} creates an @samp{x} by @samp{x} (or @samp{r}
+by @samp{c}) unit matrix. And finally, @code{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.
+
+Matrices often arise by omitting elements of another matrix. For
+instance, the submatrix @code{S} of a matrix @code{M} takes a
+rectangular block from @code{M}. The reduced matrix @code{R} is defined
+by removing one row and one column from a matrix @code{M}. (The
+determinant of a reduced matrix is called a @emph{Minor} of @code{M} and
+can be used for computing the inverse using Cramer's rule.)
+
+@cindex @code{sub_matrix()}
+@cindex @code{reduced_matrix()}
+@example
+ex sub_matrix(const matrix&m, unsigned r, unsigned nr, unsigned c, unsigned nc);
+ex reduced_matrix(const matrix& m, unsigned r, unsigned c);
+@end example
+
+The function @code{sub_matrix()} takes a row offset @code{r} and a
+column offset @code{c} and takes a block of @code{nr} rows and @code{nc}
+columns. The function @code{reduced_matrix()} has two integer arguments
+that specify which row and column to remove:
+
+@example
+@{
+ matrix m(3,3);
+ m = 11, 12, 13,
+ 21, 22, 23,
+ 31, 32, 33;
+ cout << reduced_matrix(m, 1, 1) << endl;
+ // -> [[11,13],[31,33]]
+ cout << sub_matrix(m, 1, 2, 1, 2) << endl;
+ // -> [[22,23],[32,33]]
+@}
+@end example
+
+Matrix elements can be accessed and set using the parenthesis (function call)
+operator:
+
+@example
+const ex & matrix::operator()(unsigned r, unsigned c) const;
+ex & matrix::operator()(unsigned r, unsigned c);
+@end example
+
+It is also possible to access the matrix elements in a linear fashion with
+the @code{op()} method. But C++-style subscripting with square brackets
+@samp{[]} is not available.
+
+Here are a couple of examples for constructing matrices:
+
+@example
+@{
+ 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]]
+@}
+@end example
+
+@cindex @code{is_zero_matrix()}
+The method @code{matrix::is_zero_matrix()} returns @code{true} only if
+all entries of the matrix are zeros. There is also method
+@code{ex::is_zero_matrix()} which returns @code{true} only if the
+expression is zero or a zero matrix.
+
+@cindex @code{transpose()}
+There are three ways to do arithmetic with matrices. The first (and most
+direct one) is to use the methods provided by the @code{matrix} class:
+
+@example
+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;
+@end example
+
+All of these methods return the result as a new matrix object. Here is an
+example that calculates @math{A*B-2*C} for three matrices @math{A}, @math{B}
+and @math{C}:
+
+@example
+@{
+ 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]]
+ ...
+@}
+@end example
+
+@cindex @code{evalm()}
+The second (and probably the most natural) way is to construct an expression
+containing matrices with the usual arithmetic operators and @code{pow()}.
+For efficiency reasons, expressions with sums, products and powers of
+matrices are not automatically evaluated in GiNaC. You have to call the
+method
+
+@example
+ex ex::evalm() const;
+@end example
+
+to obtain the result:
+
+@example
+@{
+ ...
+ 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]]
+ ...
+@}
+@end example
+
+The non-commutativity of the product @code{A*B} in this example is
+automatically recognized by GiNaC. There is no need to use a special
+operator here. @xref{Non-commutative objects}, for more information about
+dealing with non-commutative expressions.
+
+Finally, you can work with indexed matrices and call @code{simplify_indexed()}
+to perform the arithmetic:
+
+@example
+@{
+ ...
+ 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
+@}
+@end example
+
+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. @xref{Indexed objects}, for
+more information about using matrices with indices, and about indices in
+general.
+
+The @code{matrix} class provides a couple of additional methods for
+computing determinants, traces, characteristic polynomials and ranks:
+
+@cindex @code{determinant()}
+@cindex @code{trace()}
+@cindex @code{charpoly()}
+@cindex @code{rank()}
+@example
+ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
+ex matrix::trace() const;
+ex matrix::charpoly(const ex & lambda) const;
+unsigned matrix::rank() const;
+@end example
+
+The @samp{algo} argument of @code{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 @file{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.
+
+@cindex @code{inverse()} (matrix)
+@cindex @code{solve()}
+Matrices may also be inverted using the @code{ex matrix::inverse()}
+method and linear systems may be solved with:
+
+@example
+matrix matrix::solve(const matrix & vars, const matrix & rhs,
+ unsigned algo=solve_algo::automatic) const;
+@end example
+
+Assuming the matrix object this method is applied on is an @code{m}
+times @code{n} matrix, then @code{vars} must be a @code{n} times
+@code{p} matrix of symbolic indeterminates and @code{rhs} a @code{m}
+times @code{p} matrix. The returned matrix then has dimension @code{n}
+times @code{p} and in the case of an underdetermined system will still
+contain some of the indeterminates from @code{vars}. If the system is
+overdetermined, an exception is thrown.
+
+
+@node Indexed objects, Non-commutative objects, Matrices, Basic concepts
+@c node-name, next, previous, up
+@section 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.
+
+@cindex @code{idx} (class)
+@cindex @code{indexed} (class)
+@subsection Indexed quantities and their indices
+
+Indexed expressions in GiNaC are constructed of two special types of objects,
+@dfn{index objects} and @dfn{indexed objects}.
+
+@itemize @bullet
+
+@cindex contravariant
+@cindex covariant
+@cindex variance
+@item Index objects are of class @code{idx} or a subclass. Every index has
+a @dfn{value} and a @dfn{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 @code{varidx} have
+a @dfn{variance} (they can be co- or contravariant), and indices of class
+@code{spinidx} have a variance and can be @dfn{dotted} or @dfn{undotted}.
+
+@item Indexed objects are of class @code{indexed} or a subclass. They
+contain a @dfn{base expression} (which is the expression being indexed), and
+one or more indices.
+
+@end itemize
+
+@strong{Please notice:} when printing expressions, covariant indices and indices
+without variance are denoted @samp{.i} while contravariant indices are
+denoted @samp{~i}. Dotted indices have a @samp{*} in front of the index
+value. In the following, we are going to use that notation in the text so
+instead of @math{A^i_jk} we will write @samp{A~i.j.k}. Index dimensions are
+not visible in the output.
+
+A simple example shall illustrate the concepts:
+
+@example
+#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)
+ ...
+@end example
+
+The @code{idx} constructor takes two arguments, the index value and the
+index dimension. First we define two index objects, @code{i} and @code{j},
+both with the numeric dimension 3. The value of the index @code{i} is the
+symbol @code{i_sym} (which prints as @samp{i}) and the value of the index
+@code{j} is the symbol @code{j_sym} (which prints as @samp{j}). Next we
+construct an expression containing one indexed object, @samp{A.i.j}. It has
+the symbol @code{A} as its base expression and the two indices @code{i} and
+@code{j}.
+
+The dimensions of indices are normally not visible in the output, but one
+can request them to be printed with the @code{index_dimensions} manipulator,
+as shown above.
+
+Note the difference between the indices @code{i} and @code{j} which are of
+class @code{idx}, and the index values which are the symbols @code{i_sym}
+and @code{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:
+
+@example
+symbol i("i"), j("j");
+e = indexed(A, i, j); // ERROR: indices must be of type idx
+@end example
+
+You can have multiple indexed objects in an expression, index values can
+be numeric, and index dimensions symbolic:
+
+@example
+ ...
+ 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
+ ...
+@end example
+
+@code{B} has a 4-dimensional symbolic index @samp{k}, a 3-dimensional numeric
+index of value 2, and a symbolic index @samp{i} with the symbolic dimension
+@samp{dim}. Note that GiNaC doesn't automatically notify you that the free
+indices of @samp{A} and @samp{B} in the sum don't match (you have to call
+@code{simplify_indexed()} for that, see below).
+
+In fact, base expressions, index values and index dimensions can be
+arbitrary expressions:
+
+@example
+ ...
+ cout << indexed(A+B, idx(2*i_sym+1, dim/2)) << endl;
+ // -> (B+A).(1+2*i)
+ ...
+@end example
+
+It's also possible to construct nonsense like @samp{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.
+
+@cindex @code{get_value()}
+@cindex @code{get_dimension()}
+The methods
+
+@example
+ex idx::get_value();
+ex idx::get_dimension();
+@end example
+
+return the value and dimension of an @code{idx} object. If you have an index
+in an expression, such as returned by calling @code{.op()} on an indexed
+object, you can get a reference to the @code{idx} object with the function
+@code{ex_to<idx>()} on the expression.
+
+There are also the methods
+
+@example
+bool idx::is_numeric();
+bool idx::is_symbolic();
+bool idx::is_dim_numeric();
+bool idx::is_dim_symbolic();
+@end example
+
+for checking whether the value and dimension are numeric or symbolic
+(non-numeric). Using the @code{info()} method of an index (see @ref{Information
+about expressions}) returns information about the index value.
+
+@cindex @code{varidx} (class)
+If you need co- and contravariant indices, use the @code{varidx} class:
+
+@example
+ ...
+ 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
+ ...
+@end example
+
+A @code{varidx} is an @code{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 @code{varidx}
+constructor. The two methods
+
+@example
+bool varidx::is_covariant();
+bool varidx::is_contravariant();
+@end example
+
+allow you to check the variance of a @code{varidx} object (use @code{ex_to<varidx>()}
+to get the object reference from an expression). There's also the very useful
+method
+
+@example
+ex varidx::toggle_variance();
+@end example
+
+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.
+
+@cindex @code{spinidx} (class)
+The @code{spinidx} class provides dotted and undotted variant indices, as
+used in the Weyl-van-der-Waerden spinor formalism:
+
+@example
+ ...
+ 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
+ ...
+@end example
+
+A @code{spinidx} is a @code{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 @code{spinidx} constructor. The two
+methods
+
+@example
+bool spinidx::is_dotted();
+bool spinidx::is_undotted();
+@end example
+
+allow you to check whether or not a @code{spinidx} object is dotted (use
+@code{ex_to<spinidx>()} to get the object reference from an expression).
+Finally, the two methods
+
+@example
+ex spinidx::toggle_dot();
+ex spinidx::toggle_variance_dot();
+@end example
+
+create a new index with the same value and dimension but opposite dottedness
+and the same or opposite variance.
+
+@subsection Substituting indices
+
+@cindex @code{subs()}
+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 @code{.subs()} method, as it
+is done for symbols (see @ref{Substituting expressions}).
+
+You have two possibilities here. You can either substitute the whole index
+by another index or expression:
+
+@example
+ ...
+ 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
+ ...
+@end example
+
+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 @emph{symbol} of a symbolic index by
+another expression:
+
+@example
+ ...
+ 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
+ ...
+@end example
+
+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:
+
+@example
+ ...
+ ex e = indexed(A, mu_co);
+ cout << e << " becomes " << e.subs(A == A+B) << endl;
+ // -> A.mu becomes (B+A).mu
+ ...
+@end example
+
+@subsection Symmetries
+@cindex @code{symmetry} (class)
+@cindex @code{sy_none()}
+@cindex @code{sy_symm()}
+@cindex @code{sy_anti()}
+@cindex @code{sy_cycl()}
+
+Indexed objects can have certain symmetry properties with respect to their
+indices. Symmetries are specified as a tree of objects of class @code{symmetry}
+that is constructed with the helper functions
+
+@example
+symmetry sy_none(...);
+symmetry sy_symm(...);
+symmetry sy_anti(...);
+symmetry sy_cycl(...);
+@end example
+
+@code{sy_none()} stands for no symmetry, @code{sy_symm()} and @code{sy_anti()}
+specify fully symmetric or antisymmetric, respectively, and @code{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 @code{sy_symm()}, @code{sy_anti()}
+or @code{sy_cycl()} with no arguments specifies the respective symmetry for
+all indices.
+
+Here are some examples of symmetry definitions:
+
+@example
+ ...
+ // 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
+ ...
+@end example
+
+If you need to specify more than four indices, you have to use the
+@code{.add()} method of the @code{symmetry} class. For example, to specify
+full symmetry in the first six indices you would write
+@code{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:
+
+@example
+ ...
+ 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
+ ...
+@end example
+
+@cindex @code{get_free_indices()}
+@cindex dummy index
+@subsection Dummy indices
+
+GiNaC treats certain symbolic index pairs as @dfn{dummy indices} meaning
+that a summation over the index range is implied. Symbolic indices which are
+not dummy indices are called @dfn{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
+@samp{2*n+1} is never a dummy index). If the indices are of class
+@code{varidx} they must also be of opposite variance; if they are of class
+@code{spinidx} they must be both dotted or both undotted.
+
+The method @code{.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:
+
+@example
+@{
+ 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"
+@}
+@end example
+
+@cindex @code{expand_dummy_sum()}
+A dummy index summation like
+@tex
+$ a_i b^i$
+@end tex
+@ifnottex
+a.i b~i
+@end ifnottex
+can be expanded for indices with numeric
+dimensions (e.g. 3) into the explicit sum like
+@tex
+$a_1b^1+a_2b^2+a_3b^3 $.
+@end tex
+@ifnottex
+a.1 b~1 + a.2 b~2 + a.3 b~3.
+@end ifnottex
+This is performed by the function
+
+@example
+ ex expand_dummy_sum(const ex & e, bool subs_idx = false);
+@end example
+
+which takes an expression @code{e} and returns the expanded sum for all
+dummy indices with numeric dimensions. If the parameter @code{subs_idx}
+is set to @code{true} then all substitutions are made by @code{idx} class
+indices, i.e. without variance. In this case the above sum
+@tex
+$ a_i b^i$
+@end tex
+@ifnottex
+a.i b~i
+@end ifnottex
+will be expanded to
+@tex
+$a_1b_1+a_2b_2+a_3b_3 $.
+@end tex
+@ifnottex
+a.1 b.1 + a.2 b.2 + a.3 b.3.
+@end ifnottex
+
+
+@cindex @code{simplify_indexed()}
+@subsection 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
+
+@example
+ex ex::simplify_indexed();
+ex ex::simplify_indexed(const scalar_products & sp);
+@end example
+
+that performs some more expensive operations:
+
+@itemize @bullet
+@item it checks the consistency of free indices in sums in the same way
+ @code{get_free_indices()} does
+@item it tries to give dummy indices that appear in different terms of a sum
+ the same name to allow simplifications like @math{a_i*b_i-a_j*b_j=0}
+@item it (symbolically) calculates all possible dummy index summations/contractions
+ with the predefined tensors (this will be explained in more detail in the
+ next section)
+@item it detects contractions that vanish for symmetry reasons, for example
+ the contraction of a symmetric and a totally antisymmetric tensor
+@item as a special case of dummy index summation, it can replace scalar products
+ of two tensors with a user-defined value
+@end itemize
+
+The last point is done with the help of the @code{scalar_products} class
+which is used to store scalar products with known values (this is not an
+arithmetic class, you just pass it to @code{simplify_indexed()}):
+
+@example
+@{
+ 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
+@}
+@end example
+
+The @code{scalar_products} object @code{sp} acts as a storage for the
+scalar products added to it with the @code{.add()} method. This method
+takes three arguments: the two expressions of which the scalar product is
+taken, and the expression to replace it with.
+
+@cindex @code{expand()}
+The example above also illustrates a feature of the @code{expand()} method:
+if passed the @code{expand_indexed} option it will distribute indices
+over sums, so @samp{(A+B).i} becomes @samp{A.i+B.i}.
+
+@cindex @code{tensor} (class)
+@subsection 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).
+
+@cindex @code{delta_tensor()}
+@subsubsection Delta tensor
+
+The delta tensor takes two indices, is symmetric and has the matrix
+representation @code{diag(1, 1, 1, ...)}. It is constructed by the function
+@code{delta_tensor()}:
+
+@example
+@{
+ 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
+@}
+@end example
+
+@cindex @code{metric_tensor()}
+@subsubsection General metric tensor
+
+The function @code{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 @samp{g} in the output and if its indices are of
+mixed variance it is automatically replaced by a delta tensor:
+
+@example
+@{
+ 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
+@}
+@end example
+
+@cindex @code{lorentz_g()}
+@subsubsection Minkowski metric tensor
+
+The Minkowski metric tensor is a special metric tensor with a constant
+matrix representation which is either @code{diag(1, -1, -1, ...)} (negative
+signature, the default) or @code{diag(-1, 1, 1, ...)} (positive signature).
+It is created with the function @code{lorentz_g()} (although it is output as
+@samp{eta}):
+
+@example
+@{
+ 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
+@}
+@end example
+
+@cindex @code{spinor_metric()}
+@subsubsection Spinor metric tensor
+
+The function @code{spinor_metric()} creates an antisymmetric tensor with
+two indices that is used to raise/lower indices of 2-component spinors.
+It is output as @samp{eps}:
+
+@example
+@{
+ 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
+@}
+@end example
+
+The matrix representation of the spinor metric is @code{[[0, 1], [-1, 0]]}.
+
+@cindex @code{epsilon_tensor()}
+@cindex @code{lorentz_eps()}
+@subsubsection 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 @samp{eps.1.2.3...} (resp. @samp{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
+@samp{eps}.
+
+There are three functions defined to create epsilon tensors in 2, 3 and 4
+dimensions:
+
+@example
+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);
+@end example
+
+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 @code{bool} argument specifies whether the metric
+has negative or positive signature, as in the case of the Minkowski metric
+tensor):
+
+@example
+@{
+ 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
+@}
+@end example
+
+@subsection Linear algebra
+
+The @code{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):
+
+@example
+@{
+ 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
+@}
+@end example
+
+You can of course obtain the same results with the @code{matrix::add()},
+@code{matrix::mul()} and @code{matrix::trace()} methods (@pxref{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 @samp{F~mu~nu} and
+@samp{F.mu.nu} are different matrices. In this case you should use only
+one form for @samp{F} and explicitly multiply it with a matrix representation
+of the metric tensor.
+
+
+@node Non-commutative objects, Hash maps, Indexed objects, Basic concepts
+@c node-name, next, previous, up
+@section 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:
+
+@itemize
+@item Clifford (Dirac) algebra (class @code{clifford})
+@item su(3) Lie algebra (class @code{color})
+@item Matrices (unindexed) (class @code{matrix})
+@end itemize
+
+The @code{clifford} and @code{color} classes are subclasses of
+@code{indexed} because the elements of these algebras usually carry
+indices. The @code{matrix} class is described in more detail in
+@ref{Matrices}.
+
+Unlike most computer algebra systems, GiNaC does not primarily provide an
+operator (often denoted @samp{&*}) 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 @samp{*} 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:
+
+@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)
+ ...
+@end example
+
+As can be seen, GiNaC pulls out the overall commutative factor @samp{-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
+@emph{commutative} product with two factors that are themselves non-commutative
+products (@samp{gamma~mu*gamma~nu} and @samp{T.a*T.b}). For clarification,
+parentheses are placed around the non-commutative products in the output.
+
+@cindex @code{ncmul} (class)
+Non-commutative products are internally represented by objects of the class
+@code{ncmul}, as opposed to commutative products which are handled by the
+@code{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. Both
+symbols and user-defined functions can be specified as being non-commutative.
+
+@cindex @code{return_type()}
+@cindex @code{return_type_tinfo()}
+Information about the commutativity of an object or expression can be
+obtained with the two member functions
+
+@example
+unsigned ex::return_type() const;
+unsigned ex::return_type_tinfo() const;
+@end example
+
+The @code{return_type()} function returns one of three values (defined in
+the header file @file{flags.h}), corresponding to three categories of
+expressions in GiNaC:
+
+@itemize @bullet
+@item @code{return_types::commutative}: Commutates with everything. Most GiNaC
+ classes are of this kind.
+@item @code{return_types::noncommutative}: Non-commutative, belonging to a
+ certain class of non-commutative objects which can be determined with the
+ @code{return_type_tinfo()} method. Expressions of this category commutate
+ with everything except @code{noncommutative} expressions of the same
+ class.
+@item @code{return_types::noncommutative_composite}: Non-commutative, composed
+ of non-commutative objects of different classes. Expressions of this
+ category don't commutate with any other @code{noncommutative} or
+ @code{noncommutative_composite} expressions.
+@end itemize
+
+The value returned by the @code{return_type_tinfo()} method is valid only
+when the return type of the expression is @code{noncommutative}. It is a
+value that is unique to the class of the object, but may vary every time a
+GiNaC program is being run (it is dynamically assigned on start-up).
+
+Here are a couple of examples:
+
+@cartouche
+@multitable @columnfractions 0.33 0.33 0.34
+@item @strong{Expression} @tab @strong{@code{return_type()}} @tab @strong{@code{return_type_tinfo()}}
+@item @code{42} @tab @code{commutative} @tab -
+@item @code{2*x-y} @tab @code{commutative} @tab -
+@item @code{dirac_ONE()} @tab @code{noncommutative} @tab @code{TINFO_clifford}
+@item @code{dirac_gamma(mu)*dirac_gamma(nu)} @tab @code{noncommutative} @tab @code{TINFO_clifford}
+@item @code{2*color_T(a)} @tab @code{noncommutative} @tab @code{TINFO_color}
+@item @code{dirac_ONE()*color_T(a)} @tab @code{noncommutative_composite} @tab -
+@end multitable
+@end cartouche
+
+Note: the @code{return_type_tinfo()} of Clifford objects is only equal to
+@code{TINFO_clifford} for objects with a representation label of zero.
+Other representation labels yield a different @code{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,
+@code{pow(a*b, 2)} becomes @samp{a*b*a*b} if @samp{a} and @samp{b} are
+non-commutative expressions).
+
+
+@cindex @code{clifford} (class)
+@subsection Clifford algebra
+
+
+Clifford algebras are supported in two flavours: Dirac gamma
+matrices (more physical) and generic Clifford algebras (more
+mathematical).
+
+@cindex @code{dirac_gamma()}
+@subsubsection Dirac gamma matrices
+Dirac gamma matrices (note that GiNaC doesn't treat them
+as matrices) are designated as @samp{gamma~mu} and satisfy
+@samp{gamma~mu*gamma~nu + gamma~nu*gamma~mu = 2*eta~mu~nu} where
+@samp{eta~mu~nu} is the Minkowski metric tensor. Dirac gammas are
+constructed by the function
+
+@example
+ex dirac_gamma(const ex & mu, unsigned char rl = 0);
+@end example
+
+which takes two arguments: the index and a @dfn{representation label} in the
+range 0 to 255 which is used to distinguish elements of different Clifford
+algebras (this is also called a @dfn{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.
+
+@cindex @code{dirac_ONE()}
+The unity element of a Clifford algebra is constructed by
+
+@example
+ex dirac_ONE(unsigned char rl = 0);
+@end example
+
+@strong{Please notice:} You must always use @code{dirac_ONE()} when referring to
+multiples of the unity element, even though it's customary to omit it.
+E.g. instead of @code{dirac_gamma(mu)*(dirac_slash(q,4)+m)} you have to
+write @code{dirac_gamma(mu)*(dirac_slash(q,4)+m*dirac_ONE())}. Otherwise,
+GiNaC will complain and/or produce incorrect results.
+
+@cindex @code{dirac_gamma5()}
+There is a special element @samp{gamma5} that commutates with all other
+gammas, has a unit square, and in 4 dimensions equals
+@samp{gamma~0 gamma~1 gamma~2 gamma~3}, provided by
+
+@example
+ex dirac_gamma5(unsigned char rl = 0);
+@end example
+
+@cindex @code{dirac_gammaL()}
+@cindex @code{dirac_gammaR()}
+The chiral projectors @samp{(1+/-gamma5)/2} are also available as proper
+objects, constructed by
+
+@example
+ex dirac_gammaL(unsigned char rl = 0);
+ex dirac_gammaR(unsigned char rl = 0);
+@end example
+
+They observe the relations @samp{gammaL^2 = gammaL}, @samp{gammaR^2 = gammaR},
+and @samp{gammaL gammaR = gammaR gammaL = 0}.
+
+@cindex @code{dirac_slash()}
+Finally, the function
+
+@example
+ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
+@end example
+
+creates a term that represents a contraction of @samp{e} with the Dirac
+Lorentz vector (it behaves like a term of the form @samp{e.mu gamma~mu}
+with a unique index whose dimension is given by the @code{dim} argument).
+Such slashed expressions are printed with a trailing backslash, e.g. @samp{e\}.
+
+In products of dirac gammas, superfluous unity elements are automatically
+removed, squares are replaced by their values, and @samp{gamma5}, @samp{gammaL}
+and @samp{gammaR} are moved to the front.
+
+The @code{simplify_indexed()} function performs contractions in gamma strings,
+for example
+
+@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\
+ ...
+@}
+@end example
+
+@cindex @code{dirac_trace()}
+To calculate the trace of an expression containing strings of Dirac gammas
+you use one of the functions
+
+@example
+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);
+@end example
+
+These functions take the trace over all gammas in the specified set @code{rls}
+or list @code{rll} of representation labels, or the single label @code{rl};
+gammas with other labels are left standing. The last argument to
+@code{dirac_trace()} is the value to be returned for the trace of the unity
+element, which defaults to 4.
+
+The @code{dirac_trace()} function is a linear functional that is equal to the
+ordinary matrix trace only in @math{D = 4} dimensions. In particular, the
+functional is not cyclic in
+@tex $D \ne 4$
+@end tex
+@ifnottex
+@math{D != 4}
+@end ifnottex
+dimensions when acting on
+expressions containing @samp{gamma5}, so it's not a proper trace. This
+@samp{gamma5} scheme is described in greater detail in the article
+@cite{The Role of gamma5 in Dimensional Regularization} (@ref{Bibliography}).
+
+The value of the trace itself is also usually different in 4 and in
+@tex $D \ne 4$
+@end tex
+@ifnottex
+@math{D != 4}
+@end ifnottex
+dimensions:
+
+@example
+@{
+ // 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
+@}
+@end example
+
+Here is an example for using @code{dirac_trace()} to compute a value that
+appears in the calculation of the one-loop vacuum polarization amplitude in
+QED:
+
+@example
+@{
+ 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
+@}
+@end example
+
+The @code{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:
+
+@example
+@{
+ 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
+@}
+@end example
+
+@cindex @code{clifford_unit()}
+@subsubsection A generic Clifford algebra
+
+A generic Clifford algebra, i.e. a
+@tex $2^n$
+@end tex
+@ifnottex
+2^n
+@end ifnottex
+dimensional algebra with
+generators
+@tex $e_k$
+@end tex
+@ifnottex
+e_k
+@end ifnottex
+satisfying the identities
+@tex
+$e_i e_j + e_j e_i = M(i, j) + M(j, i)$
+@end tex
+@ifnottex
+e~i e~j + e~j e~i = M(i, j) + M(j, i)
+@end ifnottex
+for some bilinear form (@code{metric})
+@math{M(i, j)}, which may be non-symmetric (see arXiv:math.QA/9911180)
+and contain symbolic entries. Such generators are created by the
+function
+
+@example
+ ex clifford_unit(const ex & mu, const ex & metr, unsigned char rl = 0);
+@end example
+
+where @code{mu} should be a @code{idx} (or descendant) class object
+indexing the generators.
+Parameter @code{metr} defines the metric @math{M(i, j)} and can be
+represented by a square @code{matrix}, @code{tensormetric} or @code{indexed} class
+object. In fact, any expression either with two free indices or without
+indices at all is admitted as @code{metr}. In the later case an @code{indexed}
+object with two newly created indices with @code{metr} as its
+@code{op(0)} will be used.
+Optional parameter @code{rl} allows to distinguish different
+Clifford algebras, which will commute with each other.
+
+Note that the call @code{clifford_unit(mu, minkmetric())} creates
+something very close to @code{dirac_gamma(mu)}, although
+@code{dirac_gamma} have more efficient simplification mechanism.
+@cindex @code{clifford::get_metric()}
+The method @code{clifford::get_metric()} returns a metric defining this
+Clifford number.
+
+If the matrix @math{M(i, j)} is in fact symmetric you may prefer to create
+the Clifford algebra units with a call like that
+
+@example
+ ex e = clifford_unit(mu, indexed(M, sy_symm(), i, j));
+@end example
+
+since this may yield some further automatic simplifications. Again, for a
+metric defined through a @code{matrix} such a symmetry is detected
+automatically.
+
+Individual generators of a Clifford algebra can be accessed in several
+ways. For example
+
+@example
+@{
+ ...
+ idx i(symbol("i"), 4);
+ realsymbol s("s");
+ ex M = diag_matrix(lst(1, -1, 0, s));
+ ex e = clifford_unit(i, M);
+ ex e0 = e.subs(i == 0);
+ ex e1 = e.subs(i == 1);
+ ex e2 = e.subs(i == 2);
+ ex e3 = e.subs(i == 3);
+ ...
+@}
+@end example
+
+will produce four anti-commuting generators of a Clifford algebra with properties
+@tex
+$e_0^2=1 $, $e_1^2=-1$, $e_2^2=0$ and $e_3^2=s$.
+@end tex
+@ifnottex
+@code{pow(e0, 2) = 1}, @code{pow(e1, 2) = -1}, @code{pow(e2, 2) = 0} and
+@code{pow(e3, 2) = s}.
+@end ifnottex
+
+@cindex @code{lst_to_clifford()}
+A similar effect can be achieved from the function
+
+@example
+ 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);
+@end example
+
+which converts a list or vector
+@tex
+$v = (v^0, v^1, ..., v^n)$
+@end tex
+@ifnottex
+@samp{v = (v~0, v~1, ..., v~n)}
+@end ifnottex
+into the
+Clifford number
+@tex
+$v^0 e_0 + v^1 e_1 + ... + v^n e_n$
+@end tex
+@ifnottex
+@samp{v~0 e.0 + v~1 e.1 + ... + v~n e.n}
+@end ifnottex
+with @samp{e.k}
+directly supplied in the second form of the procedure. In the first form
+the Clifford unit @samp{e.k} is generated by the call of
+@code{clifford_unit(mu, metr, rl)}. The previous code may be rewritten
+with the help of @code{lst_to_clifford()} as follows
+
+@example
+@{
+ ...
+ idx i(symbol("i"), 4);
+ realsymbol s("s");
+ ex M = diag_matrix(lst(1, -1, 0, s));
+ ex e0 = lst_to_clifford(lst(1, 0, 0, 0), i, M);
+ ex e1 = lst_to_clifford(lst(0, 1, 0, 0), i, M);
+ ex e2 = lst_to_clifford(lst(0, 0, 1, 0), i, M);
+ ex e3 = lst_to_clifford(lst(0, 0, 0, 1), i, M);
+ ...
+@}
+@end example
+
+@cindex @code{clifford_to_lst()}
+There is the inverse function
+
+@example
+ lst clifford_to_lst(const ex & e, const ex & c, bool algebraic = true);
+@end example
+
+which takes an expression @code{e} and tries to find a list
+@tex
+$v = (v^0, v^1, ..., v^n)$
+@end tex
+@ifnottex
+@samp{v = (v~0, v~1, ..., v~n)}
+@end ifnottex
+such that
+@tex
+$e = v^0 c_0 + v^1 c_1 + ... + v^n c_n$
+@end tex
+@ifnottex
+@samp{e = v~0 c.0 + v~1 c.1 + ... + v~n c.n}
+@end ifnottex
+with respect to the given Clifford units @code{c} and with none of the
+@samp{v~k} containing Clifford units @code{c} (of course, this
+may be impossible). This function can use an @code{algebraic} method
+(default) or a symbolic one. With the @code{algebraic} method the @samp{v~k} are calculated as
+@tex
+$(e c_k + c_k e)/c_k^2$. If $c_k^2$
+@end tex
+@ifnottex
+@samp{(e c.k + c.k e)/pow(c.k, 2)}. If @samp{pow(c.k, 2)}
+@end ifnottex
+is zero or is not @code{numeric} for some @samp{k}
+then the method will be automatically changed to symbolic. The same effect
+is obtained by the assignment (@code{algebraic = false}) in the procedure call.
+
+@cindex @code{clifford_prime()}
+@cindex @code{clifford_star()}
+@cindex @code{clifford_bar()}
+There are several functions for (anti-)automorphisms of Clifford algebras:
+
+@example
+ 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()); @}
+@end example
+
+The automorphism of a Clifford algebra @code{clifford_prime()} simply
+changes signs of all Clifford units in the expression. The reversion
+of a Clifford algebra @code{clifford_star()} coincides with the
+@code{conjugate()} method and effectively reverses the order of Clifford
+units in any product. Finally the main anti-automorphism
+of a Clifford algebra @code{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
+@math{e'},
+@tex
+$e^*$
+@end tex
+@ifnottex
+e*
+@end ifnottex
+and
+@tex
+$\overline{e}$
+@end tex
+@ifnottex
+@code{\bar@{e@}}
+@end ifnottex
+used in Clifford algebra textbooks.
+
+@cindex @code{clifford_norm()}
+The function
+
+@example
+ ex clifford_norm(const ex & e);
+@end example
+
+@cindex @code{clifford_inverse()}
+calculates the norm of a Clifford number from the expression
+@tex
+$||e||^2 = e\overline{e}$.
+@end tex
+@ifnottex
+@code{||e||^2 = e \bar@{e@}}
+@end ifnottex
+ The inverse of a Clifford expression is returned by the function
+
+@example
+ ex clifford_inverse(const ex & e);
+@end example
+
+which calculates it as
+@tex
+$e^{-1} = \overline{e}/||e||^2$.
+@end tex
+@ifnottex
+@math{e^@{-1@} = \bar@{e@}/||e||^2}
+@end ifnottex
+ If
+@tex
+$||e|| = 0$
+@end tex
+@ifnottex
+@math{||e||=0}
+@end ifnottex
+then an exception is raised.
+
+@cindex @code{remove_dirac_ONE()}
+If a Clifford number happens to be a factor of
+@code{dirac_ONE()} then we can convert it to a ``real'' (non-Clifford)
+expression by the function
+
+@example
+ ex remove_dirac_ONE(const ex & e);
+@end example
+
+@cindex @code{canonicalize_clifford()}
+The function @code{canonicalize_clifford()} works for a
+generic Clifford algebra in a similar way as for Dirac gammas.
+
+The next provided function is
+
+@cindex @code{clifford_moebius_map()}
+@example
+ 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);
+@end example
+
+It takes a list or vector @code{v} and makes the Moebius (conformal or
+linear-fractional) transformation @samp{v -> (av+b)/(cv+d)} defined by
+the matrix @samp{M = [[a, b], [c, d]]}. The parameter @code{G} defines
+the metric of the surrounding (pseudo-)Euclidean space. This can be an
+indexed object, tensormetric, matrix or a Clifford unit, in the later
+case the optional parameter @code{rl} is ignored even if supplied.
+Depending from the type of @code{v} the returned value of this function
+is either a vector or a list holding vector's components.
+
+@cindex @code{clifford_max_label()}
+Finally the function
+
+@example
+char clifford_max_label(const ex & e, bool ignore_ONE = false);
+@end example
+
+can detect a presence of Clifford objects in the expression @code{e}: if
+such objects are found it returns the maximal
+@code{representation_label} of them, otherwise @code{-1}. The optional
+parameter @code{ignore_ONE} indicates if @code{dirac_ONE} objects should
+be ignored during the search.
+
+LaTeX output for Clifford units looks like
+@code{\clifford[1]@{e@}^@{@{\nu@}@}}, where @code{1} is the
+@code{representation_label} and @code{\nu} is the index of the
+corresponding unit. This provides a flexible typesetting with a suitable
+definition of the @code{\clifford} command. For example, the definition
+@example
+ \newcommand@{\clifford@}[1][]@{@}
+@end example
+typesets all Clifford units identically, while the alternative definition
+@example
+ \newcommand@{\clifford@}[2][]@{\ifcase #1 #2\or \tilde@{#2@} \or \breve@{#2@} \fi@}
+@end example
+prints units with @code{representation_label=0} as
+@tex
+$e$,
+@end tex
+@ifnottex
+@code{e},
+@end ifnottex
+with @code{representation_label=1} as
+@tex
+$\tilde{e}$
+@end tex
+@ifnottex
+@code{\tilde@{e@}}
+@end ifnottex
+ and with @code{representation_label=2} as
+@tex
+$\breve{e}$.
+@end tex
+@ifnottex
+@code{\breve@{e@}}.
+@end ifnottex
+
+@cindex @code{color} (class)
+@subsection Color algebra
+
+@cindex @code{color_T()}
+For computations in quantum chromodynamics, GiNaC implements the base elements
+and structure constants of the su(3) Lie algebra (color algebra). The base
+elements @math{T_a} are constructed by the function
+
+@example
+ex color_T(const ex & a, unsigned char rl = 0);
+@end example
+
+which takes two arguments: the index and a @dfn{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 @code{idx},
+not @code{varidx}.
+
+@cindex @code{color_ONE()}
+The unity element of a color algebra is constructed by
+
+@example
+ex color_ONE(unsigned char rl = 0);
+@end example
+
+@strong{Please notice:} You must always use @code{color_ONE()} when referring to
+multiples of the unity element, even though it's customary to omit it.
+E.g. instead of @code{color_T(a)*(color_T(b)*indexed(X,b)+1)} you have to
+write @code{color_T(a)*(color_T(b)*indexed(X,b)+color_ONE())}. Otherwise,
+GiNaC may produce incorrect results.
+
+@cindex @code{color_d()}
+@cindex @code{color_f()}
+The functions
+
+@example
+ex color_d(const ex & a, const ex & b, const ex & c);
+ex color_f(const ex & a, const ex & b, const ex & c);
+@end example
+
+create the symmetric and antisymmetric structure constants @math{d_abc} and
+@math{f_abc} which satisfy @math{@{T_a, T_b@} = 1/3 delta_ab + d_abc T_c}
+and @math{[T_a, T_b] = i f_abc T_c}.
+
+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.
+
+@cindex @code{color_h()}
+There's an additional function
+
+@example
+ex color_h(const ex & a, const ex & b, const ex & c);
+@end example
+
+which returns the linear combination @samp{color_d(a, b, c)+I*color_f(a, b, c)}.
+
+The function @code{simplify_indexed()} performs some simplifications on
+expressions containing color objects:
+
+@example
+@{
+ ...
+ 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
+ ...
+@end example
+
+@cindex @code{color_trace()}
+To calculate the trace of an expression containing color objects you use one
+of the functions
+
+@example
+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);
+@end example
+
+These functions take the trace over all color @samp{T} objects in the
+specified set @code{rls} or list @code{rll} of representation labels, or the
+single label @code{rl}; @samp{T}s with other labels are left standing. For
+example:
+
+@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
+@}
+@end example
+
+
+@node Hash maps, Methods and functions, Non-commutative objects, Basic concepts
+@c node-name, next, previous, up
+@section Hash Maps
+@cindex hash maps
+@cindex @code{exhashmap} (class)
+
+For your convenience, GiNaC offers the container template @code{exhashmap<T>}
+that can be used as a drop-in replacement for the STL
+@code{std::map<ex, T, ex_is_less>}, using hash tables to provide faster,
+typically constant-time, element look-up than @code{map<>}.
+
+@code{exhashmap<>} supports all @code{map<>} members and operations, with the
+following differences:
+
+@itemize @bullet
+@item
+no @code{lower_bound()} and @code{upper_bound()} methods
+@item
+no reverse iterators, no @code{rbegin()}/@code{rend()}
+@item
+no @code{operator<(exhashmap, exhashmap)}
+@item
+the comparison function object @code{key_compare} is hardcoded to
+@code{ex_is_less}
+@item
+the constructor @code{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)
+@item
+the method @code{size_t bucket_count()} returns the current size of the hash
+table
+@item
+@code{insert()} and @code{erase()} operations invalidate all iterators
+@end itemize
+
+
+@node Methods and functions, Information about expressions, Hash maps, Top
+@c node-name, next, previous, up
+@chapter Methods and functions
+@cindex polynomial
+
+In this chapter the most important algorithms provided by GiNaC will be
+described. Some of them are implemented as functions on expressions,
+others are implemented as methods provided by expression objects. If
+they are methods, there exists a wrapper function around it, so you can
+alternatively call it in a functional way as shown in the simple
+example:
+
+@example
+ ...
+ cout << "As method: " << sin(1).evalf() << endl;
+ cout << "As function: " << evalf(sin(1)) << endl;
+ ...
+@end example
+
+@cindex @code{subs()}
+The general rule is that wherever methods accept one or more parameters
+(@var{arg1}, @var{arg2}, @dots{}) the order of arguments the function
+wrapper accepts is the same but preceded by the object to act on
+(@var{object}, @var{arg1}, @var{arg2}, @dots{}). This approach is the
+most natural one in an OO model but it may lead to confusion for MapleV
+users because where they would type @code{A:=x+1; subs(x=2,A);} GiNaC
+would require @code{A=x+1; subs(A,x==2);} (after proper declaration of
+@code{A} and @code{x}). On the other hand, since MapleV returns 3 on
+@code{A:=x^2+3; coeff(A,x,0);} (GiNaC: @code{A=pow(x,2)+3;
+coeff(A,x,0);}) it is clear that MapleV is not trying to be consistent
+here. Also, users of MuPAD will in most cases feel more comfortable
+with GiNaC's convention. All function wrappers are implemented
+as simple inline functions which just call the corresponding method and
+are only provided for users uncomfortable with OO who are dead set to
+avoid method invocations. Generally, nested function wrappers are much
+harder to read than a sequence of methods and should therefore be
+avoided if possible. On the other hand, not everything in GiNaC is a
+method on class @code{ex} and sometimes calling a function cannot be
+avoided.
+
+@menu
+* Information about expressions::
+* Numerical evaluation::
+* Substituting expressions::
+* Pattern matching and advanced substitutions::
+* Applying a function on subexpressions::
+* Visitors and tree traversal::
+* Polynomial arithmetic:: Working with polynomials.
+* Rational expressions:: Working with rational functions.
+* Symbolic differentiation::
+* Series expansion:: Taylor and Laurent expansion.
+* Symmetrization::
+* Built-in functions:: List of predefined mathematical functions.
+* Multiple polylogarithms::
+* Complex expressions::
+* Solving linear systems of equations::
+* Input/output:: Input and output of expressions.
+@end menu
+
+
+@node Information about expressions, Numerical evaluation, Methods and functions, Methods and functions
+@c node-name, next, previous, up
+@section Getting information about expressions
+
+@subsection Checking expression types
+@cindex @code{is_a<@dots{}>()}
+@cindex @code{is_exactly_a<@dots{}>()}
+@cindex @code{ex_to<@dots{}>()}
+@cindex Converting @code{ex} to other classes
+@cindex @code{info()}
+@cindex @code{return_type()}
+@cindex @code{return_type_tinfo()}
+
+Sometimes it's useful to check whether a given expression is a plain number,
+a sum, a polynomial with integer coefficients, or of some other specific type.
+GiNaC provides a couple of functions for this:
+
+@example
+bool is_a<T>(const ex & e);
+bool is_exactly_a<T>(const ex & e);
+bool ex::info(unsigned flag);
+unsigned ex::return_type() const;
+unsigned ex::return_type_tinfo() const;
+@end example
+
+When the test made by @code{is_a<T>()} returns true, it is safe to call
+one of the functions @code{ex_to<T>()}, where @code{T} is one of the
+class names (@xref{The class hierarchy}, for a list of all classes). For
+example, assuming @code{e} is an @code{ex}:
+
+@example
+@{
+ @dots{}
+ if (is_a<numeric>(e))
+ numeric n = ex_to<numeric>(e);
+ @dots{}
+@}
+@end example
+
+@code{is_a<T>(e)} allows you to check whether the top-level object of
+an expression @samp{e} is an instance of the GiNaC class @samp{T}
+(@xref{The class hierarchy}, for a list of all classes). This is most useful,
+e.g., for checking whether an expression is a number, a sum, or a product:
+
+@example
+@{
+ symbol x("x");
+ ex e1 = 42;
+ ex e2 = 4*x - 3;
+ is_a<numeric>(e1); // true
+ is_a<numeric>(e2); // false
+ is_a<add>(e1); // false
+ is_a<add>(e2); // true
+ is_a<mul>(e1); // false
+ is_a<mul>(e2); // false
+@}
+@end example
+
+In contrast, @code{is_exactly_a<T>(e)} allows you to check whether the
+top-level object of an expression @samp{e} is an instance of the GiNaC
+class @samp{T}, not including parent classes.
+
+The @code{info()} method is used for checking certain attributes of
+expressions. The possible values for the @code{flag} argument are defined
+in @file{ginac/flags.h}, the most important being explained in the following
+table:
+
+@cartouche
+@multitable @columnfractions .30 .70
+@item @strong{Flag} @tab @strong{Returns true if the object is@dots{}}
+@item @code{numeric}
+@tab @dots{}a number (same as @code{is_a<numeric>(...)})
+@item @code{real}
+@tab @dots{}a real number, symbol or constant (i.e. is not complex)
+@item @code{rational}
+@tab @dots{}an exact rational number (integers are rational, too)
+@item @code{integer}
+@tab @dots{}a (non-complex) integer
+@item @code{crational}
+@tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
+@item @code{cinteger}
+@tab @dots{}a (complex) integer (such as @math{2-3*I})
+@item @code{positive}
+@tab @dots{}not complex and greater than 0
+@item @code{negative}
+@tab @dots{}not complex and less than 0
+@item @code{nonnegative}
+@tab @dots{}not complex and greater than or equal to 0
+@item @code{posint}
+@tab @dots{}an integer greater than 0
+@item @code{negint}
+@tab @dots{}an integer less than 0
+@item @code{nonnegint}
+@tab @dots{}an integer greater than or equal to 0
+@item @code{even}
+@tab @dots{}an even integer
+@item @code{odd}
+@tab @dots{}an odd integer
+@item @code{prime}
+@tab @dots{}a prime integer (probabilistic primality test)
+@item @code{relation}
+@tab @dots{}a relation (same as @code{is_a<relational>(...)})
+@item @code{relation_equal}
+@tab @dots{}a @code{==} relation
+@item @code{relation_not_equal}
+@tab @dots{}a @code{!=} relation
+@item @code{relation_less}
+@tab @dots{}a @code{<} relation
+@item @code{relation_less_or_equal}
+@tab @dots{}a @code{<=} relation
+@item @code{relation_greater}
+@tab @dots{}a @code{>} relation
+@item @code{relation_greater_or_equal}
+@tab @dots{}a @code{>=} relation
+@item @code{symbol}
+@tab @dots{}a symbol (same as @code{is_a<symbol>(...)})
+@item @code{list}
+@tab @dots{}a list (same as @code{is_a<lst>(...)})
+@item @code{polynomial}
+@tab @dots{}a polynomial (i.e. only consists of sums and products of numbers and symbols with positive integer powers)
+@item @code{integer_polynomial}
+@tab @dots{}a polynomial with (non-complex) integer coefficients
+@item @code{cinteger_polynomial}
+@tab @dots{}a polynomial with (possibly complex) integer coefficients (such as @math{2-3*I})
+@item @code{rational_polynomial}
+@tab @dots{}a polynomial with (non-complex) rational coefficients
+@item @code{crational_polynomial}
+@tab @dots{}a polynomial with (possibly complex) rational coefficients (such as @math{2/3+7/2*I})
+@item @code{rational_function}
+@tab @dots{}a rational function (@math{x+y}, @math{z/(x+y)})
+@item @code{algebraic}
+@tab @dots{}an algebraic object (@math{sqrt(2)}, @math{sqrt(x)-1})
+@end multitable
+@end cartouche
+
+To determine whether an expression is commutative or non-commutative and if
+so, with which other expressions it would commutate, you use the methods
+@code{return_type()} and @code{return_type_tinfo()}. @xref{Non-commutative objects},
+for an explanation of these.
+
+
+@subsection Accessing subexpressions
+@cindex container
+
+Many GiNaC classes, like @code{add}, @code{mul}, @code{lst}, and
+@code{function}, act as containers for subexpressions. For example, the
+subexpressions of a sum (an @code{add} object) are the individual terms,
+and the subexpressions of a @code{function} are the function's arguments.
+
+@cindex @code{nops()}
+@cindex @code{op()}
+GiNaC provides several ways of accessing subexpressions. The first way is to
+use the two methods
+
+@example
+size_t ex::nops();
+ex ex::op(size_t i);
+@end example
+
+@code{nops()} determines the number of subexpressions (operands) contained
+in the expression, while @code{op(i)} returns the @code{i}-th
+(0..@code{nops()-1}) subexpression. In the case of a @code{power} object,
+@code{op(0)} will return the basis and @code{op(1)} the exponent. For
+@code{indexed} objects, @code{op(0)} is the base expression and @code{op(i)},
+@math{i>0} are the indices.
+
+@cindex iterators
+@cindex @code{const_iterator}
+The second way to access subexpressions is via the STL-style random-access
+iterator class @code{const_iterator} and the methods
+
+@example
+const_iterator ex::begin();
+const_iterator ex::end();
+@end example
+
+@code{begin()} returns an iterator referring to the first subexpression;
+@code{end()} returns an iterator which is one-past the last subexpression.
+If the expression has no subexpressions, then @code{begin() == end()}. These
+iterators can also be used in conjunction with non-modifying STL algorithms.
+
+Here is an example that (non-recursively) prints the subexpressions of a
+given expression in three different ways:
+
+@example
+@{
+ ex e = ...
+
+ // with nops()/op()
+ for (size_t i = 0; i != e.nops(); ++i)
+ cout << e.op(i) << endl;
+
+ // with iterators
+ for (const_iterator i = e.begin(); i != e.end(); ++i)
+ cout << *i << endl;
+
+ // with iterators and STL copy()
+ std::copy(e.begin(), e.end(), std::ostream_iterator<ex>(cout, "\n"));
+@}
+@end example
+
+@cindex @code{const_preorder_iterator}
+@cindex @code{const_postorder_iterator}
+@code{op()}/@code{nops()} and @code{const_iterator} only access an
+expression's immediate children. GiNaC provides two additional iterator
+classes, @code{const_preorder_iterator} and @code{const_postorder_iterator},
+that iterate over all objects in an expression tree, in preorder or postorder,
+respectively. They are STL-style forward iterators, and are created with the
+methods
+
+@example
+const_preorder_iterator ex::preorder_begin();
+const_preorder_iterator ex::preorder_end();
+const_postorder_iterator ex::postorder_begin();
+const_postorder_iterator ex::postorder_end();
+@end example
+
+The following example illustrates the differences between
+@code{const_iterator}, @code{const_preorder_iterator}, and
+@code{const_postorder_iterator}:
+
+@example
+@{
+ symbol A("A"), B("B"), C("C");
+ ex e = lst(lst(A, B), C);
+
+ std::copy(e.begin(), e.end(),
+ std::ostream_iterator<ex>(cout, "\n"));
+ // @{A,B@}
+ // C
+
+ std::copy(e.preorder_begin(), e.preorder_end(),
+ std::ostream_iterator<ex>(cout, "\n"));
+ // @{@{A,B@},C@}
+ // @{A,B@}
+ // A
+ // B
+ // C
+
+ std::copy(e.postorder_begin(), e.postorder_end(),
+ std::ostream_iterator<ex>(cout, "\n"));
+ // A
+ // B
+ // @{A,B@}
+ // C
+ // @{@{A,B@},C@}
+@}
+@end example
+
+@cindex @code{relational} (class)
+Finally, the left-hand side and right-hand side expressions of objects of
+class @code{relational} (and only of these) can also be accessed with the
+methods
+
+@example
+ex ex::lhs();
+ex ex::rhs();
+@end example
+
+
+@subsection Comparing expressions
+@cindex @code{is_equal()}
+@cindex @code{is_zero()}
+
+Expressions can be compared with the usual C++ relational operators like
+@code{==}, @code{>}, and @code{<} but if the expressions contain symbols,
+the result is usually not determinable and the result will be @code{false},
+except in the case of the @code{!=} operator. You should also be aware that
+GiNaC will only do the most trivial test for equality (subtracting both
+expressions), so something like @code{(pow(x,2)+x)/x==x+1} will return
+@code{false}.
+
+Actually, if you construct an expression like @code{a == b}, this will be
+represented by an object of the @code{relational} class (@pxref{Relations})
+which is not evaluated until (explicitly or implicitly) cast to a @code{bool}.
+
+There are also two methods
+
+@example
+bool ex::is_equal(const ex & other);
+bool ex::is_zero();
+@end example
+
+for checking whether one expression is equal to another, or equal to zero,
+respectively. See also the method @code{ex::is_zero_matrix()},
+@pxref{Matrices}.
+
+
+@subsection Ordering expressions
+@cindex @code{ex_is_less} (class)
+@cindex @code{ex_is_equal} (class)
+@cindex @code{compare()}
+
+Sometimes it is necessary to establish a mathematically well-defined ordering
+on a set of arbitrary expressions, for example to use expressions as keys
+in a @code{std::map<>} container, or to bring a vector of expressions into
+a canonical order (which is done internally by GiNaC for sums and products).
+
+The operators @code{<}, @code{>} etc. described in the last section cannot
+be used for this, as they don't implement an ordering relation in the
+mathematical sense. In particular, they are not guaranteed to be
+antisymmetric: if @samp{a} and @samp{b} are different expressions, and
+@code{a < b} yields @code{false}, then @code{b < a} doesn't necessarily
+yield @code{true}.
+
+By default, STL classes and algorithms use the @code{<} and @code{==}
+operators to compare objects, which are unsuitable for expressions, but GiNaC
+provides two functors that can be supplied as proper binary comparison
+predicates to the STL:
+
+@example
+class ex_is_less : public std::binary_function<ex, ex, bool> @{
+public:
+ bool operator()(const ex &lh, const ex &rh) const;
+@};
+
+class ex_is_equal : public std::binary_function<ex, ex, bool> @{
+public:
+ bool operator()(const ex &lh, const ex &rh) const;
+@};
+@end example
+
+For example, to define a @code{map} that maps expressions to strings you
+have to use
+
+@example
+std::map<ex, std::string, ex_is_less> myMap;
+@end example
+
+Omitting the @code{ex_is_less} template parameter will introduce spurious
+bugs because the map operates improperly.
+
+Other examples for the use of the functors:
+
+@example
+std::vector<ex> v;
+// fill vector
+...
+
+// sort vector
+std::sort(v.begin(), v.end(), ex_is_less());
+
+// count the number of expressions equal to '1'
+unsigned num_ones = std::count_if(v.begin(), v.end(),
+ std::bind2nd(ex_is_equal(), 1));
+@end example
+
+The implementation of @code{ex_is_less} uses the member function
+
+@example
+int ex::compare(const ex & other) const;
+@end example
+
+which returns @math{0} if @code{*this} and @code{other} are equal, @math{-1}
+if @code{*this} sorts before @code{other}, and @math{1} if @code{*this} sorts
+after @code{other}.
+
+
+@node Numerical evaluation, Substituting expressions, Information about expressions, Methods and functions
+@c node-name, next, previous, up
+@section Numerical evaluation
+@cindex @code{evalf()}
+
+GiNaC keeps algebraic expressions, numbers and constants in their exact form.
+To evaluate them using floating-point arithmetic you need to call
+
+@example
+ex ex::evalf(int level = 0) const;
+@end example
+
+@cindex @code{Digits}
+The accuracy of the evaluation is controlled by the global object @code{Digits}
+which can be assigned an integer value. The default value of @code{Digits}
+is 17. @xref{Numbers}, for more information and examples.
+
+To evaluate an expression to a @code{double} floating-point number you can
+call @code{evalf()} followed by @code{numeric::to_double()}, like this:
+
+@example
+@{
+ // Approximate sin(x/Pi)
+ symbol x("x");
+ ex e = series(sin(x/Pi), x == 0, 6);
+
+ // Evaluate numerically at x=0.1
+ ex f = evalf(e.subs(x == 0.1));
+
+ // ex_to<numeric> is an unsafe cast, so check the type first
+ if (is_a<numeric>(f)) @{
+ double d = ex_to<numeric>(f).to_double();
+ cout << d << endl;
+ // -> 0.0318256
+ @} else
+ // error
+@}
+@end example
+
+
+@node Substituting expressions, Pattern matching and advanced substitutions, Numerical evaluation, Methods and functions
+@c node-name, next, previous, up
+@section Substituting expressions
+@cindex @code{subs()}
+
+Algebraic objects inside expressions can be replaced with arbitrary
+expressions via the @code{.subs()} method:
+
+@example
+ex ex::subs(const ex & e, unsigned options = 0);
+ex ex::subs(const exmap & m, unsigned options = 0);
+ex ex::subs(const lst & syms, const lst & repls, unsigned options = 0);
+@end example
+
+In the first form, @code{subs()} accepts a relational of the form
+@samp{object == expression} or a @code{lst} of such relationals:
+
+@example
+@{
+ symbol x("x"), y("y");
+
+ ex e1 = 2*x^2-4*x+3;
+ cout << "e1(7) = " << e1.subs(x == 7) << endl;
+ // -> 73
+
+ ex e2 = x*y + x;
+ cout << "e2(-2, 4) = " << e2.subs(lst(x == -2, y == 4)) << endl;
+ // -> -10
+@}
+@end example
+
+If you specify multiple substitutions, they are performed in parallel, so e.g.
+@code{subs(lst(x == y, y == x))} exchanges @samp{x} and @samp{y}.
+
+The second form of @code{subs()} takes an @code{exmap} object which is a
+pair associative container that maps expressions to expressions (currently
+implemented as a @code{std::map}). This is the most efficient one of the
+three @code{subs()} forms and should be used when the number of objects to
+be substituted is large or unknown.
+
+Using this form, the second example from above would look like this:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex e2 = x*y + x;
+
+ exmap m;
+ m[x] = -2;
+ m[y] = 4;
+ cout << "e2(-2, 4) = " << e2.subs(m) << endl;
+@}
+@end example
+
+The third form of @code{subs()} takes two lists, one for the objects to be
+replaced and one for the expressions to be substituted (both lists must
+contain the same number of elements). Using this form, you would write
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex e2 = x*y + x;
+
+ cout << "e2(-2, 4) = " << e2.subs(lst(x, y), lst(-2, 4)) << endl;
+@}
+@end example
+
+The optional last argument to @code{subs()} is a combination of
+@code{subs_options} flags. There are three options available:
+@code{subs_options::no_pattern} disables pattern matching, which makes
+large @code{subs()} operations significantly faster if you are not using
+patterns. The second option, @code{subs_options::algebraic} enables
+algebraic substitutions in products and powers.
+@xref{Pattern matching and advanced substitutions}, for more information
+about patterns and algebraic substitutions. The third option,
+@code{subs_options::no_index_renaming} disables the feature that dummy
+indices are renamed if the substitution could give a result in which a
+dummy index occurs more than two times. This is sometimes necessary if
+you want to use @code{subs()} to rename your dummy indices.
+
+@code{subs()} performs syntactic substitution of any complete algebraic
+object; it does not try to match sub-expressions as is demonstrated by the
+following example:
+
+@example
+@{
+ symbol x("x"), y("y"), z("z");
+
+ ex e1 = pow(x+y, 2);
+ cout << e1.subs(x+y == 4) << endl;
+ // -> 16
+
+ ex e2 = sin(x)*sin(y)*cos(x);
+ cout << e2.subs(sin(x) == cos(x)) << endl;
+ // -> cos(x)^2*sin(y)
+
+ ex e3 = x+y+z;
+ cout << e3.subs(x+y == 4) << endl;
+ // -> x+y+z
+ // (and not 4+z as one might expect)
+@}
+@end example
+
+A more powerful form of substitution using wildcards is described in the
+next section.
+
+
+@node Pattern matching and advanced substitutions, Applying a function on subexpressions, Substituting expressions, Methods and functions
+@c node-name, next, previous, up
+@section Pattern matching and advanced substitutions
+@cindex @code{wildcard} (class)
+@cindex Pattern matching
+
+GiNaC allows the use of patterns for checking whether an expression is of a
+certain form or contains subexpressions of a certain form, and for
+substituting expressions in a more general way.
+
+A @dfn{pattern} is an algebraic expression that optionally contains wildcards.
+A @dfn{wildcard} is a special kind of object (of class @code{wildcard}) that
+represents an arbitrary expression. Every wildcard has a @dfn{label} which is
+an unsigned integer number to allow having multiple different wildcards in a
+pattern. Wildcards are printed as @samp{$label} (this is also the way they
+are specified in @command{ginsh}). In C++ code, wildcard objects are created
+with the call
+
+@example
+ex wild(unsigned label = 0);
+@end example
+
+which is simply a wrapper for the @code{wildcard()} constructor with a shorter
+name.
+
+Some examples for patterns:
+
+@multitable @columnfractions .5 .5
+@item @strong{Constructed as} @tab @strong{Output as}
+@item @code{wild()} @tab @samp{$0}
+@item @code{pow(x,wild())} @tab @samp{x^$0}
+@item @code{atan2(wild(1),wild(2))} @tab @samp{atan2($1,$2)}
+@item @code{indexed(A,idx(wild(),3))} @tab @samp{A.$0}
+@end multitable
+
+Notes:
+
+@itemize @bullet
+@item Wildcards behave like symbols and are subject to the same algebraic
+ rules. E.g., @samp{$0+2*$0} is automatically transformed to @samp{3*$0}.
+@item As shown in the last example, to use wildcards for indices you have to
+ use them as the value of an @code{idx} object. This is because indices must
+ always be of class @code{idx} (or a subclass).
+@item Wildcards only represent expressions or subexpressions. It is not
+ possible to use them as placeholders for other properties like index
+ dimension or variance, representation labels, symmetry of indexed objects
+ etc.
+@item Because wildcards are commutative, it is not possible to use wildcards
+ as part of noncommutative products.
+@item A pattern does not have to contain wildcards. @samp{x} and @samp{x+y}
+ are also valid patterns.
+@end itemize
+
+@subsection Matching expressions
+@cindex @code{match()}
+The most basic application of patterns is to check whether an expression
+matches a given pattern. This is done by the function
+
+@example
+bool ex::match(const ex & pattern);
+bool ex::match(const ex & pattern, lst & repls);
+@end example
+
+This function returns @code{true} when the expression matches the pattern
+and @code{false} if it doesn't. If used in the second form, the actual
+subexpressions matched by the wildcards get returned in the @code{repls}
+object as a list of relations of the form @samp{wildcard == expression}.
+If @code{match()} returns false, the state of @code{repls} is undefined.
+For reproducible results, the list should be empty when passed to
+@code{match()}, but it is also possible to find similarities in multiple
+expressions by passing in the result of a previous match.
+
+The matching algorithm works as follows:
+
+@itemize
+@item A single wildcard matches any expression. If one wildcard appears
+ multiple times in a pattern, it must match the same expression in all
+ places (e.g. @samp{$0} matches anything, and @samp{$0*($0+1)} matches
+ @samp{x*(x+1)} but not @samp{x*(y+1)}).
+@item If the expression is not of the same class as the pattern, the match
+ fails (i.e. a sum only matches a sum, a function only matches a function,
+ etc.).
+@item If the pattern is a function, it only matches the same function
+ (i.e. @samp{sin($0)} matches @samp{sin(x)} but doesn't match @samp{exp(x)}).
+@item Except for sums and products, the match fails if the number of
+ subexpressions (@code{nops()}) is not equal to the number of subexpressions
+ of the pattern.
+@item If there are no subexpressions, the expressions and the pattern must
+ be equal (in the sense of @code{is_equal()}).
+@item Except for sums and products, each subexpression (@code{op()}) must
+ match the corresponding subexpression of the pattern.
+@end itemize
+
+Sums (@code{add}) and products (@code{mul}) are treated in a special way to
+account for their commutativity and associativity:
+
+@itemize
+@item If the pattern contains a term or factor that is a single wildcard,
+ this one is used as the @dfn{global wildcard}. If there is more than one
+ such wildcard, one of them is chosen as the global wildcard in a random
+ way.
+@item Every term/factor of the pattern, except the global wildcard, is
+ matched against every term of the expression in sequence. If no match is
+ found, the whole match fails. Terms that did match are not considered in
+ further matches.
+@item If there are no unmatched terms left, the match succeeds. Otherwise
+ the match fails unless there is a global wildcard in the pattern, in
+ which case this wildcard matches the remaining terms.
+@end itemize
+
+In general, having more than one single wildcard as a term of a sum or a
+factor of a product (such as @samp{a+$0+$1}) will lead to unpredictable or
+ambiguous results.
+
+Here are some examples in @command{ginsh} to demonstrate how it works (the
+@code{match()} function in @command{ginsh} returns @samp{FAIL} if the
+match fails, and the list of wildcard replacements otherwise):
+
+@example
+> match((x+y)^a,(x+y)^a);
+@{@}
+> match((x+y)^a,(x+y)^b);
+FAIL
+> match((x+y)^a,$1^$2);
+@{$1==x+y,$2==a@}
+> match((x+y)^a,$1^$1);
+FAIL
+> match((x+y)^(x+y),$1^$1);
+@{$1==x+y@}
+> match((x+y)^(x+y),$1^$2);
+@{$1==x+y,$2==x+y@}
+> match((a+b)*(a+c),($1+b)*($1+c));
+@{$1==a@}
+> match((a+b)*(a+c),(a+$1)*(a+$2));
+@{$1==b,$2==c@}
+ (Unpredictable. The result might also be [$1==c,$2==b].)
+> match((a+b)*(a+c),($1+$2)*($1+$3));
+ (The result is undefined. Due to the sequential nature of the algorithm
+ and the re-ordering of terms in GiNaC, the match for the first factor
+ may be @{$1==a,$2==b@} in which case the match for the second factor
+ succeeds, or it may be @{$1==b,$2==a@} which causes the second match to
+ fail.)
+> match(a*(x+y)+a*z+b,a*$1+$2);
+ (This is also ambiguous and may return either @{$1==z,$2==a*(x+y)+b@} or
+ @{$1=x+y,$2=a*z+b@}.)
+> match(a+b+c+d+e+f,c);
+FAIL
+> match(a+b+c+d+e+f,c+$0);
+@{$0==a+e+b+f+d@}
+> match(a+b+c+d+e+f,c+e+$0);
+@{$0==a+b+f+d@}
+> match(a+b,a+b+$0);
+@{$0==0@}
+> match(a*b^2,a^$1*b^$2);
+FAIL
+ (The matching is syntactic, not algebraic, and "a" doesn't match "a^$1"
+ even though a==a^1.)
+> match(x*atan2(x,x^2),$0*atan2($0,$0^2));
+@{$0==x@}
+> match(atan2(y,x^2),atan2(y,$0));
+@{$0==x^2@}
+@end example
+
+@subsection Matching parts of expressions
+@cindex @code{has()}
+A more general way to look for patterns in expressions is provided by the
+member function
+
+@example
+bool ex::has(const ex & pattern);
+@end example
+
+This function checks whether a pattern is matched by an expression itself or
+by any of its subexpressions.
+
+Again some examples in @command{ginsh} for illustration (in @command{ginsh},
+@code{has()} returns @samp{1} for @code{true} and @samp{0} for @code{false}):
+
+@example
+> has(x*sin(x+y+2*a),y);
+1
+> has(x*sin(x+y+2*a),x+y);
+0
+ (This is because in GiNaC, "x+y" is not a subexpression of "x+y+2*a" (which
+ has the subexpressions "x", "y" and "2*a".)
+> has(x*sin(x+y+2*a),x+y+$1);
+1
+ (But this is possible.)
+> has(x*sin(2*(x+y)+2*a),x+y);
+0
+ (This fails because "2*(x+y)" automatically gets converted to "2*x+2*y" of
+ which "x+y" is not a subexpression.)
+> has(x+1,x^$1);
+0
+ (Although x^1==x and x^0==1, neither "x" nor "1" are actually of the form
+ "x^something".)
+> has(4*x^2-x+3,$1*x);
+1
+> has(4*x^2+x+3,$1*x);
+0
+ (Another possible pitfall. The first expression matches because the term
+ "-x" has the form "(-1)*x" in GiNaC. To check whether a polynomial
+ contains a linear term you should use the coeff() function instead.)
+@end example
+
+@cindex @code{find()}
+The method
+
+@example
+bool ex::find(const ex & pattern, lst & found);
+@end example
+
+works a bit like @code{has()} but it doesn't stop upon finding the first
+match. Instead, it appends all found matches to the specified list. If there
+are multiple occurrences of the same expression, it is entered only once to
+the list. @code{find()} returns false if no matches were found (in
+@command{ginsh}, it returns an empty list):
+
+@example
+> find(1+x+x^2+x^3,x);
+@{x@}
+> find(1+x+x^2+x^3,y);
+@{@}
+> find(1+x+x^2+x^3,x^$1);
+@{x^3,x^2@}
+ (Note the absence of "x".)
+> expand((sin(x)+sin(y))*(a+b));
+sin(y)*a+sin(x)*b+sin(x)*a+sin(y)*b
+> find(%,sin($1));
+@{sin(y),sin(x)@}
+@end example
+
+@subsection Substituting expressions
+@cindex @code{subs()}
+Probably the most useful application of patterns is to use them for
+substituting expressions with the @code{subs()} method. Wildcards can be
+used in the search patterns as well as in the replacement expressions, where
+they get replaced by the expressions matched by them. @code{subs()} doesn't
+know anything about algebra; it performs purely syntactic substitutions.
+
+Some examples:
+
+@example
+> subs(a^2+b^2+(x+y)^2,$1^2==$1^3);
+b^3+a^3+(x+y)^3
+> subs(a^4+b^4+(x+y)^4,$1^2==$1^3);
+b^4+a^4+(x+y)^4
+> subs((a+b+c)^2,a+b==x);
+(a+b+c)^2
+> subs((a+b+c)^2,a+b+$1==x+$1);
+(x+c)^2
+> subs(a+2*b,a+b==x);
+a+2*b
+> subs(4*x^3-2*x^2+5*x-1,x==a);
+-1+5*a-2*a^2+4*a^3
+> subs(4*x^3-2*x^2+5*x-1,x^$0==a^$0);
+-1+5*x-2*a^2+4*a^3
+> subs(sin(1+sin(x)),sin($1)==cos($1));
+cos(1+cos(x))
+> expand(subs(a*sin(x+y)^2+a*cos(x+y)^2+b,cos($1)^2==1-sin($1)^2));
+a+b
+@end example
+
+The last example would be written in C++ in this way:
+
+@example
+@{
+ symbol a("a"), b("b"), x("x"), y("y");
+ e = a*pow(sin(x+y), 2) + a*pow(cos(x+y), 2) + b;
+ e = e.subs(pow(cos(wild()), 2) == 1-pow(sin(wild()), 2));
+ cout << e.expand() << endl;
+ // -> a+b
+@}
+@end example
+
+@subsection The option algebraic
+Both @code{has()} and @code{subs()} take an optional argument to pass them
+extra options. This section describes what happens if you give the former
+the option @code{has_options::algebraic} or the latter
+@code{subs_options::algebraic}. In that case the matching condition for
+powers and multiplications is changed in such a way that they become
+more intuitive. Intuition says that @code{x*y} is a part of @code{x*y*z}.
+If you use these options you will find that
+@code{(x*y*z).has(x*y, has_options::algebraic)} indeed returns true.
+Besides matching some of the factors of a product also powers match as
+often as is possible without getting negative exponents. For example
+@code{(x^5*y^2*z).subs(x^2*y^2==c, subs_options::algebraic)} will return
+@code{x*c^2*z}. This also works with negative powers:
+@code{(x^(-3)*y^(-2)*z).subs(1/(x*y)==c, subs_options::algebraic)} will
+return @code{x^(-1)*c^2*z}.
+
+@strong{Note:} this only works for multiplications
+and not for locating @code{x+y} within @code{x+y+z}.
+
+
+@node Applying a function on subexpressions, Visitors and tree traversal, Pattern matching and advanced substitutions, Methods and functions
+@c node-name, next, previous, up
+@section Applying a function on subexpressions
+@cindex tree traversal
+@cindex @code{map()}
+
+Sometimes you may want to perform an operation on specific parts of an
+expression while leaving the general structure of it intact. An example
+of this would be a matrix trace operation: the trace of a sum is the sum
+of the traces of the individual terms. That is, the trace should @dfn{map}
+on the sum, by applying itself to each of the sum's operands. It is possible
+to do this manually which usually results in code like this:
+
+@example
+ex calc_trace(ex e)
+@{
+ if (is_a<matrix>(e))
+ return ex_to<matrix>(e).trace();
+ else if (is_a<add>(e)) @{
+ ex sum = 0;
+ for (size_t i=0; i<e.nops(); i++)
+ sum += calc_trace(e.op(i));
+ return sum;
+ @} else if (is_a<mul>)(e)) @{
+ ...
+ @} else @{
+ ...
+ @}
+@}
+@end example
+
+This is, however, slightly inefficient (if the sum is very large it can take
+a long time to add the terms one-by-one), and its applicability is limited to
+a rather small class of expressions. If @code{calc_trace()} is called with
+a relation or a list as its argument, you will probably want the trace to
+be taken on both sides of the relation or of all elements of the list.
+
+GiNaC offers the @code{map()} method to aid in the implementation of such
+operations:
+
+@example
+ex ex::map(map_function & f) const;
+ex ex::map(ex (*f)(const ex & e)) const;
+@end example
+
+In the first (preferred) form, @code{map()} takes a function object that
+is subclassed from the @code{map_function} class. In the second form, it
+takes a pointer to a function that accepts and returns an expression.
+@code{map()} constructs a new expression of the same type, applying the
+specified function on all subexpressions (in the sense of @code{op()}),
+non-recursively.
+
+The use of a function object makes it possible to supply more arguments to
+the function that is being mapped, or to keep local state information.
+The @code{map_function} class declares a virtual function call operator
+that you can overload. Here is a sample implementation of @code{calc_trace()}
+that uses @code{map()} in a recursive fashion:
+
+@example
+struct calc_trace : public map_function @{
+ ex operator()(const ex &e)
+ @{
+ if (is_a<matrix>(e))
+ return ex_to<matrix>(e).trace();
+ else if (is_a<mul>(e)) @{
+ ...
+ @} else
+ return e.map(*this);
+ @}
+@};
+@end example
+
+This function object could then be used like this:
+
+@example
+@{
+ ex M = ... // expression with matrices
+ calc_trace do_trace;
+ ex tr = do_trace(M);
+@}
+@end example
+
+Here is another example for you to meditate over. It removes quadratic
+terms in a variable from an expanded polynomial:
+
+@example
+struct map_rem_quad : public map_function @{
+ ex var;
+ map_rem_quad(const ex & var_) : var(var_) @{@}
+
+ ex operator()(const ex & e)
+ @{
+ if (is_a<add>(e) || is_a<mul>(e))
+ return e.map(*this);
+ else if (is_a<power>(e) &&
+ e.op(0).is_equal(var) && e.op(1).info(info_flags::even))
+ return 0;
+ else
+ return e;
+ @}
+@};
+
+...
+
+@{
+ symbol x("x"), y("y");
+
+ ex e;
+ for (int i=0; i<8; i++)
+ e += pow(x, i) * pow(y, 8-i) * (i+1);
+ cout << e << endl;
+ // -> 4*y^5*x^3+5*y^4*x^4+8*y*x^7+7*y^2*x^6+2*y^7*x+6*y^3*x^5+3*y^6*x^2+y^8
+
+ map_rem_quad rem_quad(x);
+ cout << rem_quad(e) << endl;
+ // -> 4*y^5*x^3+8*y*x^7+2*y^7*x+6*y^3*x^5+y^8
+@}
+@end example
+
+@command{ginsh} offers a slightly different implementation of @code{map()}
+that allows applying algebraic functions to operands. The second argument
+to @code{map()} is an expression containing the wildcard @samp{$0} which
+acts as the placeholder for the operands:
+
+@example
+> map(a*b,sin($0));
+sin(a)*sin(b)
+> map(a+2*b,sin($0));
+sin(a)+sin(2*b)
+> map(@{a,b,c@},$0^2+$0);
+@{a^2+a,b^2+b,c^2+c@}
+@end example
+
+Note that it is only possible to use algebraic functions in the second
+argument. You can not use functions like @samp{diff()}, @samp{op()},
+@samp{subs()} etc. because these are evaluated immediately:
+
+@example
+> map(@{a,b,c@},diff($0,a));
+@{0,0,0@}
+ This is because "diff($0,a)" evaluates to "0", so the command is equivalent
+ to "map(@{a,b,c@},0)".
+@end example
+
+
+@node Visitors and tree traversal, Polynomial arithmetic, Applying a function on subexpressions, Methods and functions
+@c node-name, next, previous, up
+@section Visitors and tree traversal
+@cindex tree traversal
+@cindex @code{visitor} (class)
+@cindex @code{accept()}
+@cindex @code{visit()}
+@cindex @code{traverse()}
+@cindex @code{traverse_preorder()}
+@cindex @code{traverse_postorder()}
+
+Suppose that you need a function that returns a list of all indices appearing
+in an arbitrary expression. The indices can have any dimension, and for
+indices with variance you always want the covariant version returned.
+
+You can't use @code{get_free_indices()} because you also want to include
+dummy indices in the list, and you can't use @code{find()} as it needs
+specific index dimensions (and it would require two passes: one for indices
+with variance, one for plain ones).
+
+The obvious solution to this problem is a tree traversal with a type switch,
+such as the following:
+
+@example
+void gather_indices_helper(const ex & e, lst & l)
+@{
+ if (is_a<varidx>(e)) @{
+ const varidx & vi = ex_to<varidx>(e);
+ l.append(vi.is_covariant() ? vi : vi.toggle_variance());
+ @} else if (is_a<idx>(e)) @{
+ l.append(e);
+ @} else @{
+ size_t n = e.nops();
+ for (size_t i = 0; i < n; ++i)
+ gather_indices_helper(e.op(i), l);
+ @}
+@}
+
+lst gather_indices(const ex & e)
+@{
+ lst l;
+ gather_indices_helper(e, l);
+ l.sort();
+ l.unique();
+ return l;
+@}
+@end example
+
+This works fine but fans of object-oriented programming will feel
+uncomfortable with the type switch. One reason is that there is a possibility
+for subtle bugs regarding derived classes. If we had, for example, written
+
+@example
+ if (is_a<idx>(e)) @{
+ ...
+ @} else if (is_a<varidx>(e)) @{
+ ...
+@end example
+
+in @code{gather_indices_helper}, the code wouldn't have worked because the
+first line "absorbs" all classes derived from @code{idx}, including
+@code{varidx}, so the special case for @code{varidx} would never have been
+executed.
+
+Also, for a large number of classes, a type switch like the above can get
+unwieldy and inefficient (it's a linear search, after all).
+@code{gather_indices_helper} only checks for two classes, but if you had to
+write a function that required a different implementation for nearly
+every GiNaC class, the result would be very hard to maintain and extend.
+
+The cleanest approach to the problem would be to add a new virtual function
+to GiNaC's class hierarchy. In our example, there would be specializations
+for @code{idx} and @code{varidx} while the default implementation in
+@code{basic} performed the tree traversal. Unfortunately, in C++ it's
+impossible to add virtual member functions to existing classes without
+changing their source and recompiling everything. GiNaC comes with source,
+so you could actually do this, but for a small algorithm like the one
+presented this would be impractical.
+
+One solution to this dilemma is the @dfn{Visitor} design pattern,
+which is implemented in GiNaC (actually, Robert Martin's Acyclic Visitor
+variation, described in detail in
+@uref{http://objectmentor.com/publications/acv.pdf}). Instead of adding
+virtual functions to the class hierarchy to implement operations, GiNaC
+provides a single "bouncing" method @code{accept()} that takes an instance
+of a special @code{visitor} class and redirects execution to the one
+@code{visit()} virtual function of the visitor that matches the type of
+object that @code{accept()} was being invoked on.
+
+Visitors in GiNaC must derive from the global @code{visitor} class as well
+as from the class @code{T::visitor} of each class @code{T} they want to
+visit, and implement the member functions @code{void visit(const T &)} for
+each class.
+
+A call of
+
+@example
+void ex::accept(visitor & v) const;
+@end example
+
+will then dispatch to the correct @code{visit()} member function of the
+specified visitor @code{v} for the type of GiNaC object at the root of the
+expression tree (e.g. a @code{symbol}, an @code{idx} or a @code{mul}).
+
+Here is an example of a visitor:
+
+@example
+class my_visitor
+ : public visitor, // this is required
+ public add::visitor, // visit add objects
+ public numeric::visitor, // visit numeric objects
+ public basic::visitor // visit basic objects
+@{
+ void visit(const add & x)
+ @{ cout << "called with an add object" << endl; @}
+
+ void visit(const numeric & x)
+ @{ cout << "called with a numeric object" << endl; @}
+
+ void visit(const basic & x)
+ @{ cout << "called with a basic object" << endl; @}
+@};
+@end example
+
+which can be used as follows:
+
+@example
+...
+ symbol x("x");
+ ex e1 = 42;
+ ex e2 = 4*x-3;
+ ex e3 = 8*x;
+
+ my_visitor v;
+ e1.accept(v);
+ // prints "called with a numeric object"
+ e2.accept(v);
+ // prints "called with an add object"
+ e3.accept(v);
+ // prints "called with a basic object"
+...
+@end example
+
+The @code{visit(const basic &)} method gets called for all objects that are
+not @code{numeric} or @code{add} and acts as an (optional) default.
+
+From a conceptual point of view, the @code{visit()} methods of the visitor
+behave like a newly added virtual function of the visited hierarchy.
+In addition, visitors can store state in member variables, and they can
+be extended by deriving a new visitor from an existing one, thus building
+hierarchies of visitors.
+
+We can now rewrite our index example from above with a visitor:
+
+@example
+class gather_indices_visitor
+ : public visitor, public idx::visitor, public varidx::visitor
+@{
+ lst l;
+
+ void visit(const idx & i)
+ @{
+ l.append(i);
+ @}
+
+ void visit(const varidx & vi)
+ @{
+ l.append(vi.is_covariant() ? vi : vi.toggle_variance());
+ @}
+
+public:
+ const lst & get_result() // utility function
+ @{
+ l.sort();
+ l.unique();
+ return l;
+ @}
+@};
+@end example
+
+What's missing is the tree traversal. We could implement it in
+@code{visit(const basic &)}, but GiNaC has predefined methods for this:
+
+@example
+void ex::traverse_preorder(visitor & v) const;
+void ex::traverse_postorder(visitor & v) const;
+void ex::traverse(visitor & v) const;
+@end example
+
+@code{traverse_preorder()} visits a node @emph{before} visiting its
+subexpressions, while @code{traverse_postorder()} visits a node @emph{after}
+visiting its subexpressions. @code{traverse()} is a synonym for
+@code{traverse_preorder()}.
+
+Here is a new implementation of @code{gather_indices()} that uses the visitor
+and @code{traverse()}:
+
+@example
+lst gather_indices(const ex & e)
+@{
+ gather_indices_visitor v;
+ e.traverse(v);
+ return v.get_result();
+@}
+@end example
+
+Alternatively, you could use pre- or postorder iterators for the tree
+traversal:
+
+@example
+lst gather_indices(const ex & e)
+@{
+ gather_indices_visitor v;
+ for (const_preorder_iterator i = e.preorder_begin();
+ i != e.preorder_end(); ++i) @{
+ i->accept(v);
+ @}
+ return v.get_result();
+@}
+@end example
+
+
+@node Polynomial arithmetic, Rational expressions, Visitors and tree traversal, Methods and functions
+@c node-name, next, previous, up
+@section Polynomial arithmetic
+
+@subsection Testing whether an expression is a polynomial
+@cindex @code{is_polynomial()}
+
+Testing whether an expression is a polynomial in one or more variables
+can be done with the method
+@example
+bool ex::is_polynomial(const ex & vars) const;
+@end example
+In the case of more than
+one variable, the variables are given as a list.
+
+@example
+(x*y*sin(y)).is_polynomial(x) // Returns true.
+(x*y*sin(y)).is_polynomial(lst(x,y)) // Returns false.
+@end example
+
+@subsection Expanding and collecting
+@cindex @code{expand()}
+@cindex @code{collect()}
+@cindex @code{collect_common_factors()}
+
+A polynomial in one or more variables has many equivalent
+representations. Some useful ones serve a specific purpose. Consider
+for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
+21*y*z + 4*z^2} (written down here in output-style). It is equivalent
+to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}. Other
+representations are the recursive ones where one collects for exponents
+in one of the three variable. Since the factors are themselves
+polynomials in the remaining two variables the procedure can be
+repeated. In our example, two possibilities would be @math{(4*y + z)*x
++ 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
+x*z}.
+
+To bring an expression into expanded form, its method
+
+@example
+ex ex::expand(unsigned options = 0);
+@end example
+
+may be called. In our example above, this corresponds to @math{4*x*y +
+x*z + 20*y^2 + 21*y*z + 4*z^2}. Again, since the canonical form in
+GiNaC is not easy to guess you should be prepared to see different
+orderings of terms in such sums!
+
+Another useful representation of multivariate polynomials is as a
+univariate polynomial in one of the variables with the coefficients
+being polynomials in the remaining variables. The method
+@code{collect()} accomplishes this task:
+
+@example
+ex ex::collect(const ex & s, bool distributed = false);
+@end example
+
+The first argument to @code{collect()} can also be a list of objects in which
+case the result is either a recursively collected polynomial, or a polynomial
+in a distributed form with terms like @math{c*x1^e1*...*xn^en}, as specified
+by the @code{distributed} flag.
+
+Note that the original polynomial needs to be in expanded form (for the
+variables concerned) in order for @code{collect()} to be able to find the
+coefficients properly.
+
+The following @command{ginsh} transcript shows an application of @code{collect()}
+together with @code{find()}:
+
+@example
+> a=expand((sin(x)+sin(y))*(1+p+q)*(1+d));
+d*p*sin(x)+p*sin(x)+q*d*sin(x)+q*sin(y)+d*sin(x)+q*d*sin(y)+sin(y)+d*sin(y)
++q*sin(x)+d*sin(y)*p+sin(x)+sin(y)*p
+> collect(a,@{p,q@});
+d*sin(x)+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*p
++(d*sin(x)+sin(y)+d*sin(y)+sin(x))*q+sin(y)+d*sin(y)+sin(x)
+> collect(a,find(a,sin($1)));
+(1+q+d+q*d+d*p+p)*sin(y)+(1+q+d+q*d+d*p+p)*sin(x)
+> collect(a,@{find(a,sin($1)),p,q@});
+(1+(1+d)*p+d+q*(1+d))*sin(x)+(1+(1+d)*p+d+q*(1+d))*sin(y)
+> collect(a,@{find(a,sin($1)),d@});
+(1+q+d*(1+q+p)+p)*sin(y)+(1+q+d*(1+q+p)+p)*sin(x)
+@end example
+
+Polynomials can often be brought into a more compact form by collecting
+common factors from the terms of sums. This is accomplished by the function
+
+@example
+ex collect_common_factors(const ex & e);
+@end example
+
+This function doesn't perform a full factorization but only looks for
+factors which are already explicitly present:
+
+@example
+> collect_common_factors(a*x+a*y);
+(x+y)*a
+> collect_common_factors(a*x^2+2*a*x*y+a*y^2);
+a*(2*x*y+y^2+x^2)
+> collect_common_factors(a*(b*(a+c)*x+b*((a+c)*x+(a+c)*y)*y));
+(c+a)*a*(x*y+y^2+x)*b
+@end example
+
+@subsection Degree and coefficients
+@cindex @code{degree()}
+@cindex @code{ldegree()}