+ ...
+ 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()}
+@cindex @code{coeff()}
+
+The degree and low degree of a polynomial can be obtained using the two
+methods
+
+@example
+int ex::degree(const ex & s);
+int ex::ldegree(const ex & s);
+@end example
+
+which also work reliably on non-expanded input polynomials (they even work
+on rational functions, returning the asymptotic degree). By definition, the
+degree of zero is zero. To extract a coefficient with a certain power from
+an expanded polynomial you use
+
+@example
+ex ex::coeff(const ex & s, int n);
+@end example
+
+You can also obtain the leading and trailing coefficients with the methods
+
+@example
+ex ex::lcoeff(const ex & s);
+ex ex::tcoeff(const ex & s);
+@end example
+
+which are equivalent to @code{coeff(s, degree(s))} and @code{coeff(s, ldegree(s))},
+respectively.
+
+An application is illustrated in the next example, where a multivariate
+polynomial is analyzed:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
+ - pow(x+y,2) + 2*pow(y+2,2) - 8;
+ ex Poly = PolyInp.expand();
+
+ for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) @{
+ cout << "The x^" << i << "-coefficient is "
+ << Poly.coeff(x,i) << endl;
+ @}
+ cout << "As polynomial in y: "
+ << Poly.collect(y) << endl;
+@}
+@end example
+
+When run, it returns an output in the following fashion:
+
+@example
+The x^0-coefficient is y^2+11*y
+The x^1-coefficient is 5*y^2-2*y
+The x^2-coefficient is -1
+The x^3-coefficient is 4*y
+As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
+@end example
+
+As always, the exact output may vary between different versions of GiNaC
+or even from run to run since the internal canonical ordering is not
+within the user's sphere of influence.
+
+@code{degree()}, @code{ldegree()}, @code{coeff()}, @code{lcoeff()},
+@code{tcoeff()} and @code{collect()} can also be used to a certain degree
+with non-polynomial expressions as they not only work with symbols but with
+constants, functions and indexed objects as well:
+
+@example
+@{
+ symbol a("a"), b("b"), c("c"), x("x");
+ idx i(symbol("i"), 3);
+
+ ex e = pow(sin(x) - cos(x), 4);
+ cout << e.degree(cos(x)) << endl;
+ // -> 4
+ cout << e.expand().coeff(sin(x), 3) << endl;
+ // -> -4*cos(x)
+
+ e = indexed(a+b, i) * indexed(b+c, i);
+ e = e.expand(expand_options::expand_indexed);
+ cout << e.collect(indexed(b, i)) << endl;
+ // -> a.i*c.i+(a.i+c.i)*b.i+b.i^2
+@}
+@end example
+
+
+@subsection Polynomial division
+@cindex polynomial division
+@cindex quotient
+@cindex remainder
+@cindex pseudo-remainder
+@cindex @code{quo()}
+@cindex @code{rem()}
+@cindex @code{prem()}
+@cindex @code{divide()}
+
+The two functions
+
+@example
+ex quo(const ex & a, const ex & b, const ex & x);
+ex rem(const ex & a, const ex & b, const ex & x);
+@end example
+
+compute the quotient and remainder of univariate polynomials in the variable
+@samp{x}. The results satisfy @math{a = b*quo(a, b, x) + rem(a, b, x)}.
+
+The additional function
+
+@example
+ex prem(const ex & a, const ex & b, const ex & x);
+@end example
+
+computes the pseudo-remainder of @samp{a} and @samp{b} which satisfies
+@math{c*a = b*q + prem(a, b, x)}, where @math{c = b.lcoeff(x) ^ (a.degree(x) - b.degree(x) + 1)}.
+
+Exact division of multivariate polynomials is performed by the function
+
+@example
+bool divide(const ex & a, const ex & b, ex & q);
+@end example
+
+If @samp{b} divides @samp{a} over the rationals, this function returns @code{true}
+and returns the quotient in the variable @code{q}. Otherwise it returns @code{false}
+in which case the value of @code{q} is undefined.
+
+
+@subsection Unit, content and primitive part
+@cindex @code{unit()}
+@cindex @code{content()}
+@cindex @code{primpart()}
+@cindex @code{unitcontprim()}
+
+The methods
+
+@example
+ex ex::unit(const ex & x);
+ex ex::content(const ex & x);
+ex ex::primpart(const ex & x);
+ex ex::primpart(const ex & x, const ex & c);
+@end example
+
+return the unit part, content part, and primitive polynomial of a multivariate
+polynomial with respect to the variable @samp{x} (the unit part being the sign
+of the leading coefficient, the content part being the GCD of the coefficients,
+and the primitive polynomial being the input polynomial divided by the unit and
+content parts). The second variant of @code{primpart()} expects the previously
+calculated content part of the polynomial in @code{c}, which enables it to
+work faster in the case where the content part has already been computed. The
+product of unit, content, and primitive part is the original polynomial.
+
+Additionally, the method
+
+@example
+void ex::unitcontprim(const ex & x, ex & u, ex & c, ex & p);
+@end example
+
+computes the unit, content, and primitive parts in one go, returning them
+in @code{u}, @code{c}, and @code{p}, respectively.
+
+
+@subsection GCD, LCM and resultant
+@cindex GCD
+@cindex LCM
+@cindex @code{gcd()}
+@cindex @code{lcm()}
+
+The functions for polynomial greatest common divisor and least common
+multiple have the synopsis
+
+@example
+ex gcd(const ex & a, const ex & b);
+ex lcm(const ex & a, const ex & b);
+@end example
+
+The functions @code{gcd()} and @code{lcm()} accept two expressions
+@code{a} and @code{b} as arguments and return a new expression, their
+greatest common divisor or least common multiple, respectively. If the
+polynomials @code{a} and @code{b} are coprime @code{gcd(a,b)} returns 1
+and @code{lcm(a,b)} returns the product of @code{a} and @code{b}. Note that all
+the coefficients must be rationals.
+
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+int main()
+@{
+ symbol x("x"), y("y"), z("z");
+ ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
+ ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
+
+ ex P_gcd = gcd(P_a, P_b);
+ // x + 5*y + 4*z
+ ex P_lcm = lcm(P_a, P_b);
+ // 4*x*y^2 + 13*y*x*z + 20*y^3 + 81*y^2*z + 67*y*z^2 + 3*x*z^2 + 12*z^3
+@}
+@end example
+
+@cindex resultant
+@cindex @code{resultant()}
+
+The resultant of two expressions only makes sense with polynomials.
+It is always computed with respect to a specific symbol within the
+expressions. The function has the interface
+
+@example
+ex resultant(const ex & a, const ex & b, const ex & s);
+@end example
+
+Resultants are symmetric in @code{a} and @code{b}. The following example
+computes the resultant of two expressions with respect to @code{x} and
+@code{y}, respectively:
+
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+int main()
+@{
+ symbol x("x"), y("y");
+
+ ex e1 = x+pow(y,2), e2 = 2*pow(x,3)-1; // x+y^2, 2*x^3-1
+ ex r;
+
+ r = resultant(e1, e2, x);
+ // -> 1+2*y^6
+ r = resultant(e1, e2, y);
+ // -> 1-4*x^3+4*x^6
+@}
+@end example
+
+@subsection Square-free decomposition
+@cindex square-free decomposition
+@cindex factorization
+@cindex @code{sqrfree()}
+
+GiNaC still lacks proper factorization support. Some form of
+factorization is, however, easily implemented by noting that factors
+appearing in a polynomial with power two or more also appear in the
+derivative and hence can easily be found by computing the GCD of the
+original polynomial and its derivatives. Any decent system has an
+interface for this so called square-free factorization. So we provide
+one, too:
+@example
+ex sqrfree(const ex & a, const lst & l = lst());
+@end example
+Here is an example that by the way illustrates how the exact form of the
+result may slightly depend on the order of differentiation, calling for
+some care with subsequent processing of the result:
+@example
+ ...
+ symbol x("x"), y("y");
+ ex BiVarPol = expand(pow(2-2*y,3) * pow(1+x*y,2) * pow(x-2*y,2) * (x+y));
+
+ cout << sqrfree(BiVarPol, lst(x,y)) << endl;
+ // -> 8*(1-y)^3*(y*x^2-2*y+x*(1-2*y^2))^2*(y+x)
+
+ cout << sqrfree(BiVarPol, lst(y,x)) << endl;
+ // -> 8*(1-y)^3*(-y*x^2+2*y+x*(-1+2*y^2))^2*(y+x)
+
+ cout << sqrfree(BiVarPol) << endl;
+ // -> depending on luck, any of the above
+ ...
+@end example
+Note also, how factors with the same exponents are not fully factorized
+with this method.
+
+
+@node Rational expressions, Symbolic differentiation, Polynomial arithmetic, Methods and functions
+@c node-name, next, previous, up
+@section Rational expressions
+
+@subsection The @code{normal} method
+@cindex @code{normal()}
+@cindex simplification
+@cindex temporary replacement
+
+Some basic form of simplification of expressions is called for frequently.
+GiNaC provides the method @code{.normal()}, which converts a rational function
+into an equivalent rational function of the form @samp{numerator/denominator}
+where numerator and denominator are coprime. If the input expression is already
+a fraction, it just finds the GCD of numerator and denominator and cancels it,
+otherwise it performs fraction addition and multiplication.
+
+@code{.normal()} can also be used on expressions which are not rational functions
+as it will replace all non-rational objects (like functions or non-integer
+powers) by temporary symbols to bring the expression to the domain of rational
+functions before performing the normalization, and re-substituting these
+symbols afterwards. This algorithm is also available as a separate method
+@code{.to_rational()}, described below.
+
+This means that both expressions @code{t1} and @code{t2} are indeed
+simplified in this little code snippet:
+
+@example
+@{
+ symbol x("x");
+ ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
+ ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
+ std::cout << "t1 is " << t1.normal() << std::endl;
+ std::cout << "t2 is " << t2.normal() << std::endl;
+@}
+@end example
+
+Of course this works for multivariate polynomials too, so the ratio of
+the sample-polynomials from the section about GCD and LCM above would be
+normalized to @code{P_a/P_b} = @code{(4*y+z)/(y+3*z)}.
+
+
+@subsection Numerator and denominator
+@cindex numerator
+@cindex denominator
+@cindex @code{numer()}
+@cindex @code{denom()}
+@cindex @code{numer_denom()}
+
+The numerator and denominator of an expression can be obtained with
+
+@example
+ex ex::numer();
+ex ex::denom();
+ex ex::numer_denom();
+@end example
+
+These functions will first normalize the expression as described above and
+then return the numerator, denominator, or both as a list, respectively.
+If you need both numerator and denominator, calling @code{numer_denom()} is
+faster than using @code{numer()} and @code{denom()} separately.
+
+
+@subsection Converting to a polynomial or rational expression
+@cindex @code{to_polynomial()}
+@cindex @code{to_rational()}
+
+Some of the methods described so far only work on polynomials or rational
+functions. GiNaC provides a way to extend the domain of these functions to
+general expressions by using the temporary replacement algorithm described
+above. You do this by calling
+
+@example
+ex ex::to_polynomial(exmap & m);
+ex ex::to_polynomial(lst & l);
+@end example
+or
+@example
+ex ex::to_rational(exmap & m);
+ex ex::to_rational(lst & l);
+@end example
+
+on the expression to be converted. The supplied @code{exmap} or @code{lst}
+will be filled with the generated temporary symbols and their replacement
+expressions in a format that can be used directly for the @code{subs()}
+method. It can also already contain a list of replacements from an earlier
+application of @code{.to_polynomial()} or @code{.to_rational()}, so it's
+possible to use it on multiple expressions and get consistent results.
+
+The difference between @code{.to_polynomial()} and @code{.to_rational()}
+is probably best illustrated with an example:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex a = 2*x/sin(x) - y/(3*sin(x));
+ cout << a << endl;
+
+ lst lp;
+ ex p = a.to_polynomial(lp);
+ cout << " = " << p << "\n with " << lp << endl;
+ // = symbol3*symbol2*y+2*symbol2*x
+ // with @{symbol2==sin(x)^(-1),symbol3==-1/3@}
+
+ lst lr;
+ ex r = a.to_rational(lr);
+ cout << " = " << r << "\n with " << lr << endl;
+ // = -1/3*symbol4^(-1)*y+2*symbol4^(-1)*x
+ // with @{symbol4==sin(x)@}
+@}
+@end example
+
+The following more useful example will print @samp{sin(x)-cos(x)}:
+
+@example
+@{
+ symbol x("x");
+ ex a = pow(sin(x), 2) - pow(cos(x), 2);
+ ex b = sin(x) + cos(x);
+ ex q;
+ exmap m;
+ divide(a.to_polynomial(m), b.to_polynomial(m), q);
+ cout << q.subs(m) << endl;
+@}
+@end example
+
+
+@node Symbolic differentiation, Series expansion, Rational expressions, Methods and functions
+@c node-name, next, previous, up
+@section Symbolic differentiation
+@cindex differentiation
+@cindex @code{diff()}
+@cindex chain rule
+@cindex product rule
+
+GiNaC's objects know how to differentiate themselves. Thus, a
+polynomial (class @code{add}) knows that its derivative is the sum of
+the derivatives of all the monomials:
+
+@example
+@{
+ symbol x("x"), y("y"), z("z");
+ ex P = pow(x, 5) + pow(x, 2) + y;
+
+ cout << P.diff(x,2) << endl;
+ // -> 20*x^3 + 2
+ cout << P.diff(y) << endl; // 1
+ // -> 1
+ cout << P.diff(z) << endl; // 0
+ // -> 0
+@}
+@end example
+
+If a second integer parameter @var{n} is given, the @code{diff} method
+returns the @var{n}th derivative.
+
+If @emph{every} object and every function is told what its derivative
+is, all derivatives of composed objects can be calculated using the
+chain rule and the product rule. Consider, for instance the expression
+@code{1/cosh(x)}. Since the derivative of @code{cosh(x)} is
+@code{sinh(x)} and the derivative of @code{pow(x,-1)} is
+@code{-pow(x,-2)}, GiNaC can readily compute the composition. It turns
+out that the composition is the generating function for Euler Numbers,
+i.e. the so called @var{n}th Euler number is the coefficient of
+@code{x^n/n!} in the expansion of @code{1/cosh(x)}. We may use this
+identity to code a function that generates Euler numbers in just three
+lines:
+
+@cindex Euler numbers
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+ex EulerNumber(unsigned n)
+@{
+ symbol x;
+ const ex generator = pow(cosh(x),-1);
+ return generator.diff(x,n).subs(x==0);
+@}
+
+int main()
+@{
+ for (unsigned i=0; i<11; i+=2)
+ std::cout << EulerNumber(i) << std::endl;
+ return 0;
+@}
+@end example
+
+When you run it, it produces the sequence @code{1}, @code{-1}, @code{5},
+@code{-61}, @code{1385}, @code{-50521}. We increment the loop variable
+@code{i} by two since all odd Euler numbers vanish anyways.
+
+
+@node Series expansion, Symmetrization, Symbolic differentiation, Methods and functions
+@c node-name, next, previous, up
+@section Series expansion
+@cindex @code{series()}
+@cindex Taylor expansion
+@cindex Laurent expansion
+@cindex @code{pseries} (class)
+@cindex @code{Order()}
+
+Expressions know how to expand themselves as a Taylor series or (more
+generally) a Laurent series. As in most conventional Computer Algebra
+Systems, no distinction is made between those two. There is a class of
+its own for storing such series (@code{class pseries}) and a built-in
+function (called @code{Order}) for storing the order term of the series.
+As a consequence, if you want to work with series, i.e. multiply two
+series, you need to call the method @code{ex::series} again to convert
+it to a series object with the usual structure (expansion plus order
+term). A sample application from special relativity could read:
+
+@example
+#include <ginac/ginac.h>
+using namespace std;
+using namespace GiNaC;
+
+int main()
+@{
+ symbol v("v"), c("c");
+
+ ex gamma = 1/sqrt(1 - pow(v/c,2));
+ ex mass_nonrel = gamma.series(v==0, 10);
+
+ cout << "the relativistic mass increase with v is " << endl
+ << mass_nonrel << endl;
+
+ cout << "the inverse square of this series is " << endl
+ << pow(mass_nonrel,-2).series(v==0, 10) << endl;
+@}
+@end example
+
+Only calling the series method makes the last output simplify to
+@math{1-v^2/c^2+O(v^10)}, without that call we would just have a long
+series raised to the power @math{-2}.
+
+@cindex Machin's formula
+As another instructive application, let us calculate the numerical
+value of Archimedes' constant
+@tex
+$\pi$
+@end tex
+@ifnottex
+@math{Pi}
+@end ifnottex
+(for which there already exists the built-in constant @code{Pi})
+using John Machin's amazing formula
+@tex
+$\pi=16$~atan~$\!\left(1 \over 5 \right)-4$~atan~$\!\left(1 \over 239 \right)$.
+@end tex
+@ifnottex
+@math{Pi==16*atan(1/5)-4*atan(1/239)}.
+@end ifnottex
+This equation (and similar ones) were used for over 200 years for
+computing digits of pi (see @cite{Pi Unleashed}). We may expand the
+arcus tangent around @code{0} and insert the fractions @code{1/5} and
+@code{1/239}. However, as we have seen, a series in GiNaC carries an
+order term with it and the question arises what the system is supposed
+to do when the fractions are plugged into that order term. The solution
+is to use the function @code{series_to_poly()} to simply strip the order
+term off:
+
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+ex machin_pi(int degr)
+@{
+ symbol x;
+ ex pi_expansion = series_to_poly(atan(x).series(x,degr));
+ ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
+ -4*pi_expansion.subs(x==numeric(1,239));
+ return pi_approx;
+@}
+
+int main()
+@{
+ using std::cout; // just for fun, another way of...
+ using std::endl; // ...dealing with this namespace std.
+ ex pi_frac;
+ for (int i=2; i<12; i+=2) @{
+ pi_frac = machin_pi(i);
+ cout << i << ":\t" << pi_frac << endl
+ << "\t" << pi_frac.evalf() << endl;
+ @}
+ return 0;
+@}
+@end example
+
+Note how we just called @code{.series(x,degr)} instead of
+@code{.series(x==0,degr)}. This is a simple shortcut for @code{ex}'s
+method @code{series()}: if the first argument is a symbol the expression
+is expanded in that symbol around point @code{0}. When you run this
+program, it will type out:
+
+@example
+2: 3804/1195
+ 3.1832635983263598326
+4: 5359397032/1706489875
+ 3.1405970293260603143
+6: 38279241713339684/12184551018734375
+ 3.141621029325034425
+8: 76528487109180192540976/24359780855939418203125
+ 3.141591772182177295
+10: 327853873402258685803048818236/104359128170408663038552734375
+ 3.1415926824043995174
+@end example
+
+
+@node Symmetrization, Built-in functions, Series expansion, Methods and functions
+@c node-name, next, previous, up
+@section Symmetrization
+@cindex @code{symmetrize()}
+@cindex @code{antisymmetrize()}
+@cindex @code{symmetrize_cyclic()}
+
+The three methods
+
+@example
+ex ex::symmetrize(const lst & l);
+ex ex::antisymmetrize(const lst & l);
+ex ex::symmetrize_cyclic(const lst & l);
+@end example
+
+symmetrize an expression by returning the sum over all symmetric,
+antisymmetric or cyclic permutations of the specified list of objects,
+weighted by the number of permutations.
+
+The three additional methods
+
+@example
+ex ex::symmetrize();
+ex ex::antisymmetrize();
+ex ex::symmetrize_cyclic();
+@end example
+
+symmetrize or antisymmetrize an expression over its free indices.
+
+Symmetrization is most useful with indexed expressions but can be used with
+almost any kind of object (anything that is @code{subs()}able):
+
+@example
+@{
+ idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
+ symbol A("A"), B("B"), a("a"), b("b"), c("c");
+
+ cout << indexed(A, i, j).symmetrize() << endl;
+ // -> 1/2*A.j.i+1/2*A.i.j
+ cout << indexed(A, i, j, k).antisymmetrize(lst(i, j)) << endl;
+ // -> -1/2*A.j.i.k+1/2*A.i.j.k
+ cout << lst(a, b, c).symmetrize_cyclic(lst(a, b, c)) << endl;
+ // -> 1/3*@{a,b,c@}+1/3*@{b,c,a@}+1/3*@{c,a,b@}
+@}
+@end example
+
+@page
+
+@node Built-in functions, Multiple polylogarithms, Symmetrization, Methods and functions
+@c node-name, next, previous, up
+@section Predefined mathematical functions
+@c
+@subsection Overview
+
+GiNaC contains the following predefined mathematical functions:
+
+@cartouche
+@multitable @columnfractions .30 .70
+@item @strong{Name} @tab @strong{Function}
+@item @code{abs(x)}
+@tab absolute value
+@cindex @code{abs()}
+@item @code{step(x)}
+@tab step function
+@cindex @code{step()}
+@item @code{csgn(x)}
+@tab complex sign
+@cindex @code{conjugate()}
+@item @code{conjugate(x)}
+@tab complex conjugation
+@cindex @code{real_part()}
+@item @code{real_part(x)}
+@tab real part
+@cindex @code{imag_part()}
+@item @code{imag_part(x)}
+@tab imaginary part
+@item @code{sqrt(x)}
+@tab square root (not a GiNaC function, rather an alias for @code{pow(x, numeric(1, 2))})
+@cindex @code{sqrt()}
+@item @code{sin(x)}
+@tab sine
+@cindex @code{sin()}
+@item @code{cos(x)}
+@tab cosine
+@cindex @code{cos()}
+@item @code{tan(x)}
+@tab tangent
+@cindex @code{tan()}
+@item @code{asin(x)}
+@tab inverse sine
+@cindex @code{asin()}
+@item @code{acos(x)}
+@tab inverse cosine
+@cindex @code{acos()}
+@item @code{atan(x)}
+@tab inverse tangent
+@cindex @code{atan()}
+@item @code{atan2(y, x)}
+@tab inverse tangent with two arguments
+@item @code{sinh(x)}
+@tab hyperbolic sine
+@cindex @code{sinh()}
+@item @code{cosh(x)}
+@tab hyperbolic cosine
+@cindex @code{cosh()}
+@item @code{tanh(x)}
+@tab hyperbolic tangent
+@cindex @code{tanh()}
+@item @code{asinh(x)}
+@tab inverse hyperbolic sine
+@cindex @code{asinh()}
+@item @code{acosh(x)}
+@tab inverse hyperbolic cosine
+@cindex @code{acosh()}
+@item @code{atanh(x)}
+@tab inverse hyperbolic tangent
+@cindex @code{atanh()}
+@item @code{exp(x)}
+@tab exponential function
+@cindex @code{exp()}
+@item @code{log(x)}
+@tab natural logarithm
+@cindex @code{log()}
+@item @code{Li2(x)}
+@tab dilogarithm
+@cindex @code{Li2()}
+@item @code{Li(m, x)}
+@tab classical polylogarithm as well as multiple polylogarithm
+@cindex @code{Li()}
+@item @code{G(a, y)}
+@tab multiple polylogarithm
+@cindex @code{G()}
+@item @code{G(a, s, y)}
+@tab multiple polylogarithm with explicit signs for the imaginary parts
+@cindex @code{G()}
+@item @code{S(n, p, x)}
+@tab Nielsen's generalized polylogarithm
+@cindex @code{S()}
+@item @code{H(m, x)}
+@tab harmonic polylogarithm
+@cindex @code{H()}
+@item @code{zeta(m)}
+@tab Riemann's zeta function as well as multiple zeta value
+@cindex @code{zeta()}
+@item @code{zeta(m, s)}
+@tab alternating Euler sum
+@cindex @code{zeta()}
+@item @code{zetaderiv(n, x)}
+@tab derivatives of Riemann's zeta function
+@item @code{tgamma(x)}
+@tab gamma function
+@cindex @code{tgamma()}
+@cindex gamma function
+@item @code{lgamma(x)}
+@tab logarithm of gamma function
+@cindex @code{lgamma()}
+@item @code{beta(x, y)}
+@tab beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
+@cindex @code{beta()}
+@item @code{psi(x)}
+@tab psi (digamma) function
+@cindex @code{psi()}
+@item @code{psi(n, x)}
+@tab derivatives of psi function (polygamma functions)
+@item @code{factorial(n)}
+@tab factorial function @math{n!}
+@cindex @code{factorial()}
+@item @code{binomial(n, k)}
+@tab binomial coefficients
+@cindex @code{binomial()}
+@item @code{Order(x)}
+@tab order term function in truncated power series
+@cindex @code{Order()}
+@end multitable
+@end cartouche
+
+@cindex branch cut
+For functions that have a branch cut in the complex plane GiNaC follows
+the conventions for C++ as defined in the ANSI standard as far as
+possible. In particular: the natural logarithm (@code{log}) and the
+square root (@code{sqrt}) both have their branch cuts running along the
+negative real axis where the points on the axis itself belong to the
+upper part (i.e. continuous with quadrant II). The inverse
+trigonometric and hyperbolic functions are not defined for complex
+arguments by the C++ standard, however. In GiNaC we follow the
+conventions used by CLN, which in turn follow the carefully designed
+definitions in the Common Lisp standard. It should be noted that this
+convention is identical to the one used by the C99 standard and by most
+serious CAS. It is to be expected that future revisions of the C++
+standard incorporate these functions in the complex domain in a manner
+compatible with C99.
+
+@node Multiple polylogarithms, Complex expressions, Built-in functions, Methods and functions
+@c node-name, next, previous, up
+@subsection Multiple polylogarithms
+
+@cindex polylogarithm
+@cindex Nielsen's generalized polylogarithm
+@cindex harmonic polylogarithm
+@cindex multiple zeta value
+@cindex alternating Euler sum
+@cindex multiple polylogarithm
+
+The multiple polylogarithm is the most generic member of a family of functions,
+to which others like the harmonic polylogarithm, Nielsen's generalized
+polylogarithm and the multiple zeta value belong.
+Everyone of these functions can also be written as a multiple polylogarithm with specific
+parameters. This whole family of functions is therefore often referred to simply as
+multiple polylogarithms, containing @code{Li}, @code{G}, @code{H}, @code{S} and @code{zeta}.
+The multiple polylogarithm itself comes in two variants: @code{Li} and @code{G}. While
+@code{Li} and @code{G} in principle represent the same function, the different
+notations are more natural to the series representation or the integral
+representation, respectively.
+
+To facilitate the discussion of these functions we distinguish between indices and
+arguments as parameters. In the table above indices are printed as @code{m}, @code{s},
+@code{n} or @code{p}, whereas arguments are printed as @code{x}, @code{a} and @code{y}.
+
+To define a @code{Li}, @code{H} or @code{zeta} with a depth greater than one, you have to
+pass a GiNaC @code{lst} for the indices @code{m} and @code{s}, and in the case of @code{Li}
+for the argument @code{x} as well. The parameter @code{a} of @code{G} must always be a @code{lst} containing
+the arguments in expanded form. If @code{G} is used with a third parameter @code{s}, @code{s} must
+have the same length as @code{a}. It contains then the signs of the imaginary parts of the arguments. If
+@code{s} is not given, the signs default to +1.
+Note that @code{Li} and @code{zeta} are polymorphic in this respect. They can stand in for
+the classical polylogarithm and Riemann's zeta function (if depth is one), as well as for
+the multiple polylogarithm and the multiple zeta value, respectively. Note also, that
+GiNaC doesn't check whether the @code{lst}s for two parameters do have the same length.
+It is up to the user to ensure this, otherwise evaluating will result in undefined behavior.
+
+The functions print in LaTeX format as
+@tex
+${\rm Li\;\!}_{m_1,m_2,\ldots,m_k}(x_1,x_2,\ldots,x_k)$,
+@end tex
+@tex
+${\rm S}_{n,p}(x)$,
+@end tex
+@tex
+${\rm H\;\!}_{m_1,m_2,\ldots,m_k}(x)$ and
+@end tex
+@tex
+$\zeta(m_1,m_2,\ldots,m_k)$.
+@end tex
+@ifnottex
+@command{\mbox@{Li@}_@{m_1,m_2,...,m_k@}(x_1,x_2,...,x_k)},
+@command{\mbox@{S@}_@{n,p@}(x)},
+@command{\mbox@{H@}_@{m_1,m_2,...,m_k@}(x)} and
+@command{\zeta(m_1,m_2,...,m_k)} (with the dots replaced by actual parameters).
+@end ifnottex
+If @code{zeta} is an alternating zeta sum, i.e. @code{zeta(m,s)}, the indices with negative sign
+are printed with a line above, e.g.
+@tex
+$\zeta(5,\overline{2})$.
+@end tex
+@ifnottex
+@command{\zeta(5,\overline@{2@})}.
+@end ifnottex
+The order of indices and arguments in the GiNaC @code{lst}s and in the output is the same.
+
+Definitions and analytical as well as numerical properties of multiple polylogarithms
+are too numerous to be covered here. Instead, the user is referred to the publications listed at the
+end of this section. The implementation in GiNaC adheres to the definitions and conventions therein,
+except for a few differences which will be explicitly stated in the following.
+
+One difference is about the order of the indices and arguments. For GiNaC we adopt the convention
+that the indices and arguments are understood to be in the same order as in which they appear in
+the series representation. This means
+@tex
+${\rm Li\;\!}_{m_1,m_2,m_3}(x,1,1) = {\rm H\;\!}_{m_1,m_2,m_3}(x)$ and
+@end tex
+@tex
+${\rm Li\;\!}_{2,1}(1,1) = \zeta(2,1) = \zeta(3)$, but
+@end tex
+@tex
+$\zeta(1,2)$ evaluates to infinity.
+@end tex
+@ifnottex
+@code{Li_@{m_1,m_2,m_3@}(x,1,1) = H_@{m_1,m_2,m_3@}(x)} and
+@code{Li_@{2,1@}(1,1) = zeta(2,1) = zeta(3)}, but
+@code{zeta(1,2)} evaluates to infinity.
+@end ifnottex
+So in comparison to the older ones of the referenced publications the order of
+indices and arguments for @code{Li} is reversed.
+
+The functions only evaluate if the indices are integers greater than zero, except for the indices
+@code{s} in @code{zeta} and @code{G} as well as @code{m} in @code{H}. Since @code{s}
+will be interpreted as the sequence of signs for the corresponding indices
+@code{m} or the sign of the imaginary part for the
+corresponding arguments @code{a}, it must contain 1 or -1, e.g.
+@code{zeta(lst(3,4), lst(-1,1))} means
+@tex
+$\zeta(\overline{3},4)$
+@end tex
+@ifnottex
+@command{zeta(\overline@{3@},4)}
+@end ifnottex
+and
+@code{G(lst(a,b), lst(-1,1), c)} means
+@tex
+$G(a-0\epsilon,b+0\epsilon;c)$.
+@end tex
+@ifnottex
+@command{G(a-0\epsilon,b+0\epsilon;c)}.
+@end ifnottex
+The definition of @code{H} allows indices to be 0, 1 or -1 (in expanded notation) or equally to
+be any integer (in compact notation). With GiNaC expanded and compact notation can be mixed,
+e.g. @code{lst(0,0,-1,0,1,0,0)}, @code{lst(0,0,-1,2,0,0)} and @code{lst(-3,2,0,0)} are equivalent as
+indices. The anonymous evaluator @code{eval()} tries to reduce the functions, if possible, to
+the least-generic multiple polylogarithm. If all arguments are unit, it returns @code{zeta}.
+Arguments equal to zero get considered, too. Riemann's zeta function @code{zeta} (with depth one)
+evaluates also for negative integers and positive even integers. For example:
+
+@example
+> Li(@{3,1@},@{x,1@});
+S(2,2,x)
+> H(@{-3,2@},1);
+-zeta(@{3,2@},@{-1,-1@})
+> S(3,1,1);
+1/90*Pi^4
+@end example
+
+It is easy to tell for a given function into which other function it can be rewritten, may
+it be a less-generic or a more-generic one, except for harmonic polylogarithms @code{H}
+with negative indices or trailing zeros (the example above gives a hint). Signs can
+quickly be messed up, for example. Therefore GiNaC offers a C++ function
+@code{convert_H_to_Li()} to deal with the upgrade of a @code{H} to a multiple polylogarithm
+@code{Li} (@code{eval()} already cares for the possible downgrade):
+
+@example
+> convert_H_to_Li(@{0,-2,-1,3@},x);
+Li(@{3,1,3@},@{-x,1,-1@})
+> convert_H_to_Li(@{2,-1,0@},x);
+-Li(@{2,1@},@{x,-1@})*log(x)+2*Li(@{3,1@},@{x,-1@})+Li(@{2,2@},@{x,-1@})
+@end example
+
+Every function can be numerically evaluated for
+arbitrary real or complex arguments. The precision is arbitrary and can be set through the
+global variable @code{Digits}:
+
+@example
+> Digits=100;
+100
+> evalf(zeta(@{3,1,3,1@}));
+0.005229569563530960100930652283899231589890420784634635522547448972148869544...
+@end example
+
+Note that the convention for arguments on the branch cut in GiNaC as stated above is
+different from the one Remiddi and Vermaseren have chosen for the harmonic polylogarithm.
+
+If a function evaluates to infinity, no exceptions are raised, but the function is returned
+unevaluated, e.g.
+@tex
+$\zeta(1)$.
+@end tex
+@ifnottex
+@command{zeta(1)}.
+@end ifnottex
+In long expressions this helps a lot with debugging, because you can easily spot
+the divergencies. But on the other hand, you have to make sure for yourself, that no illegal
+cancellations of divergencies happen.
+
+Useful publications:
+
+@cite{Nested Sums, Expansion of Transcendental Functions and Multi-Scale Multi-Loop Integrals},
+S.Moch, P.Uwer, S.Weinzierl, hep-ph/0110083
+
+@cite{Harmonic Polylogarithms},
+E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
+
+@cite{Special Values of Multiple Polylogarithms},
+J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001), pp. 907-941
+
+@cite{Numerical Evaluation of Multiple Polylogarithms},
+J.Vollinga, S.Weinzierl, hep-ph/0410259
+
+@node Complex expressions, Solving linear systems of equations, Multiple polylogarithms, Methods and functions
+@c node-name, next, previous, up
+@section Complex expressions
+@c
+@cindex @code{conjugate()}
+
+For dealing with complex expressions there are the methods
+
+@example
+ex ex::conjugate();
+ex ex::real_part();
+ex ex::imag_part();
+@end example
+
+that return respectively the complex conjugate, the real part and the
+imaginary part of an expression. Complex conjugation works as expected
+for all built-in functions and objects. Taking real and imaginary
+parts has not yet been implemented for all built-in functions. In cases where
+it is not known how to conjugate or take a real/imaginary part one
+of the functions @code{conjugate}, @code{real_part} or @code{imag_part}
+is returned. For instance, in case of a complex symbol @code{x}
+(symbols are complex by default), one could not simplify
+@code{conjugate(x)}. In the case of strings of gamma matrices,
+the @code{conjugate} method takes the Dirac conjugate.
+
+For example,
+@example
+@{
+ varidx a(symbol("a"), 4), b(symbol("b"), 4);
+ symbol x("x");
+ realsymbol y("y");
+
+ cout << (3*I*x*y + sin(2*Pi*I*y)).conjugate() << endl;
+ // -> -3*I*conjugate(x)*y+sin(-2*I*Pi*y)
+ cout << (dirac_gamma(a)*dirac_gamma(b)*dirac_gamma5()).conjugate() << endl;
+ // -> -gamma5*gamma~b*gamma~a
+@}
+@end example
+
+If you declare your own GiNaC functions, then they will conjugate themselves
+by conjugating their arguments. This is the default strategy. If you want to
+change this behavior, you have to supply a specialized conjugation method
+for your function (see @ref{Symbolic functions} and the GiNaC source-code
+for @code{abs} as an example). Also, specialized methods can be provided
+to take real and imaginary parts of user-defined functions.
+
+@node Solving linear systems of equations, Input/output, Complex expressions, Methods and functions
+@c node-name, next, previous, up
+@section Solving linear systems of equations
+@cindex @code{lsolve()}
+
+The function @code{lsolve()} provides a convenient wrapper around some
+matrix operations that comes in handy when a system of linear equations
+needs to be solved:
+
+@example
+ex lsolve(const ex & eqns, const ex & symbols,
+ unsigned options = solve_algo::automatic);
+@end example
+
+Here, @code{eqns} is a @code{lst} of equalities (i.e. class
+@code{relational}) while @code{symbols} is a @code{lst} of
+indeterminates. (@xref{The class hierarchy}, for an exposition of class
+@code{lst}).
+
+It returns the @code{lst} of solutions as an expression. As an example,
+let us solve the two equations @code{a*x+b*y==3} and @code{x-y==b}:
+
+@example
+@{
+ symbol a("a"), b("b"), x("x"), y("y");
+ lst eqns, vars;
+ eqns = a*x+b*y==3, x-y==b;
+ vars = x, y;
+ cout << lsolve(eqns, vars) << endl;
+ // -> @{x==(3+b^2)/(b+a),y==(3-b*a)/(b+a)@}
+@end example
+
+When the linear equations @code{eqns} are underdetermined, the solution
+will contain one or more tautological entries like @code{x==x},
+depending on the rank of the system. When they are overdetermined, the
+solution will be an empty @code{lst}. Note the third optional parameter
+to @code{lsolve()}: it accepts the same parameters as
+@code{matrix::solve()}. This is because @code{lsolve} is just a wrapper
+around that method.
+
+
+@node Input/output, Extending GiNaC, Solving linear systems of equations, Methods and functions
+@c node-name, next, previous, up
+@section Input and output of expressions
+@cindex I/O
+
+@subsection Expression output
+@cindex printing
+@cindex output of expressions
+
+Expressions can simply be written to any stream:
+
+@example
+@{
+ symbol x("x");
+ ex e = 4.5*I+pow(x,2)*3/2;
+ cout << e << endl; // prints '4.5*I+3/2*x^2'
+ // ...
+@end example
+
+The default output format is identical to the @command{ginsh} input syntax and
+to that used by most computer algebra systems, but not directly pastable
+into a GiNaC C++ program (note that in the above example, @code{pow(x,2)}
+is printed as @samp{x^2}).
+
+It is possible to print expressions in a number of different formats with
+a set of stream manipulators;
+
+@example
+std::ostream & dflt(std::ostream & os);
+std::ostream & latex(std::ostream & os);
+std::ostream & tree(std::ostream & os);
+std::ostream & csrc(std::ostream & os);
+std::ostream & csrc_float(std::ostream & os);
+std::ostream & csrc_double(std::ostream & os);
+std::ostream & csrc_cl_N(std::ostream & os);
+std::ostream & index_dimensions(std::ostream & os);
+std::ostream & no_index_dimensions(std::ostream & os);
+@end example
+
+The @code{tree}, @code{latex} and @code{csrc} formats are also available in
+@command{ginsh} via the @code{print()}, @code{print_latex()} and
+@code{print_csrc()} functions, respectively.
+
+@cindex @code{dflt}
+All manipulators affect the stream state permanently. To reset the output
+format to the default, use the @code{dflt} manipulator:
+
+@example
+ // ...
+ cout << latex; // all output to cout will be in LaTeX format from
+ // now on
+ cout << e << endl; // prints '4.5 i+\frac@{3@}@{2@} x^@{2@}'
+ cout << sin(x/2) << endl; // prints '\sin(\frac@{1@}@{2@} x)'
+ cout << dflt; // revert to default output format
+ cout << e << endl; // prints '4.5*I+3/2*x^2'
+ // ...
+@end example
+
+If you don't want to affect the format of the stream you're working with,
+you can output to a temporary @code{ostringstream} like this:
+
+@example
+ // ...
+ ostringstream s;
+ s << latex << e; // format of cout remains unchanged
+ cout << s.str() << endl; // prints '4.5 i+\frac@{3@}@{2@} x^@{2@}'
+ // ...
+@end example
+
+@anchor{csrc printing}
+@cindex @code{csrc}
+@cindex @code{csrc_float}
+@cindex @code{csrc_double}
+@cindex @code{csrc_cl_N}
+The @code{csrc} (an alias for @code{csrc_double}), @code{csrc_float},
+@code{csrc_double} and @code{csrc_cl_N} manipulators set the output to a
+format that can be directly used in a C or C++ program. The three possible
+formats select the data types used for numbers (@code{csrc_cl_N} uses the
+classes provided by the CLN library):
+
+@example
+ // ...
+ cout << "f = " << csrc_float << e << ";\n";
+ cout << "d = " << csrc_double << e << ";\n";
+ cout << "n = " << csrc_cl_N << e << ";\n";
+ // ...
+@end example
+
+The above example will produce (note the @code{x^2} being converted to
+@code{x*x}):
+
+@example
+f = (3.0/2.0)*(x*x)+std::complex<float>(0.0,4.5000000e+00);
+d = (3.0/2.0)*(x*x)+std::complex<double>(0.0,4.5000000000000000e+00);
+n = cln::cl_RA("3/2")*(x*x)+cln::complex(cln::cl_I("0"),cln::cl_F("4.5_17"));
+@end example
+
+@cindex @code{tree}
+The @code{tree} manipulator allows dumping the internal structure of an
+expression for debugging purposes:
+
+@example
+ // ...
+ cout << tree << e;
+@}
+@end example
+
+produces
+
+@example
+add, hash=0x0, flags=0x3, nops=2
+ power, hash=0x0, flags=0x3, nops=2
+ x (symbol), serial=0, hash=0xc8d5bcdd, flags=0xf
+ 2 (numeric), hash=0x6526b0fa, flags=0xf
+ 3/2 (numeric), hash=0xf9828fbd, flags=0xf
+ -----
+ overall_coeff
+ 4.5L0i (numeric), hash=0xa40a97e0, flags=0xf
+ =====
+@end example
+
+@cindex @code{latex}
+The @code{latex} output format is for LaTeX parsing in mathematical mode.
+It is rather similar to the default format but provides some braces needed
+by LaTeX for delimiting boxes and also converts some common objects to
+conventional LaTeX names. It is possible to give symbols a special name for
+LaTeX output by supplying it as a second argument to the @code{symbol}
+constructor.
+
+For example, the code snippet
+
+@example
+@{
+ symbol x("x", "\\circ");
+ ex e = lgamma(x).series(x==0,3);
+ cout << latex << e << endl;
+@}
+@end example
+
+will print
+
+@example
+ @{(-\ln(\circ))@}+@{(-\gamma_E)@} \circ+@{(\frac@{1@}@{12@} \pi^@{2@})@} \circ^@{2@}
+ +\mathcal@{O@}(\circ^@{3@})
+@end example
+
+@cindex @code{index_dimensions}
+@cindex @code{no_index_dimensions}
+Index dimensions are normally hidden in the output. To make them visible, use
+the @code{index_dimensions} manipulator. The dimensions will be written in
+square brackets behind each index value in the default and LaTeX output
+formats:
+
+@example
+@{
+ symbol x("x"), y("y");
+ varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
+ ex e = indexed(x, mu) * indexed(y, nu);
+
+ cout << e << endl;
+ // prints 'x~mu*y~nu'
+ cout << index_dimensions << e << endl;
+ // prints 'x~mu[4]*y~nu[4]'
+ cout << no_index_dimensions << e << endl;
+ // prints 'x~mu*y~nu'
+@}
+@end example
+
+
+@cindex Tree traversal
+If you need any fancy special output format, e.g. for interfacing GiNaC
+with other algebra systems or for producing code for different
+programming languages, you can always traverse the expression tree yourself:
+
+@example
+static void my_print(const ex & e)
+@{
+ if (is_a<function>(e))
+ cout << ex_to<function>(e).get_name();
+ else
+ cout << ex_to<basic>(e).class_name();
+ cout << "(";
+ size_t n = e.nops();
+ if (n)
+ for (size_t i=0; i<n; i++) @{
+ my_print(e.op(i));
+ if (i != n-1)
+ cout << ",";
+ @}
+ else
+ cout << e;
+ cout << ")";
+@}
+
+int main()
+@{
+ my_print(pow(3, x) - 2 * sin(y / Pi)); cout << endl;
+ return 0;
+@}
+@end example
+
+This will produce
+
+@example
+add(power(numeric(3),symbol(x)),mul(sin(mul(power(constant(Pi),numeric(-1)),
+symbol(y))),numeric(-2)))
+@end example
+
+If you need an output format that makes it possible to accurately
+reconstruct an expression by feeding the output to a suitable parser or
+object factory, you should consider storing the expression in an
+@code{archive} object and reading the object properties from there.
+See the section on archiving for more information.
+
+
+@subsection Expression input
+@cindex input of expressions
+
+GiNaC provides no way to directly read an expression from a stream because
+you will usually want the user to be able to enter something like @samp{2*x+sin(y)}
+and have the @samp{x} and @samp{y} correspond to the symbols @code{x} and
+@code{y} you defined in your program and there is no way to specify the
+desired symbols to the @code{>>} stream input operator.
+
+Instead, GiNaC lets you construct an expression from a string, specifying the
+list of symbols to be used:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex e("2*x+sin(y)", lst(x, y));
+@}
+@end example
+
+The input syntax is the same as that used by @command{ginsh} and the stream
+output operator @code{<<}. The symbols in the string are matched by name to
+the symbols in the list and if GiNaC encounters a symbol not specified in
+the list it will throw an exception.
+
+With this constructor, it's also easy to implement interactive GiNaC programs:
+
+@example
+#include <iostream>
+#include <string>
+#include <stdexcept>
+#include <ginac/ginac.h>
+using namespace std;
+using namespace GiNaC;
+
+int main()
+@{
+ symbol x("x");
+ string s;
+
+ cout << "Enter an expression containing 'x': ";
+ getline(cin, s);
+
+ try @{
+ ex e(s, lst(x));
+ cout << "The derivative of " << e << " with respect to x is ";
+ cout << e.diff(x) << ".\n";
+ @} catch (exception &p) @{
+ cerr << p.what() << endl;
+ @}
+@}
+@end example
+
+@subsection Compiling expressions to C function pointers
+@cindex compiling expressions
+
+Numerical evaluation of algebraic expressions is seamlessly integrated into
+GiNaC by help of the CLN library. While CLN allows for very fast arbitrary
+precision numerics, which is more than sufficient for most users, sometimes only
+the speed of built-in floating point numbers is fast enough, e.g. for Monte
+Carlo integration. The only viable option then is the following: print the
+expression in C syntax format, manually add necessary C code, compile that
+program and run is as a separate application. This is not only cumbersome and
+involves a lot of manual intervention, but it also separates the algebraic and
+the numerical evaluation into different execution stages.
+
+GiNaC offers a couple of functions that help to avoid these inconveniences and
+problems. The functions automatically perform the printing of a GiNaC expression
+and the subsequent compiling of its associated C code. The created object code
+is then dynamically linked to the currently running program. A function pointer
+to the C function that performs the numerical evaluation is returned and can be
+used instantly. This all happens automatically, no user intervention is needed.
+
+The following example demonstrates the use of @code{compile_ex}:
+
+@example
+ // ...
+ symbol x("x");
+ ex myexpr = sin(x) / x;
+
+ FUNCP_1P fp;
+ compile_ex(myexpr, x, fp);
+
+ cout << fp(3.2) << endl;
+ // ...
+@end example
+
+The function @code{compile_ex} is called with the expression to be compiled and
+its only free variable @code{x}. Upon successful completion the third parameter
+contains a valid function pointer to the corresponding C code module. If called
+like in the last line only built-in double precision numerics is involved.
+
+@cindex FUNCP_1P
+@cindex FUNCP_2P
+@cindex FUNCP_CUBA
+The function pointer has to be defined in advance. GiNaC offers three function
+pointer types at the moment:
+
+@example
+ typedef double (*FUNCP_1P) (double);
+ typedef double (*FUNCP_2P) (double, double);
+ typedef void (*FUNCP_CUBA) (const int*, const double[], const int*, double[]);
+@end example
+
+@cindex CUBA library
+@cindex Monte Carlo integration
+@code{FUNCP_2P} allows for two variables in the expression. @code{FUNCP_CUBA} is
+the correct type to be used with the CUBA library
+(@uref{http://www.feynarts/cuba}) for numerical integrations. The details for the
+parameters of @code{FUNCP_CUBA} are explained in the CUBA manual.
+
+@cindex compile_ex
+For every function pointer type there is a matching @code{compile_ex} available:
+
+@example
+ void compile_ex(const ex& expr, const symbol& sym, FUNCP_1P& fp,
+ const std::string filename = "");
+ void compile_ex(const ex& expr, const symbol& sym1, const symbol& sym2,
+ FUNCP_2P& fp, const std::string filename = "");
+ void compile_ex(const lst& exprs, const lst& syms, FUNCP_CUBA& fp,
+ const std::string filename = "");
+@end example
+
+When the last parameter @code{filename} is not supplied, @code{compile_ex} will
+choose a unique random name for the intermediate source and object files it
+produces. On program termination these files will be deleted. If one wishes to
+keep the C code and the object files, one can supply the @code{filename}
+parameter. The intermediate files will use that filename and will not be
+deleted.
+
+@cindex link_ex
+@code{link_ex} is a function that allows to dynamically link an existing object
+file and to make it available via a function pointer. This is useful if you
+have already used @code{compile_ex} on an expression and want to avoid the
+compilation step to be performed over and over again when you restart your
+program. The precondition for this is of course, that you have chosen a
+filename when you did call @code{compile_ex}. For every above mentioned
+function pointer type there exists a corresponding @code{link_ex} function:
+
+@example
+ void link_ex(const std::string filename, FUNCP_1P& fp);
+ void link_ex(const std::string filename, FUNCP_2P& fp);
+ void link_ex(const std::string filename, FUNCP_CUBA& fp);
+@end example
+
+The complete filename (including the suffix @code{.so}) of the object file has
+to be supplied.
+
+The function
+
+@cindex unlink_ex
+@example
+ void unlink_ex(const std::string filename);
+@end example
+
+is supplied for the rare cases when one wishes to close the dynamically linked
+object files directly and have the intermediate files (only if filename has not
+been given) deleted. Normally one doesn't need this function, because all the
+clean-up will be done automatically upon (regular) program termination.
+
+All the described functions will throw an exception in case they cannot perform
+correctly, like for example when writing the file or starting the compiler
+fails. Since internally the same printing methods as described in section
+@ref{csrc printing} are used, only functions and objects that are available in
+standard C will compile successfully (that excludes polylogarithms for example
+at the moment). Another precondition for success is, of course, that it must be
+possible to evaluate the expression numerically. No free variables despite the
+ones supplied to @code{compile_ex} should appear in the expression.
+
+@cindex ginac-excompiler
+@code{compile_ex} uses the shell script @code{ginac-excompiler} to start the C
+compiler and produce the object files. This shell script comes with GiNaC and
+will be installed together with GiNaC in the configured @code{$PREFIX/bin}
+directory.
+
+@subsection Archiving
+@cindex @code{archive} (class)
+@cindex archiving
+
+GiNaC allows creating @dfn{archives} of expressions which can be stored
+to or retrieved from files. To create an archive, you declare an object
+of class @code{archive} and archive expressions in it, giving each
+expression a unique name:
+
+@example
+#include <fstream>
+using namespace std;
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+int main()
+@{
+ symbol x("x"), y("y"), z("z");
+
+ ex foo = sin(x + 2*y) + 3*z + 41;
+ ex bar = foo + 1;
+
+ archive a;
+ a.archive_ex(foo, "foo");
+ a.archive_ex(bar, "the second one");
+ // ...
+@end example
+
+The archive can then be written to a file:
+
+@example
+ // ...
+ ofstream out("foobar.gar");
+ out << a;
+ out.close();
+ // ...
+@end example
+
+The file @file{foobar.gar} contains all information that is needed to
+reconstruct the expressions @code{foo} and @code{bar}.
+
+@cindex @command{viewgar}
+The tool @command{viewgar} that comes with GiNaC can be used to view
+the contents of GiNaC archive files:
+
+@example
+$ viewgar foobar.gar
+foo = 41+sin(x+2*y)+3*z
+the second one = 42+sin(x+2*y)+3*z
+@end example
+
+The point of writing archive files is of course that they can later be
+read in again:
+
+@example
+ // ...
+ archive a2;
+ ifstream in("foobar.gar");
+ in >> a2;
+ // ...