]> www.ginac.de Git - ginac.git/blobdiff - doc/tutorial/ginac.texi
* Added fsolve() numerical univariate real-valued function solver.
[ginac.git] / doc / tutorial / ginac.texi
index 3f98183162437c3a8d7c4bcf03c9cc977c52e2eb..0d337cdb74275d10cb7961edb3d2c0a87f9e830a 100644 (file)
@@ -23,7 +23,7 @@
 This is a tutorial that documents GiNaC @value{VERSION}, an open
 framework for symbolic computation within the C++ programming language.
 
-Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
+Copyright (C) 1999-2005 Johannes Gutenberg University Mainz, Germany
 
 Permission is granted to make and distribute verbatim copies of
 this manual provided the copyright notice and this permission notice
@@ -52,7 +52,7 @@ notice identical to this one.
 
 @page
 @vskip 0pt plus 1filll
-Copyright @copyright{} 1999-2004 Johannes Gutenberg University Mainz, Germany
+Copyright @copyright{} 1999-2005 Johannes Gutenberg University Mainz, Germany
 @sp 2
 Permission is granted to make and distribute verbatim copies of
 this manual provided the copyright notice and this permission notice
@@ -135,7 +135,7 @@ the near future.
 
 @section License
 The GiNaC framework for symbolic computation within the C++ programming
-language is Copyright @copyright{} 1999-2004 Johannes Gutenberg
+language is Copyright @copyright{} 1999-2005 Johannes Gutenberg
 University Mainz, Germany.
 
 This program is free software; you can redistribute it and/or
@@ -150,8 +150,8 @@ General Public License for more details.
 
 You should have received a copy of the GNU General Public License
 along with this program; see the file COPYING.  If not, write to the
-Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
-MA 02111-1307, USA.
+Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston,
+MA 02110-1301, USA.
 
 
 @node A Tour of GiNaC, How to use it from within C++, Introduction, Top
@@ -417,6 +417,27 @@ x^(-1)-0.5772156649015328606+(0.9890559953279725555)*x
 Here we have made use of the @command{ginsh}-command @code{%} to pop the
 previously evaluated element from @command{ginsh}'s internal stack.
 
+Often, functions don't have roots in closed form.  Nevertheless, it's
+quite easy to compute a solution numerically, to arbitrary precision:
+
+@cindex fsolve
+@example
+> Digits=50:
+> fsolve(cos(x)-x,x,0,2);
+0.7390851332151606416553120876738734040134117589007574649658
+> f=exp(sin(x))-x:
+> X=fsolve(f,x,-10,10);
+2.2191071489137460325957851882042901681753665565320678854155
+> subs(f,x==X);
+-6.372367644529809108115521591070847222364418220770475144296E-58
+@end example
+
+Notice how the final result above differs slightly from zero by about
+@math{6*10^(-58)}.  This is because with 50 decimal digits precision the
+root cannot be represented more accurately than @code{X}.  Such
+inaccuracies are to be expected when computing with finite floating
+point values.
+
 If you ever wanted to convert units in C or C++ and found this is
 cumbersome, here is the solution.  Symbolic types can always be used as
 tags for different types of objects.  Converting from wrong units to the
@@ -602,6 +623,17 @@ insane amounts of memory and CPU time.  Feel free to kill them if your
 machine catches fire.  Another quite important intent is to allow people
 to fiddle around with optimization.
 
+By default, the only documentation that will be built is this tutorial
+in @file{.info} format. To build the GiNaC tutorial and reference manual
+in HTML, DVI, PostScript, or PDF formats, use one of
+
+@example
+$ make html
+$ make dvi
+$ make ps
+$ make pdf
+@end example
+
 Generally, the top-level Makefile runs recursively to the
 subdirectories.  It is therefore safe to go into any subdirectory
 (@code{doc/}, @code{ginsh/}, @dots{}) and simply type @code{make}
@@ -637,7 +669,7 @@ All the header files will be installed into @file{@var{PREFIX}/include/ginac/}
 (or @file{@var{INCLUDEDIR}/ginac/}, if specified).
 
 @item
-All documentation (HTML and Postscript) will be stuffed into
+All documentation (info) will be stuffed into
 @file{@var{PREFIX}/share/doc/GiNaC/} (or
 @file{@var{DATADIR}/doc/GiNaC/}, if @var{DATADIR} was specified).
 
@@ -682,6 +714,7 @@ meta-class for storing all mathematical objects.
 * Lists::                        Lists of expressions.
 * Mathematical functions::       Mathematical functions.
 * Relations::                    Equality, Inequality and all that.
+* Integrals::                    Symbolic integrals.
 * Matrices::                     Matrices.
 * Indexed objects::              Handling indexed quantities.
 * Non-commutative objects::      Algebras with non-commutative products.
@@ -1796,7 +1829,7 @@ point number of class @code{numeric} you should call
 wrapped inside an @code{ex}.
 
 
-@node Relations, Matrices, Mathematical functions, Basic Concepts
+@node Relations, Integrals, Mathematical functions, Basic Concepts
 @c    node-name, next, previous, up
 @section Relations
 @cindex @code{relational} (class)
@@ -1822,8 +1855,68 @@ 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)
 
-@node Matrices, Indexed objects, Relations, Basic Concepts
+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 integrant 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)
@@ -1871,7 +1964,8 @@ 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);
+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
@@ -2013,15 +2107,17 @@ 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, and characteristic polynomials:
+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
@@ -2039,7 +2135,8 @@ 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;
+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}
@@ -2090,7 +2187,7 @@ one or more indices.
 
 @end itemize
 
-@strong{Note:} when printing expressions, covariant indices and indices
+@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
@@ -2491,6 +2588,47 @@ of a sum are consistent:
 @}
 @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
 
@@ -2582,7 +2720,7 @@ representation @code{diag(1, 1, 1, ...)}. It is constructed by the function
         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) << endl;
+         * delta_tensor(i, k) * delta_tensor(j, l);
     cout << e.simplify_indexed() << endl;
      // -> B.i.j*A.i.j
 
@@ -2710,7 +2848,8 @@ 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);
+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
@@ -2809,7 +2948,7 @@ 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 commute and will group the factors
+figuring out by itself which objects commutate and will group the factors
 by their class. Consider this example:
 
 @example
@@ -2825,7 +2964,7 @@ by their class. Consider this 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 commute with color objects). The resulting expression is a
+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.
@@ -2843,7 +2982,7 @@ than in other computer algebra systems; they can, for example, automatically
 canonicalize themselves according to rules specified in the implementation
 of the non-commutative classes. The drawback is that to work with other than
 the built-in algebras you have to implement new classes yourself. Symbols
-always commute and it's not possible to construct non-commutative products
+always commutate and it's not possible to construct non-commutative products
 using symbols to represent the algebra elements or generators. User-defined
 functions can, however, be specified as being non-commutative.
 
@@ -2862,16 +3001,16 @@ the header file @file{flags.h}), corresponding to three categories of
 expressions in GiNaC:
 
 @itemize
-@item @code{return_types::commutative}: Commutes with everything. Most GiNaC
+@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 commute
+  @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 commute with any other @code{noncommutative} or
+  category don't commutate with any other @code{noncommutative} or
   @code{noncommutative_composite} expressions.
 @end itemize
 
@@ -2909,11 +3048,18 @@ 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()}
-Clifford algebra elements (also called Dirac gamma matrices, although 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
+@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);
@@ -2922,7 +3068,7 @@ ex dirac_gamma(const ex & mu, unsigned char rl = 0);
 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 commute with each other. The dimension of the index can be 4 or (in
+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.
 
@@ -2933,14 +3079,14 @@ The unity element of a Clifford algebra is constructed by
 ex dirac_ONE(unsigned char rl = 0);
 @end example
 
-@strong{Note:} You must always use @code{dirac_ONE()} when referring to
+@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 commutes with all other
+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
 
@@ -3000,24 +3146,35 @@ for example
 
 @cindex @code{dirac_trace()}
 To calculate the trace of an expression containing strings of Dirac gammas
-you use the function
+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
 
-This function takes the trace of all gammas with the specified representation
-label; gammas with other labels are left standing. The last argument to
+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 usual trace only in @math{D = 4} dimensions.
-In particular, the functional is not cyclic in @math{D != 4} 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
+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
+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
 @cite{The Role of gamma5 in Dimensional Regularization}.
 
 The value of the trace itself is also usually different in 4 and in
-@math{D != 4} dimensions:
+@tex $D \ne 4$
+@end tex
+dimensions:
 
 @example
 @{
@@ -3081,6 +3238,333 @@ You can use this to compare two expressions or for further simplifications:
 @}
 @end example
 
+@cindex @code{clifford_unit()}
+@subsubsection A generic Clifford algebra
+
+A generic Clifford algebra, i.e. a
+@tex
+$2^n$
+@end tex
+dimensional algebra with
+generators 
+@tex $e_k$
+@end tex 
+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, 
+                                bool anticommuting = false);    
+@end example
+
+where @code{mu} should be a @code{varidx} class object indexing the
+generators, an index @code{mu} with a numeric value may be of type
+@code{idx} as well.
+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. Optional parameter @code{rl} allows to distinguish different
+Clifford algebras, which will commute with each other. The last
+optional parameter @code{anticommuting} defines if the anticommuting
+assumption (i.e.
+@tex
+$e_i e_j + e_j e_i = 0$)
+@end tex
+@ifnottex
+e~i e~j + e~j e~i = 0)
+@end ifnottex
+will be used for contraction of Clifford units. If the @code{metric} is
+supplied by a @code{matrix} object, then the value of
+@code{anticommuting} is calculated automatically and the supplied one
+will be ignored. One can overcome this by giving @code{metric} through
+matrix wrapped into an @code{indexed} object.
+
+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.
+@cindex @code{clifford::is_anticommuting()}
+The method @code{clifford::is_anticommuting()} returns the
+@code{anticommuting} property of a unit.
+
+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
+@{
+    ... 
+    varidx nu(symbol("nu"), 4);
+    realsymbol s("s");
+    ex M = diag_matrix(lst(1, -1, 0, s));
+    ex e = clifford_unit(nu, M);
+    ex e0 = e.subs(nu == 0);
+    ex e1 = e.subs(nu == 1);
+    ex e2 = e.subs(nu == 2);
+    ex e3 = e.subs(nu == 3);
+    ...
+@}
+@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, bool anticommuting = false);
+    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, anticommuting)}. The previous code may be rewritten
+with the help of @code{lst_to_clifford()} as follows
+
+@example
+@{
+    ...
+    varidx nu(symbol("nu"), 4);
+    realsymbol s("s");
+    ex M = diag_matrix(lst(1, -1, 0, s));
+    ex e0 = lst_to_clifford(lst(1, 0, 0, 0), nu, M);
+    ex e1 = lst_to_clifford(lst(0, 1, 0, 0), nu, M);
+    ex e2 = lst_to_clifford(lst(0, 0, 1, 0), nu, M);
+    ex e3 = lst_to_clifford(lst(0, 0, 0, 1), nu, M);
+  ...
+@}
+@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, bool anticommuting = false);
+    ex clifford_moebius_map(const ex & M, const ex & v, const ex & G,
+                            unsigned char rl = 0, bool anticommuting = false);
+@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 parameters @code{rl} and @code{anticommuting} are ignored
+even if supplied.  The returned value of this function is a list of
+components of the resulting vector.
+
+@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
+defintion 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
@@ -3096,7 +3580,7 @@ ex color_T(const ex & a, unsigned char rl = 0);
 
 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 commute with each other. The
+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}.
 
@@ -3107,7 +3591,7 @@ The unity element of a color algebra is constructed by
 ex color_ONE(unsigned char rl = 0);
 @end example
 
-@strong{Note:} You must always use @code{color_ONE()} when referring to
+@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,
@@ -3126,6 +3610,11 @@ 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
 
@@ -3175,16 +3664,19 @@ expressions containing color objects:
 @end example
 
 @cindex @code{color_trace()}
-To calculate the trace of an expression containing color objects you use the
-function
+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
 
-This function takes the trace of all color @samp{T} objects with the
-specified representation label; @samp{T}s with other labels are left
-standing. For 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
     ...
@@ -3428,7 +3920,7 @@ table:
 @end cartouche
 
 To determine whether an expression is commutative or non-commutative and if
-so, with which other expressions it would commute, you use the methods
+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.
 
@@ -4541,9 +5033,11 @@ 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
+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)
+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@});
@@ -4704,6 +5198,7 @@ in which case the value of @code{q} is undefined.
 @cindex @code{unit()}
 @cindex @code{content()}
 @cindex @code{primpart()}
+@cindex @code{unitcontprim()}
 
 The methods
 
@@ -4711,17 +5206,29 @@ The methods
 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 product of unit, content, and primitive part is the
-original polynomial.
+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 and LCM
+@subsection GCD, LCM and resultant
 @cindex GCD
 @cindex LCM
 @cindex @code{gcd()}
@@ -4739,7 +5246,8 @@ 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}.
+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>
@@ -4758,6 +5266,38 @@ int main()
 @}
 @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
@@ -5220,6 +5760,12 @@ GiNaC contains the following predefined mathematical functions:
 @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()}
@@ -5293,15 +5839,22 @@ 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{H}, @code{S} and @code{zeta}.
+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{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.
+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
@@ -5349,11 +5902,18 @@ So in comparison to the referenced publications the order of indices and argumen
 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{m} in @code{H}. Since @code{s} will be interpreted as the sequence
-of signs for the corresponding indices @code{m}, it must contain 1 or -1, e.g.
+@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)$.
+$\zeta(\overline{3},4)$
+@end tex
+and
+@code{G(lst(a,b), lst(-1,1), c)} means
+@tex
+$G(a-0\epsilon,b+0\epsilon;c)$.
 @end tex
 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,
@@ -5386,14 +5946,9 @@ Li(@{3,1,3@},@{-x,1,-1@})
 -Li(@{2,1@},@{x,-1@})*log(x)+2*Li(@{3,1@},@{x,-1@})+Li(@{2,2@},@{x,-1@})
 @end example
 
-Every function apart from the multiple polylogarithm @code{Li} can be numerically evaluated for
-arbitrary real or complex arguments. @code{Li} only evaluates if for all arguments
-@tex
-$x_i$ the condition
-@end tex
-@tex
-$x_1x_2\cdots x_i < 1$ holds.
-@end tex
+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;
@@ -5425,6 +5980,9 @@ 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 Conjugation, Solving Linear Systems of Equations, Multiple polylogarithms, Methods and Functions
 @c    node-name, next, previous, up
 @section Complex Conjugation
@@ -5469,7 +6027,8 @@ 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);
+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
@@ -5548,7 +6107,8 @@ 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 << 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
@@ -5639,7 +6199,8 @@ For example, the code snippet
 will print
 
 @example
-    @{(-\ln(\circ))@}+@{(-\gamma_E)@} \circ+@{(\frac@{1@}@{12@} \pi^@{2@})@} \circ^@{2@}+\mathcal@{O@}(\circ^@{3@})
+    @{(-\ln(\circ))@}+@{(-\gamma_E)@} \circ+@{(\frac@{1@}@{12@} \pi^@{2@})@} \circ^@{2@}
+    +\mathcal@{O@}(\circ^@{3@})
 @end example
 
 @cindex @code{index_dimensions}
@@ -6273,6 +6834,24 @@ print_func<C>(<C++ function>)
 
 option which is explained in the next section.
 
+@subsection Functions with a variable number of arguments
+
+The @code{DECLARE_FUNCTION} and @code{REGISTER_FUNCTION} macros define
+functions with a fixed number of arguments. Sometimes, though, you may need
+to have a function that accepts a variable number of expressions. One way to
+accomplish this is to pass variable-length lists as arguments. The
+@code{Li()} function uses this method for multiple polylogarithms.
+
+It is also possible to define functions that accept a different number of
+parameters under the same function name, such as the @code{psi()} function
+which can be called either as @code{psi(z)} (the digamma function) or as
+@code{psi(n, z)} (polygamma functions). These are actually two different
+functions in GiNaC that, however, have the same name. Defining such
+functions is not possible with the macros but requires manually fiddling
+with GiNaC internals. If you are interested, please consult the GiNaC source
+code for the @code{psi()} function (@file{inifcns.h} and
+@file{inifcns_gamma.cpp}).
+
 
 @node Printing, Structures, Symbolic functions, Extending GiNaC
 @c    node-name, next, previous, up
@@ -6449,7 +7028,8 @@ int main()
     // our own one
     set_print_func<power, print_latex>(my_print_power_as_latex);
 
-    // this prints "-1+@{@{(y+x)@}@}\uparrow@{2@}-3 \frac@{@{x@}\uparrow@{3@}@}@{@{y@}\uparrow@{2@}@}"
+    // this prints "-1+@{@{(y+x)@}@}\uparrow@{2@}-3 \frac@{@{x@}\uparrow@{3@}@}@{@{y@}
+    //              \uparrow@{2@}@}"
     cout << e << endl;
 @}
 @end example
@@ -6793,7 +7373,8 @@ inline bool operator==(const sprod_s & lhs, const sprod_s & rhs)
 
 inline bool operator<(const sprod_s & lhs, const sprod_s & rhs)
 @{
-    return lhs.left.compare(rhs.left) < 0 ? true : lhs.right.compare(rhs.right) < 0;
+    return lhs.left.compare(rhs.left) < 0
+           ? true : lhs.right.compare(rhs.right) < 0;
 @}
 @end example
 
@@ -7117,7 +7698,7 @@ constructor.
 by GiNaC to establish a canonical sort order for terms. It returns 0, +1 or
 -1, depending on the relative order of this object and the @code{other}
 object. If it returns 0, the objects are considered equal.
-@strong{Note:} This has nothing to do with the (numeric) ordering
+@strong{Please notice:} This has nothing to do with the (numeric) ordering
 relationship expressed by @code{<}, @code{>=} etc (which cannot be defined
 for non-numeric classes). For example, @code{numeric(1).compare_same_type(numeric(2))}
 may return +1 even though 1 is clearly smaller than 2. Every GiNaC class
@@ -7851,7 +8432,8 @@ For packages configured using GNU automake, GiNaC also provides
 a macro to automate the process of checking for GiNaC.
 
 @example
-AM_PATH_GINAC([@var{MINIMUM-VERSION}, [@var{ACTION-IF-FOUND} [, @var{ACTION-IF-NOT-FOUND}]]])
+AM_PATH_GINAC([@var{MINIMUM-VERSION}, [@var{ACTION-IF-FOUND}
+              [, @var{ACTION-IF-NOT-FOUND}]]])
 @end example
 
 This macro: