]> www.ginac.de Git - ginac.git/blobdiff - doc/tutorial/ginac.texi
- Fixed typos and layout problems.
[ginac.git] / doc / tutorial / ginac.texi
index 602e2fc542521b18090f9d80501e3a60b7af6bbe..51341f4dd2f2f16075b676378d9f091ebec35cfd 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-2003 Johannes Gutenberg University Mainz, Germany
+Copyright (C) 1999-2007 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
@@ -47,12 +47,11 @@ notice identical to this one.
 @title GiNaC @value{VERSION}
 @subtitle An open framework for symbolic computation within the C++ programming language
 @subtitle @value{UPDATED}
-@author The GiNaC Group:
-@author Christian Bauer, Alexander Frink, Richard Kreckel
+@author @uref{http://www.ginac.de}
 
 @page
 @vskip 0pt plus 1filll
-Copyright @copyright{} 1999-2003 Johannes Gutenberg University Mainz, Germany
+Copyright @copyright{} 1999-2007 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
@@ -79,20 +78,20 @@ framework for symbolic computation within the C++ programming language.
 
 @menu
 * Introduction::                 GiNaC's purpose.
-* A Tour of GiNaC::              A quick tour of the library.
+* A tour of GiNaC::              A quick tour of the library.
 * Installation::                 How to install the package.
-* Basic Concepts::               Description of fundamental classes.
-* Methods and Functions::        Algorithms for symbolic manipulations.
+* Basic concepts::               Description of fundamental classes.
+* Methods and functions::        Algorithms for symbolic manipulations.
 * Extending GiNaC::              How to extend the library.
-* A Comparison With Other CAS::  Compares GiNaC to traditional CAS.
-* Internal Structures::          Description of some internal structures.
-* Package Tools::                Configuring packages to work with GiNaC.
+* A comparison with other CAS::  Compares GiNaC to traditional CAS.
+* Internal structures::          Description of some internal structures.
+* Package tools::                Configuring packages to work with GiNaC.
 * Bibliography::
-* Concept Index::
+* Concept index::
 @end menu
 
 
-@node Introduction, A Tour of GiNaC, Top, Top
+@node Introduction, A tour of GiNaC, Top, Top
 @c    node-name, next, previous, up
 @chapter Introduction
 @cindex history of GiNaC
@@ -135,7 +134,7 @@ the near future.
 
 @section License
 The GiNaC framework for symbolic computation within the C++ programming
-language is Copyright @copyright{} 1999-2003 Johannes Gutenberg
+language is Copyright @copyright{} 1999-2007 Johannes Gutenberg
 University Mainz, Germany.
 
 This program is free software; you can redistribute it and/or
@@ -150,11 +149,11 @@ 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
+@node A tour of GiNaC, How to use it from within C++, Introduction, Top
 @c    node-name, next, previous, up
 @chapter A Tour of GiNaC
 
@@ -168,7 +167,7 @@ leaves many open questions.
 @end menu
 
 
-@node How to use it from within C++, What it can do for you, A Tour of GiNaC, A Tour of GiNaC
+@node How to use it from within C++, What it can do for you, A tour of GiNaC, A tour of GiNaC
 @c    node-name, next, previous, up
 @section How to use it from within C++
 
@@ -206,7 +205,7 @@ $ ./hello
 355687428096000*x*y+20922789888000*y^2+6402373705728000*x^2
 @end example
 
-(@xref{Package Tools}, for tools that help you when creating a software
+(@xref{Package tools}, for tools that help you when creating a software
 package that uses GiNaC.)
 
 @cindex Hermite polynomial
@@ -256,7 +255,7 @@ the @command{ginsh}, a simple GiNaC interactive shell that provides a
 convenient window into GiNaC's capabilities.
 
 
-@node What it can do for you, Installation, How to use it from within C++, A Tour of GiNaC
+@node What it can do for you, Installation, How to use it from within C++, A tour of GiNaC
 @c    node-name, next, previous, up
 @section What it can do for you
 
@@ -417,6 +416,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
@@ -458,14 +478,13 @@ met.  First of all, you need to have a C++-compiler adhering to the
 ANSI-standard @cite{ISO/IEC 14882:1998(E)}.  We used GCC for development
 so if you have a different compiler you are on your own.  For the
 configuration to succeed you need a Posix compliant shell installed in
-@file{/bin/sh}, GNU @command{bash} is fine.  Perl is needed by the built
-process as well, since some of the source files are automatically
-generated by Perl scripts.  Last but not least, Bruno Haible's library
-CLN is extensively used and needs to be installed on your system.
-Please get it either from @uref{ftp://ftp.santafe.edu/pub/gnu/}, from
-@uref{ftp://ftpthep.physik.uni-mainz.de/pub/gnu/, GiNaC's FTP site} or
-from @uref{ftp://ftp.ilog.fr/pub/Users/haible/gnu/, Bruno Haible's FTP
-site} (it is covered by GPL) and install it prior to trying to install
+@file{/bin/sh}, GNU @command{bash} is fine. The pkg-config utility is
+required for the configuration, it can be downloaded from
+@uref{http://pkg-config.freedesktop.org}.
+Last but not least, the CLN library
+is used extensively and needs to be installed on your system.
+Please get it from @uref{ftp://ftpthep.physik.uni-mainz.de/pub/gnu/}
+(it is covered by GPL) and install it prior to trying to install
 GiNaC.  The configure script checks if it can find it and if it cannot
 it will refuse to continue.
 
@@ -535,7 +554,7 @@ warning levels.  If omitted, it defaults to @option{-g
 the file @file{configure.ac}.  It is only distributed in packaged
 releases of GiNaC.  If you got the naked sources, e.g. from CVS, you
 must generate @command{configure} along with the various
-@file{Makefile.in} by using the @command{autogen.sh} script.  This will
+@file{Makefile.in} by using the @command{autoreconf} utility.  This will
 require a fair amount of support from your local toolchain, though.}
 
 The whole process is illustrated in the following two
@@ -602,13 +621,24 @@ 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}
 @var{target} there in case something went wrong.
 
 
-@node Installing GiNaC, Basic Concepts, Building GiNaC, Installation
+@node Installing GiNaC, Basic concepts, Building GiNaC, Installation
 @c    node-name, next, previous, up
 @section Installing GiNaC
 @cindex installation
@@ -637,7 +667,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).
 
@@ -661,9 +691,9 @@ do it by hand since you now know where all the files went during
 installation.}.
 
 
-@node Basic Concepts, Expressions, Installing GiNaC, Top
+@node Basic concepts, Expressions, Installing GiNaC, Top
 @c    node-name, next, previous, up
-@chapter Basic Concepts
+@chapter Basic concepts
 
 This chapter will describe the different fundamental objects that can be
 handled by GiNaC.  But before doing so, it is worthwhile introducing you
@@ -674,7 +704,7 @@ meta-class for storing all mathematical objects.
 * Expressions::                  The fundamental GiNaC class.
 * Automatic evaluation::         Evaluation and canonicalization.
 * Error handling::               How the library reports errors.
-* The Class Hierarchy::          Overview of GiNaC's classes.
+* The class hierarchy::          Overview of GiNaC's classes.
 * Symbols::                      Symbolic objects.
 * Numbers::                      Numerical objects.
 * Constants::                    Pre-defined constants.
@@ -682,13 +712,15 @@ 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.
+* Hash maps::                    A faster alternative to std::map<>.
 @end menu
 
 
-@node Expressions, Automatic evaluation, Basic Concepts, Basic Concepts
+@node Expressions, Automatic evaluation, Basic concepts, Basic concepts
 @c    node-name, next, previous, up
 @section Expressions
 @cindex expression (class @code{ex})
@@ -710,7 +742,7 @@ ex MyEx5 = MyEx4 + 1;               // similar to above
 
 Expressions are handles to other more fundamental objects, that often
 contain other expressions thus creating a tree of expressions
-(@xref{Internal Structures}, for particular examples).  Most methods on
+(@xref{Internal structures}, for particular examples).  Most methods on
 @code{ex} therefore run top-down through such an expression tree.  For
 example, the method @code{has()} scans recursively for occurrences of
 something inside an expression.  Thus, if you have declared @code{MyEx4}
@@ -721,8 +753,27 @@ The next sections will outline the general picture of GiNaC's class
 hierarchy and describe the classes of objects that are handled by
 @code{ex}.
 
+@subsection Note: Expressions and STL containers
+
+GiNaC expressions (@code{ex} objects) have value semantics (they can be
+assigned, reassigned and copied like integral types) but the operator
+@code{<} doesn't provide a well-defined ordering on them. In STL-speak,
+expressions are @samp{Assignable} but not @samp{LessThanComparable}.
+
+This implies that in order to use expressions in sorted containers such as
+@code{std::map<>} and @code{std::set<>} you have to supply a suitable
+comparison predicate. GiNaC provides such a predicate, called
+@code{ex_is_less}. For example, a set of expressions should be defined
+as @code{std::set<ex, ex_is_less>}.
+
+Unsorted containers such as @code{std::vector<>} and @code{std::list<>}
+don't pose a problem. A @code{std::vector<ex>} works as expected.
 
-@node Automatic evaluation, Error handling, Expressions, Basic Concepts
+@xref{Information about expressions}, for more about comparing and ordering
+expressions.
+
+
+@node Automatic evaluation, Error handling, Expressions, Basic concepts
 @c    node-name, next, previous, up
 @section Automatic evaluation and canonicalization of expressions
 @cindex evaluation
@@ -742,7 +793,13 @@ evaluation}. GiNaC only performs transformations that are
 
 @itemize @bullet
 @item
-at most of complexity @math{O(n log n)}
+at most of complexity
+@tex
+$O(n\log n)$
+@end tex
+@ifnottex
+@math{O(n log n)}
+@end ifnottex
 @item
 algebraically correct, possibly except for a set of measure zero (e.g.
 @math{x/x} is transformed to @math{1} although this is incorrect for @math{x=0})
@@ -756,7 +813,7 @@ behave in an entirely obvious way at first glance:
 The terms of sums and products (and some other things like the arguments of
 symmetric functions, the indices of symmetric tensors etc.) are re-ordered
 into a canonical form that is deterministic, but not lexicographical or in
-any other way easily guessable (it almost always depends on the number and
+any other way easy to guess (it almost always depends on the number and
 order of the symbols you define). However, constructing the same expression
 twice, either implicitly or explicitly, will always result in the same
 canonical form.
@@ -790,7 +847,7 @@ transform expressions, like @code{subs()} or @code{normal()}, automatically
 re-evaluate their results.
 
 
-@node Error handling, The Class Hierarchy, Automatic evaluation, Basic Concepts
+@node Error handling, The class hierarchy, Automatic evaluation, Basic concepts
 @c    node-name, next, previous, up
 @section Error handling
 @cindex exceptions
@@ -814,7 +871,7 @@ int pole_error::degree() const;
 that returns the order of the singularity (or 0 when the pole is
 logarithmic or the order is undefined).
 
-When using GiNaC it is useful to arrange for exceptions to be catched in
+When using GiNaC it is useful to arrange for exceptions to be caught in
 the main program even if you don't want to do any special error handling.
 Otherwise whenever an error occurs in GiNaC, it will be delegated to the
 default exception handler of your C++ compiler's run-time system which
@@ -846,9 +903,9 @@ int main()
 @end example
 
 
-@node The Class Hierarchy, Symbols, Error handling, Basic Concepts
+@node The class hierarchy, Symbols, Error handling, Basic concepts
 @c    node-name, next, previous, up
-@section The Class Hierarchy
+@section The class hierarchy
 
 GiNaC's class hierarchy consists of several classes representing
 mathematical objects, all of which (except for @code{ex} and some
@@ -872,7 +929,7 @@ features.  An example is @code{expairseq}, a container for a sequence of
 pairs each consisting of one expression and a number (@code{numeric}).
 What @emph{is} visible to the user are the derived classes @code{add}
 and @code{mul}, representing sums and products.  @xref{Internal
-Structures}, where these two classes are described in more detail.  The
+structures}, where these two classes are described in more detail.  The
 following table shortly summarizes what kinds of mathematical objects
 are stored in the different classes:
 
@@ -899,7 +956,13 @@ $\sqrt{2}$
 @end ifnottex
 @dots{}
 @item @code{pseries} @tab Power Series, e.g. @math{x-1/6*x^3+1/120*x^5+O(x^7)}
-@item @code{function} @tab A symbolic function like @math{sin(2*x)}
+@item @code{function} @tab A symbolic function like
+@tex
+$\sin 2x$
+@end tex
+@ifnottex
+@math{sin(2*x)}
+@end ifnottex
 @item @code{lst} @tab Lists of expressions @{@math{x}, @math{2*y}, @math{3+z}@}
 @item @code{matrix} @tab @math{m}x@math{n} matrices of expressions
 @item @code{relational} @tab A relation like the identity @math{x}@code{==}@math{y}
@@ -914,46 +977,197 @@ $\sqrt{2}$
 @end cartouche
 
 
-@node Symbols, Numbers, The Class Hierarchy, Basic Concepts
+@node Symbols, Numbers, The class hierarchy, Basic concepts
 @c    node-name, next, previous, up
 @section Symbols
 @cindex @code{symbol} (class)
 @cindex hierarchy of classes
 
 @cindex atom
-Symbols are for symbolic manipulation what atoms are for chemistry.  You
-can declare objects of class @code{symbol} as any other object simply by
-saying @code{symbol x,y;}.  There is, however, a catch in here having to
-do with the fact that C++ is a compiled language.  The information about
-the symbol's name is thrown away by the compiler but at a later stage
-you may want to print expressions holding your symbols.  In order to
-avoid confusion GiNaC's symbols are able to know their own name.  This
-is accomplished by declaring its name for output at construction time in
-the fashion @code{symbol x("x");}.  If you declare a symbol using the
-default constructor (i.e. without string argument) the system will deal
-out a unique name.  That name may not be suitable for printing but for
-internal routines when no output is desired it is often enough.  We'll
-come across examples of such symbols later in this tutorial.
-
-This implies that the strings passed to symbols at construction time may
-not be used for comparing two of them.  It is perfectly legitimate to
-write @code{symbol x("x"),y("x");} but it is likely to lead into
-trouble.  Here, @code{x} and @code{y} are different symbols and
-statements like @code{x-y} will not be simplified to zero although the
-output @code{x-x} looks funny.  Such output may also occur when there
-are two different symbols in two scopes, for instance when you call a
-function that declares a symbol with a name already existent in a symbol
-in the calling function.  Again, comparing them (using @code{operator==}
-for instance) will always reveal their difference.  Watch out, please.
+Symbolic indeterminates, or @dfn{symbols} for short, are for symbolic
+manipulation what atoms are for chemistry.
 
-@cindex @code{subs()}
-Although symbols can be assigned expressions for internal reasons, you
-should not do it (and we are not going to tell you how it is done).  If
-you want to replace a symbol with something else in an expression, you
-can use the expression's @code{.subs()} method (@pxref{Substituting Expressions}).
+A typical symbol definition looks like this:
+@example
+symbol x("x");
+@end example
+
+This definition actually contains three very different things:
+@itemize
+@item a C++ variable named @code{x}
+@item a @code{symbol} object stored in this C++ variable; this object
+  represents the symbol in a GiNaC expression
+@item the string @code{"x"} which is the name of the symbol, used (almost)
+  exclusively for printing expressions holding the symbol
+@end itemize
+
+Symbols have an explicit name, supplied as a string during construction,
+because in C++, variable names can't be used as values, and the C++ compiler
+throws them away during compilation.
+
+It is possible to omit the symbol name in the definition:
+@example
+symbol x;
+@end example
+
+In this case, GiNaC will assign the symbol an internal, unique name of the
+form @code{symbolNNN}. This won't affect the usability of the symbol but
+the output of your calculations will become more readable if you give your
+symbols sensible names (for intermediate expressions that are only used
+internally such anonymous symbols can be quite useful, however).
+
+Now, here is one important property of GiNaC that differentiates it from
+other computer algebra programs you may have used: GiNaC does @emph{not} use
+the names of symbols to tell them apart, but a (hidden) serial number that
+is unique for each newly created @code{symbol} object. If you want to use
+one and the same symbol in different places in your program, you must only
+create one @code{symbol} object and pass that around. If you create another
+symbol, even if it has the same name, GiNaC will treat it as a different
+indeterminate.
+
+Observe:
+@example
+ex f(int n)
+@{
+    symbol x("x");
+    return pow(x, n);
+@}
+
+int main()
+@{
+    symbol x("x");
+    ex e = f(6);
+
+    cout << e << endl;
+     // prints "x^6" which looks right, but...
+
+    cout << e.degree(x) << endl;
+     // ...this doesn't work. The symbol "x" here is different from the one
+     // in f() and in the expression returned by f(). Consequently, it
+     // prints "0".
+@}
+@end example
+
+One possibility to ensure that @code{f()} and @code{main()} use the same
+symbol is to pass the symbol as an argument to @code{f()}:
+@example
+ex f(int n, const ex & x)
+@{
+    return pow(x, n);
+@}
+
+int main()
+@{
+    symbol x("x");
+
+    // Now, f() uses the same symbol.
+    ex e = f(6, x);
+
+    cout << e.degree(x) << endl;
+     // prints "6", as expected
+@}
+@end example
+
+Another possibility would be to define a global symbol @code{x} that is used
+by both @code{f()} and @code{main()}. If you are using global symbols and
+multiple compilation units you must take special care, however. Suppose
+that you have a header file @file{globals.h} in your program that defines
+a @code{symbol x("x");}. In this case, every unit that includes
+@file{globals.h} would also get its own definition of @code{x} (because
+header files are just inlined into the source code by the C++ preprocessor),
+and hence you would again end up with multiple equally-named, but different,
+symbols. Instead, the @file{globals.h} header should only contain a
+@emph{declaration} like @code{extern symbol x;}, with the definition of
+@code{x} moved into a C++ source file such as @file{globals.cpp}.
+
+A different approach to ensuring that symbols used in different parts of
+your program are identical is to create them with a @emph{factory} function
+like this one:
+@example
+const symbol & get_symbol(const string & s)
+@{
+    static map<string, symbol> directory;
+    map<string, symbol>::iterator i = directory.find(s);
+    if (i != directory.end())
+        return i->second;
+    else
+        return directory.insert(make_pair(s, symbol(s))).first->second;
+@}
+@end example
+
+This function returns one newly constructed symbol for each name that is
+passed in, and it returns the same symbol when called multiple times with
+the same name. Using this symbol factory, we can rewrite our example like
+this:
+@example
+ex f(int n)
+@{
+    return pow(get_symbol("x"), n);
+@}
+
+int main()
+@{
+    ex e = f(6);
+
+    // Both calls of get_symbol("x") yield the same symbol.
+    cout << e.degree(get_symbol("x")) << endl;
+     // prints "6"
+@}
+@end example
+
+Instead of creating symbols from strings we could also have
+@code{get_symbol()} take, for example, an integer number as its argument.
+In this case, we would probably want to give the generated symbols names
+that include this number, which can be accomplished with the help of an
+@code{ostringstream}.
+
+In general, if you're getting weird results from GiNaC such as an expression
+@samp{x-x} that is not simplified to zero, you should check your symbol
+definitions.
+
+As we said, the names of symbols primarily serve for purposes of expression
+output. But there are actually two instances where GiNaC uses the names for
+identifying symbols: When constructing an expression from a string, and when
+recreating an expression from an archive (@pxref{Input/output}).
 
+In addition to its name, a symbol may contain a special string that is used
+in LaTeX output:
+@example
+symbol x("x", "\\Box");
+@end example
 
-@node Numbers, Constants, Symbols, Basic Concepts
+This creates a symbol that is printed as "@code{x}" in normal output, but
+as "@code{\Box}" in LaTeX code (@xref{Input/output}, for more
+information about the different output formats of expressions in GiNaC).
+GiNaC automatically creates proper LaTeX code for symbols having names of
+greek letters (@samp{alpha}, @samp{mu}, etc.).
+
+@cindex @code{subs()}
+Symbols in GiNaC can't be assigned values. If you need to store results of
+calculations and give them a name, use C++ variables of type @code{ex}.
+If you want to replace a symbol in an expression with something else, you
+can invoke the expression's @code{.subs()} method
+(@pxref{Substituting expressions}).
+
+@cindex @code{realsymbol()}
+By default, symbols are expected to stand in for complex values, i.e. they live
+in the complex domain.  As a consequence, operations like complex conjugation,
+for example (@pxref{Complex expressions}), do @emph{not} evaluate if applied
+to such symbols. Likewise @code{log(exp(x))} does not evaluate to @code{x},
+because of the unknown imaginary part of @code{x}.
+On the other hand, if you are sure that your symbols will hold only real
+values, you would like to have such functions evaluated. Therefore GiNaC
+allows you to specify
+the domain of the symbol. Instead of @code{symbol x("x");} you can write
+@code{realsymbol x("x");} to tell GiNaC that @code{x} stands in for real values.
+
+@cindex @code{possymbol()}
+Furthermore, it is also possible to declare a symbol as positive. This will,
+for instance, enable the automatic simplification of @code{abs(x)} into 
+@code{x}. This is done by declaring the symbol as @code{possymbol x("x");}.
+
+
+@node Numbers, Constants, Symbols, Basic concepts
 @c    node-name, next, previous, up
 @section Numbers
 @cindex @code{numeric} (class)
@@ -1121,7 +1335,7 @@ int main()
 @end example
 
 Note that the variable @code{answer} is constructed here as an integer
-by @code{numeric}'s copy constructor but in an intermediate step it
+by @code{numeric}'s copy constructor, but in an intermediate step it
 holds a rational number represented as integer numerator and integer
 denominator.  When multiplied by 10, the denominator becomes unity and
 the result is automatically converted to a pure integer again.
@@ -1139,6 +1353,8 @@ can be applied is listed in the following table.
 @tab @dots{}equal to zero
 @item @code{.is_positive()}
 @tab @dots{}not complex and greater than 0
+@item @code{.is_negative()}
+@tab @dots{}not complex and smaller than 0
 @item @code{.is_integer()}
 @tab @dots{}a (non-complex) integer
 @item @code{.is_pos_integer()}
@@ -1162,8 +1378,150 @@ can be applied is listed in the following table.
 @end multitable
 @end cartouche
 
+@page
+
+@subsection Numeric functions
+
+The following functions can be applied to @code{numeric} objects and will be
+evaluated immediately:
+
+@cartouche
+@multitable @columnfractions .30 .70
+@item @strong{Name} @tab @strong{Function}
+@item @code{inverse(z)}
+@tab returns @math{1/z}
+@cindex @code{inverse()} (numeric)
+@item @code{pow(a, b)}
+@tab exponentiation @math{a^b}
+@item @code{abs(z)}
+@tab absolute value
+@item @code{real(z)}
+@tab real part
+@cindex @code{real()}
+@item @code{imag(z)}
+@tab imaginary part
+@cindex @code{imag()}
+@item @code{csgn(z)}
+@tab complex sign (returns an @code{int})
+@item @code{step(x)}
+@tab step function (returns an @code{numeric})
+@item @code{numer(z)}
+@tab numerator of rational or complex rational number
+@item @code{denom(z)}
+@tab denominator of rational or complex rational number
+@item @code{sqrt(z)}
+@tab square root
+@item @code{isqrt(n)}
+@tab integer square root
+@cindex @code{isqrt()}
+@item @code{sin(z)}
+@tab sine
+@item @code{cos(z)}
+@tab cosine
+@item @code{tan(z)}
+@tab tangent
+@item @code{asin(z)}
+@tab inverse sine
+@item @code{acos(z)}
+@tab inverse cosine
+@item @code{atan(z)}
+@tab inverse tangent
+@item @code{atan(y, x)}
+@tab inverse tangent with two arguments
+@item @code{sinh(z)}
+@tab hyperbolic sine
+@item @code{cosh(z)}
+@tab hyperbolic cosine
+@item @code{tanh(z)}
+@tab hyperbolic tangent
+@item @code{asinh(z)}
+@tab inverse hyperbolic sine
+@item @code{acosh(z)}
+@tab inverse hyperbolic cosine
+@item @code{atanh(z)}
+@tab inverse hyperbolic tangent
+@item @code{exp(z)}
+@tab exponential function
+@item @code{log(z)}
+@tab natural logarithm
+@item @code{Li2(z)}
+@tab dilogarithm
+@item @code{zeta(z)}
+@tab Riemann's zeta function
+@item @code{tgamma(z)}
+@tab gamma function
+@item @code{lgamma(z)}
+@tab logarithm of gamma function
+@item @code{psi(z)}
+@tab psi (digamma) function
+@item @code{psi(n, z)}
+@tab derivatives of psi function (polygamma functions)
+@item @code{factorial(n)}
+@tab factorial function @math{n!}
+@item @code{doublefactorial(n)}
+@tab double factorial function @math{n!!}
+@cindex @code{doublefactorial()}
+@item @code{binomial(n, k)}
+@tab binomial coefficients
+@item @code{bernoulli(n)}
+@tab Bernoulli numbers
+@cindex @code{bernoulli()}
+@item @code{fibonacci(n)}
+@tab Fibonacci numbers
+@cindex @code{fibonacci()}
+@item @code{mod(a, b)}
+@tab modulus in positive representation (in the range @code{[0, abs(b)-1]} with the sign of b, or zero)
+@cindex @code{mod()}
+@item @code{smod(a, b)}
+@tab modulus in symmetric representation (in the range @code{[-iquo(abs(b)-1, 2), iquo(abs(b), 2)]})
+@cindex @code{smod()}
+@item @code{irem(a, b)}
+@tab integer remainder (has the sign of @math{a}, or is zero)
+@cindex @code{irem()}
+@item @code{irem(a, b, q)}
+@tab integer remainder and quotient, @code{irem(a, b, q) == a-q*b}
+@item @code{iquo(a, b)}
+@tab integer quotient
+@cindex @code{iquo()}
+@item @code{iquo(a, b, r)}
+@tab integer quotient and remainder, @code{r == a-iquo(a, b)*b}
+@item @code{gcd(a, b)}
+@tab greatest common divisor
+@item @code{lcm(a, b)}
+@tab least common multiple
+@end multitable
+@end cartouche
+
+Most of these functions are also available as symbolic functions that can be
+used in expressions (@pxref{Mathematical functions}) or, like @code{gcd()},
+as polynomial algorithms.
+
+@subsection Converting numbers
+
+Sometimes it is desirable to convert a @code{numeric} object back to a
+built-in arithmetic type (@code{int}, @code{double}, etc.). The @code{numeric}
+class provides a couple of methods for this purpose:
+
+@cindex @code{to_int()}
+@cindex @code{to_long()}
+@cindex @code{to_double()}
+@cindex @code{to_cl_N()}
+@example
+int numeric::to_int() const;
+long numeric::to_long() const;
+double numeric::to_double() const;
+cln::cl_N numeric::to_cl_N() const;
+@end example
+
+@code{to_int()} and @code{to_long()} only work when the number they are
+applied on is an exact integer. Otherwise the program will halt with a
+message like @samp{Not a 32-bit integer}. @code{to_double()} applied on a
+rational number will return a floating-point approximation. Both
+@code{to_int()/to_long()} and @code{to_double()} discard the imaginary
+part of complex numbers.
+
 
-@node Constants, Fundamental containers, Numbers, Basic Concepts
+@node Constants, Fundamental containers, Numbers, Basic concepts
 @c    node-name, next, previous, up
 @section Constants
 @cindex @code{constant} (class)
@@ -1178,7 +1536,7 @@ specific number when the method @code{.evalf()} is called.
 The predefined known constants are:
 
 @cartouche
-@multitable @columnfractions .14 .30 .56
+@multitable @columnfractions .14 .32 .54
 @item @strong{Name} @tab @strong{Common Name} @tab @strong{Numerical Value (to 35 digits)}
 @item @code{Pi}
 @tab Archimedes' constant
@@ -1193,7 +1551,7 @@ The predefined known constants are:
 @end cartouche
 
 
-@node Fundamental containers, Lists, Constants, Basic Concepts
+@node Fundamental containers, Lists, Constants, Basic concepts
 @c    node-name, next, previous, up
 @section Sums, products and powers
 @cindex polynomial
@@ -1264,7 +1622,7 @@ and safe simplifications are carried out like transforming
 @code{3*x+4-x} to @code{2*x+4}.
 
 
-@node Lists, Mathematical functions, Fundamental containers, Basic Concepts
+@node Lists, Mathematical functions, Fundamental containers, Basic concepts
 @c    node-name, next, previous, up
 @section Lists of expressions
 @cindex @code{lst} (class)
@@ -1280,17 +1638,30 @@ and safe simplifications are carried out like transforming
 The GiNaC class @code{lst} serves for holding a @dfn{list} of arbitrary
 expressions. They are not as ubiquitous as in many other computer algebra
 packages, but are sometimes used to supply a variable number of arguments of
-the same type to GiNaC methods such as @code{subs()} and @code{to_rational()},
-so you should have a basic understanding of them.
+the same type to GiNaC methods such as @code{subs()} and some @code{matrix}
+constructors, so you should have a basic understanding of them.
 
-Lists of up to 16 expressions can be directly constructed from single
+Lists can be constructed by assigning a comma-separated sequence of
 expressions:
 
 @example
 @{
     symbol x("x"), y("y");
-    lst l(x, 2, y, x+y);
-    // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y'
+    lst l;
+    l = x, 2, y, x+y;
+    // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y',
+    // in that order
+    ...
+@end example
+
+There are also constructors that allow direct creation of lists of up to
+16 expressions, which is often more convenient but slightly less efficient:
+
+@example
+    ...
+    // This produces the same list 'l' as above:
+    // lst l(x, 2, y, x+y);
+    // lst l = lst(x, 2, y, x+y);
     ...
 @end example
 
@@ -1345,10 +1716,10 @@ the C++ standard library:
 @example
     ...
     // print the elements of the list (requires #include <iterator>)
-    copy(l.begin(), l.end(), ostream_iterator<ex>(cout, "\n"));
+    std::copy(l.begin(), l.end(), ostream_iterator<ex>(cout, "\n"));
 
     // sum up the elements of the list (requires #include <numeric>)
-    ex sum = accumulate(l.begin(), l.end(), ex(0));
+    ex sum = std::accumulate(l.begin(), l.end(), ex(0));
     cout << sum << endl;  // prints '2+2*x+2*y'
     ...
 @end example
@@ -1395,8 +1766,9 @@ You can bring the elements of a list into a canonical order with @code{sort()}:
 
 @example
     ...
-    lst l1(x, 2, y, x+y);
-    lst l2(2, x+y, x, y);
+    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
@@ -1408,13 +1780,14 @@ elements with @code{unique()}:
 
 @example
     ...
-    lst l3(x, 2, 2, 2, y, x+y, y+x);
+    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
+@node Mathematical functions, Relations, Lists, Basic concepts
 @c    node-name, next, previous, up
 @section Mathematical functions
 @cindex @code{function} (class)
@@ -1423,7 +1796,7 @@ elements with @code{unique()}:
 
 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).
+(@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
@@ -1466,7 +1839,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)
@@ -1492,8 +1865,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)
+
+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, Relations, Basic Concepts
+@node Matrices, Indexed objects, Integrals, Basic concepts
 @c    node-name, next, previous, up
 @section Matrices
 @cindex @code{matrix} (class)
@@ -1504,33 +1937,84 @@ matrix with @math{m} rows and @math{n} columns are accessed with two
 second one in the range 0@dots{}@math{n-1}.
 
 There are a couple of ways to construct matrices, with or without preset
-elements:
+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
-matrix::matrix(unsigned r, unsigned c);
-matrix::matrix(unsigned r, unsigned c, const lst & l);
-ex lst_to_matrix(const lst & l);
 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
+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 first two functions are @code{matrix} constructors which create a matrix
-with @samp{r} rows and @samp{c} columns. The matrix elements can be
-initialized from a (flat) list of expressions @samp{l}. Otherwise they are
-all set to zero. The @code{lst_to_matrix()} function constructs a matrix
-from a list of lists, each list representing a matrix row. @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.
+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:
@@ -1544,18 +2028,24 @@ 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 of constructing matrices:
+Here are a couple of examples for constructing matrices:
 
 @example
 @{
     symbol a("a"), b("b");
 
     matrix M(2, 2);
-    M(0, 0) = a;
-    M(1, 1) = b;
+    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]]
 
@@ -1573,6 +2063,12 @@ Here are a couple of examples of constructing matrices:
 @}
 @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:
@@ -1592,9 +2088,13 @@ and @math{C}:
 
 @example
 @{
-    matrix A(2, 2, lst(1, 2, 3, 4));
-    matrix B(2, 2, lst(-1, 0, 2, 1));
-    matrix C(2, 2, lst(8, 4, 2, 1));
+    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;
@@ -1655,15 +2155,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 symbol & lambda) const;
+ex matrix::charpoly(const ex & lambda) const;
+unsigned matrix::rank() const;
 @end example
 
 The @samp{algo} argument of @code{determinant()} allows to select
@@ -1675,13 +2177,14 @@ 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()}
+@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;
+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}
@@ -1693,7 +2196,7 @@ 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
+@node Indexed objects, Non-commutative objects, Matrices, Basic concepts
 @c    node-name, next, previous, up
 @section Indexed objects
 
@@ -1732,7 +2235,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
@@ -1842,7 +2345,7 @@ bool idx::is_dim_symbolic();
 
 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.
+about expressions}) returns information about the index value.
 
 @cindex @code{varidx} (class)
 If you need co- and contravariant indices, use the @code{varidx} class:
@@ -1933,7 +2436,7 @@ and the same or opposite variance.
 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}).
+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:
@@ -2133,6 +2636,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
 
@@ -2148,7 +2692,7 @@ ex ex::simplify_indexed(const scalar_products & sp);
 
 that performs some more expensive operations:
 
-@itemize
+@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
@@ -2189,11 +2733,7 @@ arithmetic class, you just pass it to @code{simplify_indexed()}):
 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. After @code{sp.add(A, B, 0)},
-@code{simplify_indexed()} will replace all scalar products of indexed
-objects that have the symbols @code{A} and @code{B} as base expressions
-with the single value 0. The number, type and dimension of the indices
-don't matter; @samp{A~mu~nu*B.mu.nu} would also be replaced by 0.
+taken, and the expression to replace it with.
 
 @cindex @code{expand()}
 The example above also illustrates a feature of the @code{expand()} method:
@@ -2224,7 +2764,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
 
@@ -2352,7 +2892,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
@@ -2393,7 +2934,10 @@ and scalar products):
     symbol x("x"), y("y");
 
     // A is a 2x2 matrix, X is a 2x1 vector
-    matrix A(2, 2, lst(1, 2, 3, 4)), X(2, 1, lst(x, y));
+    matrix A(2, 2), X(2, 1);
+    A = 1, 2,
+        3, 4;
+    X = x, y;
 
     cout << indexed(A, i, i) << endl;
      // -> 5
@@ -2424,7 +2968,7 @@ one form for @samp{F} and explicitly multiply it with a matrix representation
 of the metric tensor.
 
 
-@node Non-commutative objects, Methods and Functions, Indexed objects, Basic Concepts
+@node Non-commutative objects, Hash maps, Indexed objects, Basic concepts
 @c    node-name, next, previous, up
 @section Non-commutative objects
 
@@ -2448,7 +2992,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
@@ -2464,7 +3008,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.
@@ -2481,10 +3025,8 @@ expressions. Also, non-commutative products in GiNaC are more intelligent
 than in other computer algebra systems; they can, for example, automatically
 canonicalize themselves according to rules specified in the implementation
 of the non-commutative classes. The drawback is that to work with other than
-the built-in algebras you have to implement new classes yourself. Symbols
-always commute 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.
+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()}
@@ -2500,24 +3042,24 @@ 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
-@item @code{return_types::commutative}: Commutes with everything. Most 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 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
 
 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 and usually one of the
-constants in @file{tinfos.h}, or derived therefrom.
+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:
 
@@ -2548,11 +3090,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);
@@ -2561,7 +3110,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.
 
@@ -2572,14 +3121,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
 
@@ -2639,24 +3188,41 @@ 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
-@cite{The Role of gamma5 in Dimensional Regularization}.
+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
-@math{D != 4} dimensions:
+@tex $D \ne 4$
+@end tex
+@ifnottex
+@math{D != 4}
+@end ifnottex
+dimensions:
 
 @example
 @{
@@ -2716,57 +3282,380 @@ You can use this to compare two expressions or for further simplifications:
 
     e = canonicalize_clifford(e);
     cout << e << endl;
-     // -> 2*eta~mu~nu
+     // -> 2*ONE*eta~mu~nu
 @}
 @end example
 
+@cindex @code{clifford_unit()}
+@subsubsection A generic Clifford algebra
 
-@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
+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 color_T(const ex & a, unsigned char rl = 0);
+    ex clifford_unit(const ex & mu, const ex & metr, 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 commute with each other. The
-dimension of the index must be exactly 8 and it should be of class @code{idx},
-not @code{varidx}.
+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. 
 
-@cindex @code{color_ONE()}
-The unity element of a color algebra is constructed by
+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 color_ONE(unsigned char rl = 0);
+    ex e = clifford_unit(mu, indexed(M, sy_symm(), i, j));
 @end example
 
-@strong{Note:} 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.
+since this may yield some further automatic simplifications. Again, for a
+metric defined through a @code{matrix} such a symmetry is detected
+automatically. 
 
-@cindex @code{color_d()}
-@cindex @code{color_f()}
-The functions
+Individual generators of a Clifford algebra can be accessed in several
+ways. For example 
 
 @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);
+@{
+    ... 
+    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
 
-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}.
-
-@cindex @code{color_h()}
-There's an additional function
+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);
@@ -2814,16 +3703,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
     ...
@@ -2834,9 +3726,45 @@ standing. For example:
 @end example
 
 
-@node Methods and Functions, Information About Expressions, Non-commutative objects, Top
+@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
+@chapter Methods and functions
 @cindex polynomial
 
 In this chapter the most important algorithms provided by GiNaC will be
@@ -2875,23 +3803,26 @@ method on class @code{ex} and sometimes calling a function cannot be
 avoided.
 
 @menu
-* Information About Expressions::
-* 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.
+* 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.
-* Solving Linear Systems of Equations::
-* Input/Output::                    Input and output of expressions.
+* 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, Substituting Expressions, Methods and Functions, Methods and Functions
+@node Information about expressions, Numerical evaluation, Methods and functions, Methods and functions
 @c    node-name, next, previous, up
 @section Getting information about expressions
 
@@ -2918,7 +3849,7 @@ unsigned ex::return_type_tinfo() const;
 
 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
+class names (@xref{The class hierarchy}, for a list of all classes). For
 example, assuming @code{e} is an @code{ex}:
 
 @example
@@ -2932,7 +3863,7 @@ example, assuming @code{e} is an @code{ex}:
 
 @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,
+(@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
@@ -2962,9 +3893,9 @@ table:
 @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_<numeric>(...)})
+@tab @dots{}a number (same as @code{is_a<numeric>(...)})
 @item @code{real}
-@tab @dots{}a real integer, rational or float (i.e. is not complex)
+@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}
@@ -3027,34 +3958,123 @@ 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.
 
 
 @subsection Accessing subexpressions
-@cindex @code{nops()}
-@cindex @code{op()}
 @cindex container
-@cindex @code{relational} (class)
 
-GiNaC provides the two methods
+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
 
-for accessing the subexpressions in the container-like GiNaC classes like
-@code{add}, @code{mul}, @code{lst}, and @code{function}. @code{nops()}
-determines the number of subexpressions (@samp{operands}) contained, while
-@code{op()} 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.
+@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
 
-The left-hand and right-hand side expressions of objects of class
-@code{relational} (and only of these) can also be accessed 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();
@@ -3086,7 +4106,8 @@ bool ex::is_zero();
 @end example
 
 for checking whether one expression is equal to another, or equal to zero,
-respectively.
+respectively. See also the method @code{ex::is_zero_matrix()}, 
+@pxref{Matrices}. 
 
 
 @subsection Ordering expressions
@@ -3159,7 +4180,47 @@ if @code{*this} sorts before @code{other}, and @math{1} if @code{*this} sorts
 after @code{other}.
 
 
-@node Substituting Expressions, Pattern Matching and Advanced Substitutions, Information About Expressions, Methods and Functions
+@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()}
@@ -3169,6 +4230,7 @@ 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
 
@@ -3192,19 +4254,51 @@ In the first form, @code{subs()} accepts a relational of the form
 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 two lists, one for the objects to be
+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
-@code{subs(lst(x, y), lst(y, x))} to exchange @samp{x} and @samp{y}.
+
+@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 two options available:
+@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.
-@ref{Pattern Matching and Advanced Substitutions}, for more information
-about patterns and algebraic substitutions.
+@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
@@ -3233,7 +4327,7 @@ 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
+@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)
@@ -3270,7 +4364,7 @@ Some examples for patterns:
 
 Notes:
 
-@itemize
+@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
@@ -3367,7 +4461,7 @@ FAIL
 > match((a+b)*(a+c),($1+b)*($1+c));
 @{$1==a@}
 > match((a+b)*(a+c),(a+$1)*(a+$2));
-@{$1==c,$2==b@}
+@{$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
@@ -3508,61 +4602,29 @@ The last example would be written in C++ in this way:
 @}
 @end example
 
-@subsection Algebraic substitutions
-Supplying the @code{subs_options::algebraic} option to @code{subs()}
-enables smarter, algebraic substitutions in products and powers. If you want
-to substitute some factors of a product, you only need to list these factors
-in your pattern. Furthermore, if an (integer) power of some expression occurs
-in your pattern and in the expression that you want the substitution to occur
-in, it can be substituted as many times as possible, without getting negative
-powers.
-
-An example clarifies it all (hopefully):
-
-@example
-cout << (a*a*a*a+b*b*b*b+pow(x+y,4)).subs(wild()*wild()==pow(wild(),3),
-                                        subs_options::algebraic) << endl;
-// --> (y+x)^6+b^6+a^6
-
-cout << ((a+b+c)*(a+b+c)).subs(a+b==x,subs_options::algebraic) << endl;
-// --> (c+b+a)^2
-// Powers and products are smart, but addition is just the same.
-
-cout << ((a+b+c)*(a+b+c)).subs(a+b+wild()==x+wild(), subs_options::algebraic)
-                                                                      << endl;
-// --> (x+c)^2
-// As I said: addition is just the same.
-
-cout << (pow(a,5)*pow(b,7)+2*b).subs(b*b*a==x,subs_options::algebraic) << endl;
-// --> x^3*b*a^2+2*b
-
-cout << (pow(a,-5)*pow(b,-7)+2*b).subs(1/(b*b*a)==x,subs_options::algebraic)
-                                                                       << endl;
-// --> 2*b+x^3*b^(-1)*a^(-2)
-
-cout << (4*x*x*x-2*x*x+5*x-1).subs(x==a,subs_options::algebraic) << endl;
-// --> -1-2*a^2+4*a^3+5*a
-
-cout << (4*x*x*x-2*x*x+5*x-1).subs(pow(x,wild())==pow(a,wild()),
-                                subs_options::algebraic) << endl;
-// --> -1+5*x+4*x^3-2*x^2
-// You should not really need this kind of patterns very often now.
-// But perhaps this it's-not-a-bug-it's-a-feature (c/sh)ould still change.
-
-cout << ex(sin(1+sin(x))).subs(sin(wild())==cos(wild()),
-                                subs_options::algebraic) << endl;
-// --> cos(1+cos(x))
+@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}. 
 
-cout << expand((a*sin(x+y)*sin(x+y)+a*cos(x+y)*cos(x+y)+b)
-        .subs((pow(cos(wild()),2)==1-pow(sin(wild()),2)),
-                                subs_options::algebraic)) << endl;
-// --> b+a
-@end example
+@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
+@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
+@section Applying a function on subexpressions
 @cindex tree traversal
 @cindex @code{map()}
 
@@ -3705,9 +4767,9 @@ argument. You can not use functions like @samp{diff()}, @samp{op()},
 @end example
 
 
-@node Visitors and Tree Traversal, Polynomial Arithmetic, Applying a Function on Subexpressions, Methods and Functions
+@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
+@section Visitors and tree traversal
 @cindex tree traversal
 @cindex @code{visitor} (class)
 @cindex @code{accept()}
@@ -3911,11 +4973,42 @@ lst gather_indices(const ex & e)
 @}
 @end example
 
+Alternatively, you could use pre- or postorder iterators for the tree
+traversal:
 
-@node Polynomial Arithmetic, Rational Expressions, Visitors and Tree Traversal, Methods and Functions
+@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()}
@@ -3941,7 +5034,7 @@ ex ex::expand(unsigned options = 0);
 
 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 easily guessable you should be prepared to see different
+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
@@ -3967,9 +5060,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@});
@@ -4011,8 +5106,9 @@ 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). To extract
-a coefficient with a certain power from an expanded polynomial you use
+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);
@@ -4068,7 +5164,7 @@ constants, functions and indexed objects as well:
 
 @example
 @{
-    symbol a("a"), b("b"), c("c");
+    symbol a("a"), b("b"), c("c"), x("x");
     idx i(symbol("i"), 3);
 
     ex e = pow(sin(x) - cos(x), 4);
@@ -4098,8 +5194,8 @@ constants, functions and indexed objects as well:
 The two functions
 
 @example
-ex quo(const ex & a, const ex & b, const symbol & x);
-ex rem(const ex & a, const ex & b, const symbol & x);
+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
@@ -4108,7 +5204,7 @@ compute the quotient and remainder of univariate polynomials in the variable
 The additional function
 
 @example
-ex prem(const ex & a, const ex & b, const symbol & x);
+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
@@ -4129,24 +5225,37 @@ in which case the value of @code{q} is undefined.
 @cindex @code{unit()}
 @cindex @code{content()}
 @cindex @code{primpart()}
+@cindex @code{unitcontprim()}
 
 The methods
 
 @example
-ex ex::unit(const symbol & x);
-ex ex::content(const symbol & x);
-ex ex::primpart(const symbol & x);
+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
 
-@subsection GCD and LCM
+@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()}
@@ -4164,7 +5273,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>
@@ -4183,6 +5293,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
@@ -4221,7 +5363,7 @@ 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
+@node Rational expressions, Symbolic differentiation, Polynomial arithmetic, Methods and functions
 @c    node-name, next, previous, up
 @section Rational expressions
 
@@ -4293,21 +5435,23 @@ general expressions by using the temporary replacement algorithm described
 above. You do this by calling
 
 @example
-ex ex::to_polynomial(lst &l);
+ex ex::to_polynomial(exmap & m);
+ex ex::to_polynomial(lst & l);
 @end example
 or
 @example
-ex ex::to_rational(lst &l);
+ex ex::to_rational(exmap & m);
+ex ex::to_rational(lst & l);
 @end example
 
-on the expression to be converted. The supplied @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.
+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 betwerrn @code{.to_polynomial()} and @code{.to_rational()}
+The difference between @code{.to_polynomial()} and @code{.to_rational()}
 is probably best illustrated with an example:
 
 @example
@@ -4338,14 +5482,14 @@ The following more useful example will print @samp{sin(x)-cos(x)}:
     ex a = pow(sin(x), 2) - pow(cos(x), 2);
     ex b = sin(x) + cos(x);
     ex q;
-    lst l;
-    divide(a.to_polynomial(l), b.to_polynomial(l), q);
-    cout << q.subs(l) << endl;
+    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
+@node Symbolic differentiation, Series expansion, Rational expressions, Methods and functions
 @c    node-name, next, previous, up
 @section Symbolic differentiation
 @cindex differentiation
@@ -4411,7 +5555,7 @@ When you run it, it produces the sequence @code{1}, @code{-1}, @code{5},
 @code{i} by two since all odd Euler numbers vanish anyways.
 
 
-@node Series Expansion, Symmetrization, Symbolic Differentiation, Methods and Functions
+@node Series expansion, Symmetrization, Symbolic differentiation, Methods and functions
 @c    node-name, next, previous, up
 @section Series expansion
 @cindex @code{series()}
@@ -4460,6 +5604,9 @@ 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
@@ -4524,7 +5671,7 @@ program, it will type out:
 @end example
 
 
-@node Symmetrization, Built-in Functions, Series Expansion, Methods and Functions
+@node Symmetrization, Built-in functions, Series expansion, Methods and functions
 @c    node-name, next, previous, up
 @section Symmetrization
 @cindex @code{symmetrize()}
@@ -4570,10 +5717,13 @@ almost any kind of object (anything that is @code{subs()}able):
 @}
 @end example
 
+@page
 
-@node Built-in Functions, Solving Linear Systems of Equations, Symmetrization, Methods and Functions
+@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:
 
@@ -4583,9 +5733,20 @@ GiNaC contains the following predefined mathematical functions:
 @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{csgn()}
+@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()}
@@ -4634,22 +5795,40 @@ GiNaC contains the following predefined mathematical functions:
 @tab natural logarithm
 @cindex @code{log()}
 @item @code{Li2(x)}
-@tab Dilogarithm
+@tab dilogarithm
 @cindex @code{Li2()}
-@item @code{zeta(x)}
-@tab Riemann's zeta function
+@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(n, x)}
+@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
+@tab gamma function
 @cindex @code{tgamma()}
-@cindex Gamma function
+@cindex gamma function
 @item @code{lgamma(x)}
-@tab logarithm of Gamma function
+@tab logarithm of gamma function
 @cindex @code{lgamma()}
 @item @code{beta(x, y)}
-@tab Beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
+@tab beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
 @cindex @code{beta()}
 @item @code{psi(x)}
 @tab psi (digamma) function
@@ -4657,9 +5836,9 @@ GiNaC contains the following predefined mathematical functions:
 @item @code{psi(n, x)}
 @tab derivatives of psi function (polygamma functions)
 @item @code{factorial(n)}
-@tab factorial function
+@tab factorial function @math{n!}
 @cindex @code{factorial()}
-@item @code{binomial(n, m)}
+@item @code{binomial(n, k)}
 @tab binomial coefficients
 @cindex @code{binomial()}
 @item @code{Order(x)}
@@ -4684,10 +5863,238 @@ 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
 
-@node Solving Linear Systems of Equations, Input/Output, Built-in Functions, Methods and Functions
+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
+@section Solving linear systems of equations
 @cindex @code{lsolve()}
 
 The function @code{lsolve()} provides a convenient wrapper around some
@@ -4695,12 +6102,13 @@ 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
 @code{relational}) while @code{symbols} is a @code{lst} of
-indeterminates.  (@xref{The Class Hierarchy}, for an exposition of class
+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,
@@ -4709,12 +6117,11 @@ 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;
-    eqns.append(a*x+b*y==3).append(x-y==b);
-    lst vars;
-    vars.append(x).append(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)@}
+     // -> @{x==(3+b^2)/(b+a),y==(3-b*a)/(b+a)@}
 @end example
 
 When the linear equations @code{eqns} are underdetermined, the solution
@@ -4726,7 +6133,7 @@ to @code{lsolve()}: it accepts the same parameters as
 around that method.
 
 
-@node Input/Output, Extending GiNaC, Solving Linear Systems of Equations, Methods and Functions
+@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
@@ -4775,7 +6182,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
@@ -4794,6 +6202,7 @@ you can output to a temporary @code{ostringstream} like this:
     // ...
 @end example
 
+@anchor{csrc printing}
 @cindex @code{csrc}
 @cindex @code{csrc_float}
 @cindex @code{csrc_double}
@@ -4866,7 +6275,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}
@@ -4903,7 +6313,7 @@ static void my_print(const ex & e)
     if (is_a<function>(e))
         cout << ex_to<function>(e).get_name();
     else
-        cout << e.bp->class_name();
+        cout << ex_to<basic>(e).class_name();
     cout << "(";
     size_t n = e.nops();
     if (n)
@@ -4962,34 +6372,155 @@ 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:
+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
-#include <iostream>
-#include <string>
-#include <stdexcept>
-#include <ginac/ginac.h>
-using namespace std;
-using namespace GiNaC;
+    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
 
-int main()
-@{
-    symbol x("x");
-    string s;
+The complete filename (including the suffix @code{.so}) of the object file has
+to be supplied.
 
-    cout << "Enter an expression containing 'x': ";
-    getline(cin, s);
+The function
 
-    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;
-    @}
-@}
+@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)
@@ -5057,7 +6588,8 @@ And the stored expressions can be retrieved by their name:
 
 @example
     // ...
-    lst syms(x, y);
+    lst syms;
+    syms = x, y;
 
     ex ex1 = a2.unarchive_ex(syms, "foo");
     ex ex2 = a2.unarchive_ex(syms, "the second one");
@@ -5165,12 +6697,12 @@ Be warned, however, that the set of properties and their meaning for each
 class may change between GiNaC versions.
 
 
-@node Extending GiNaC, What does not belong into GiNaC, Input/Output, Top
+@node Extending GiNaC, What does not belong into GiNaC, Input/output, Top
 @c    node-name, next, previous, up
 @chapter Extending GiNaC
 
 By reading so far you should have gotten a fairly good understanding of
-GiNaC's design-patterns.  From here on you should start reading the
+GiNaC's design patterns.  From here on you should start reading the
 sources.  All we can do now is issue some recommendations how to tackle
 GiNaC's many loose ends in order to fulfill everybody's dreams.  If you
 develop some useful extension please don't hesitate to contact the GiNaC
@@ -5179,6 +6711,7 @@ authors---they will happily incorporate them into future versions.
 @menu
 * What does not belong into GiNaC::  What to avoid.
 * Symbolic functions::               Implementing symbolic functions.
+* Printing::                         Adding new output formats.
 * Structures::                       Defining new algebraic classes (the easy way).
 * Adding classes::                   Defining new algebraic classes (the hard way).
 @end menu
@@ -5208,7 +6741,7 @@ inefficient.  For this purpose, the underlying foundation classes
 provided by CLN are much better suited.
 
 
-@node Symbolic functions, Structures, What does not belong into GiNaC, Extending GiNaC
+@node Symbolic functions, Printing, What does not belong into GiNaC, Extending GiNaC
 @c    node-name, next, previous, up
 @section Symbolic functions
 
@@ -5249,12 +6782,7 @@ that is not further evaluated:
 @example
 DECLARE_FUNCTION_2P(myfcn)
 
-static ex myfcn_eval(const ex & x, const ex & y)
-@{
-    return myfcn(x, y).hold();
-@}
-
-REGISTER_FUNCTION(myfcn, eval_func(myfcn_eval))
+REGISTER_FUNCTION(myfcn, dummy())
 @end example
 
 Any code that has seen the @code{DECLARE_FUNCTION} line can use @code{myfcn()}
@@ -5264,33 +6792,20 @@ in algebraic expressions:
 @{
     ...
     symbol x("x");
-    ex e = 2*myfcn(42, 3*x+1) - x;
-     // this calls myfcn_eval(42, 3*x+1), and inserts its return value into
-     // the actual expression
+    ex e = 2*myfcn(42, 1+3*x) - x;
     cout << e << endl;
      // prints '2*myfcn(42,1+3*x)-x'
     ...
 @}
 @end example
 
-@cindex @code{hold()}
-@cindex evaluation
-The @code{eval_func()} option specifies the C++ function that implements
-the @code{eval()} method, GiNaC's anonymous evaluator. This function takes
-the same number of arguments as the associated symbolic function (two in this
-case) and returns the (possibly transformed or in some way simplified)
-symbolically evaluated function (@xref{Automatic evaluation}, for a description
-of the automatic evaluation process). If no (further) evaluation is to take
-place, the @code{eval_func()} function must return the original function
-with @code{.hold()}, to avoid a potential infinite recursion. If your
-symbolic functions produce a segmentation fault or stack overflow when
-using them in expressions, you are probably missing a @code{.hold()}
-somewhere.
+The @code{dummy()} option in the @code{REGISTER_FUNCTION} line signifies
+"no options". A function with no options specified merely acts as a kind of
+container for its arguments. It is a pure "dummy" function with no associated
+logic (which is, however, sometimes perfectly sufficient).
 
-There is not much you can do with the @code{myfcn} function. It merely acts
-as a kind of container for its arguments (which is, however, sometimes
-perfectly sufficient). Let's have a look at the implementation of GiNaC's
-cosine function.
+Let's now have a look at the implementation of GiNaC's cosine function for an
+example of how to make an "intelligent" function.
 
 @subsection The cosine function
 
@@ -5304,24 +6819,53 @@ which declares to all programs using GiNaC that there is a function @samp{cos}
 that takes one @code{ex} as an argument. This is all they need to know to use
 this function in expressions.
 
-The implementation of the cosine function is in @file{inifcns_trans.cpp}. The
-@code{eval_func()} function looks something like this (actually, it doesn't
-look like this at all, but it should give you an idea what is going on):
+The implementation of the cosine function is in @file{inifcns_trans.cpp}. Here
+is its @code{REGISTER_FUNCTION} line:
+
+@example
+REGISTER_FUNCTION(cos, eval_func(cos_eval).
+                       evalf_func(cos_evalf).
+                       derivative_func(cos_deriv).
+                       latex_name("\\cos"));
+@end example
+
+There are four options defined for the cosine function. One of them
+(@code{latex_name}) gives the function a proper name for LaTeX output; the
+other three indicate the C++ functions in which the "brains" of the cosine
+function are defined.
+
+@cindex @code{hold()}
+@cindex evaluation
+The @code{eval_func()} option specifies the C++ function that implements
+the @code{eval()} method, GiNaC's anonymous evaluator. This function takes
+the same number of arguments as the associated symbolic function (one in this
+case) and returns the (possibly transformed or in some way simplified)
+symbolically evaluated function (@xref{Automatic evaluation}, for a description
+of the automatic evaluation process). If no (further) evaluation is to take
+place, the @code{eval_func()} function must return the original function
+with @code{.hold()}, to avoid a potential infinite recursion. If your
+symbolic functions produce a segmentation fault or stack overflow when
+using them in expressions, you are probably missing a @code{.hold()}
+somewhere.
+
+The @code{eval_func()} function for the cosine looks something like this
+(actually, it doesn't look like this at all, but it should give you an idea
+what is going on):
 
 @example
 static ex cos_eval(const ex & x)
 @{
-    if (<x is a multiple of 2*Pi>)
+    if ("x is a multiple of 2*Pi")
         return 1;
-    else if (<x is a multiple of Pi>)
+    else if ("x is a multiple of Pi")
         return -1;
-    else if (<x is a multiple of Pi/2>)
+    else if ("x is a multiple of Pi/2")
         return 0;
     // more rules...
 
-    else if (<x has the form 'acos(y)'>)
+    else if ("x has the form 'acos(y)'")
         return y;
-    else if (<x has the form 'asin(y)'>)
+    else if ("x has the form 'asin(y)'")
         return sqrt(1-y^2);
     // more rules...
 
@@ -5330,6 +6874,20 @@ static ex cos_eval(const ex & x)
 @}
 @end example
 
+This function is called every time the cosine is used in a symbolic expression:
+
+@example
+@{
+    ...
+    e = cos(Pi);
+     // this calls cos_eval(Pi), and inserts its return value into
+     // the actual expression
+    cout << e << endl;
+     // prints '-1'
+    ...
+@}
+@end example
+
 In this way, @code{cos(4*Pi)} automatically becomes @math{1},
 @code{cos(asin(a+b))} becomes @code{sqrt(1-(a+b)^2)}, etc. If no reasonable
 symbolic transformation can be done, the unmodified function is returned
@@ -5385,7 +6943,7 @@ static ex tan_series(const ex & x, const relational & rel,
     // Find the actual expansion point
     const ex x_pt = x.subs(rel);
 
-    if (<x_pt is not an odd multiple of Pi/2>)
+    if ("x_pt is not an odd multiple of Pi/2")
         throw do_taylor();  // tell function::series() to do Taylor expansion
 
     // On a pole, expand sin()/cos()
@@ -5396,32 +6954,20 @@ static ex tan_series(const ex & x, const relational & rel,
 The @code{series()} implementation of a function @emph{must} return a
 @code{pseries} object, otherwise your code will crash.
 
-Now that all the ingredients have been set up, the @code{REGISTER_FUNCTION}
-macro is used to tell the system how the @code{cos()} function behaves:
-
-@example
-REGISTER_FUNCTION(cos, eval_func(cos_eval).
-                       evalf_func(cos_evalf).
-                       derivative_func(cos_deriv).
-                       latex_name("\\cos"));
-@end example
-
-This registers the @code{cos_eval()}, @code{cos_evalf()} and
-@code{cos_deriv()} C++ functions with the @code{cos()} function, and also
-gives it a proper LaTeX name.
-
 @subsection Function options
 
 GiNaC functions understand several more options which are always
 specified as @code{.option(params)}. None of them are required, but you
-need to specify at least one option to @code{REGISTER_FUNCTION()} (usually
-the @code{eval()} method).
+need to specify at least one option to @code{REGISTER_FUNCTION()}. There
+is a do-nothing option called @code{dummy()} which you can use to define
+functions without any special options.
 
 @example
 eval_func(<C++ function>)
 evalf_func(<C++ function>)
 derivative_func(<C++ function>)
 series_func(<C++ function>)
+conjugate_func(<C++ function>)
 @end example
 
 These specify the C++ functions that implement symbolic evaluation,
@@ -5473,8 +7019,369 @@ arguments. @xref{Indexed objects}, for an explanation of symmetry
 specifications. GiNaC will automatically rearrange the arguments of
 symmetric functions into a canonical order.
 
+Sometimes you may want to have finer control over how functions are
+displayed in the output. For example, the @code{abs()} function prints
+itself as @samp{abs(x)} in the default output format, but as @samp{|x|}
+in LaTeX mode, and @code{fabs(x)} in C source output. This is achieved
+with the
+
+@example
+print_func<C>(<C++ function>)
+@end example
+
+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
+@section GiNaC's expression output system
+
+GiNaC allows the output of expressions in a variety of different formats
+(@pxref{Input/output}). This section will explain how expression output
+is implemented internally, and how to define your own output formats or
+change the output format of built-in algebraic objects. You will also want
+to read this section if you plan to write your own algebraic classes or
+functions.
+
+@cindex @code{print_context} (class)
+@cindex @code{print_dflt} (class)
+@cindex @code{print_latex} (class)
+@cindex @code{print_tree} (class)
+@cindex @code{print_csrc} (class)
+All the different output formats are represented by a hierarchy of classes
+rooted in the @code{print_context} class, defined in the @file{print.h}
+header file:
+
+@table @code
+@item print_dflt
+the default output format
+@item print_latex
+output in LaTeX mathematical mode
+@item print_tree
+a dump of the internal expression structure (for debugging)
+@item print_csrc
+the base class for C source output
+@item print_csrc_float
+C source output using the @code{float} type
+@item print_csrc_double
+C source output using the @code{double} type
+@item print_csrc_cl_N
+C source output using CLN types
+@end table
+
+The @code{print_context} base class provides two public data members:
+
+@example
+class print_context
+@{
+    ...
+public:
+    std::ostream & s;
+    unsigned options;
+@};
+@end example
+
+@code{s} is a reference to the stream to output to, while @code{options}
+holds flags and modifiers. Currently, there is only one flag defined:
+@code{print_options::print_index_dimensions} instructs the @code{idx} class
+to print the index dimension which is normally hidden.
+
+When you write something like @code{std::cout << e}, where @code{e} is
+an object of class @code{ex}, GiNaC will construct an appropriate
+@code{print_context} object (of a class depending on the selected output
+format), fill in the @code{s} and @code{options} members, and call
+
+@cindex @code{print()}
+@example
+void ex::print(const print_context & c, unsigned level = 0) const;
+@end example
+
+which in turn forwards the call to the @code{print()} method of the
+top-level algebraic object contained in the expression.
+
+Unlike other methods, GiNaC classes don't usually override their
+@code{print()} method to implement expression output. Instead, the default
+implementation @code{basic::print(c, level)} performs a run-time double
+dispatch to a function selected by the dynamic type of the object and the
+passed @code{print_context}. To this end, GiNaC maintains a separate method
+table for each class, similar to the virtual function table used for ordinary
+(single) virtual function dispatch.
+
+The method table contains one slot for each possible @code{print_context}
+type, indexed by the (internally assigned) serial number of the type. Slots
+may be empty, in which case GiNaC will retry the method lookup with the
+@code{print_context} object's parent class, possibly repeating the process
+until it reaches the @code{print_context} base class. If there's still no
+method defined, the method table of the algebraic object's parent class
+is consulted, and so on, until a matching method is found (eventually it
+will reach the combination @code{basic/print_context}, which prints the
+object's class name enclosed in square brackets).
+
+You can think of the print methods of all the different classes and output
+formats as being arranged in a two-dimensional matrix with one axis listing
+the algebraic classes and the other axis listing the @code{print_context}
+classes.
+
+Subclasses of @code{basic} can, of course, also overload @code{basic::print()}
+to implement printing, but then they won't get any of the benefits of the
+double dispatch mechanism (such as the ability for derived classes to
+inherit only certain print methods from its parent, or the replacement of
+methods at run-time).
+
+@subsection Print methods for classes
+
+The method table for a class is set up either in the definition of the class,
+by passing the appropriate @code{print_func<C>()} option to
+@code{GINAC_IMPLEMENT_REGISTERED_CLASS_OPT()} (@xref{Adding classes}, for
+an example), or at run-time using @code{set_print_func<T, C>()}. The latter
+can also be used to override existing methods dynamically.
+
+The argument to @code{print_func<C>()} and @code{set_print_func<T, C>()} can
+be a member function of the class (or one of its parent classes), a static
+member function, or an ordinary (global) C++ function. The @code{C} template
+parameter specifies the appropriate @code{print_context} type for which the
+method should be invoked, while, in the case of @code{set_print_func<>()}, the
+@code{T} parameter specifies the algebraic class (for @code{print_func<>()},
+the class is the one being implemented by
+@code{GINAC_IMPLEMENT_REGISTERED_CLASS_OPT}).
+
+For print methods that are member functions, their first argument must be of
+a type convertible to a @code{const C &}, and the second argument must be an
+@code{unsigned}.
+
+For static members and global functions, the first argument must be of a type
+convertible to a @code{const T &}, the second argument must be of a type
+convertible to a @code{const C &}, and the third argument must be an
+@code{unsigned}. A global function will, of course, not have access to
+private and protected members of @code{T}.
+
+The @code{unsigned} argument of the print methods (and of @code{ex::print()}
+and @code{basic::print()}) is used for proper parenthesizing of the output
+(and by @code{print_tree} for proper indentation). It can be used for similar
+purposes if you write your own output formats.
+
+The explanations given above may seem complicated, but in practice it's
+really simple, as shown in the following example. Suppose that we want to
+display exponents in LaTeX output not as superscripts but with little
+upwards-pointing arrows. This can be achieved in the following way:
+
+@example
+void my_print_power_as_latex(const power & p,
+                             const print_latex & c,
+                             unsigned level)
+@{
+    // get the precedence of the 'power' class
+    unsigned power_prec = p.precedence();
+
+    // if the parent operator has the same or a higher precedence
+    // we need parentheses around the power
+    if (level >= power_prec)
+        c.s << '(';
+
+    // print the basis and exponent, each enclosed in braces, and
+    // separated by an uparrow
+    c.s << '@{';
+    p.op(0).print(c, power_prec);
+    c.s << "@}\\uparrow@{";
+    p.op(1).print(c, power_prec);
+    c.s << '@}';
+
+    // don't forget the closing parenthesis
+    if (level >= power_prec)
+        c.s << ')';
+@}
+                                                                                
+int main()
+@{
+    // a sample expression
+    symbol x("x"), y("y");
+    ex e = -3*pow(x, 3)*pow(y, -2) + pow(x+y, 2) - 1;
+
+    // switch to LaTeX mode
+    cout << latex;
+
+    // this prints "-1+@{(y+x)@}^@{2@}-3 \frac@{x^@{3@}@}@{y^@{2@}@}"
+    cout << e << endl;
+
+    // now we replace the method for the LaTeX output of powers with
+    // 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@}@}"
+    cout << e << endl;
+@}
+@end example
+
+Some notes:
+
+@itemize
+
+@item
+The first argument of @code{my_print_power_as_latex} could also have been
+a @code{const basic &}, the second one a @code{const print_context &}.
+
+@item
+The above code depends on @code{mul} objects converting their operands to
+@code{power} objects for the purpose of printing.
+
+@item
+The output of products including negative powers as fractions is also
+controlled by the @code{mul} class.
+
+@item
+The @code{power/print_latex} method provided by GiNaC prints square roots
+using @code{\sqrt}, but the above code doesn't.
+
+@end itemize
+
+It's not possible to restore a method table entry to its previous or default
+value. Once you have called @code{set_print_func()}, you can only override
+it with another call to @code{set_print_func()}, but you can't easily go back
+to the default behavior again (you can, of course, dig around in the GiNaC
+sources, find the method that is installed at startup
+(@code{power::do_print_latex} in this case), and @code{set_print_func} that
+one; that is, after you circumvent the C++ member access control@dots{}).
+
+@subsection Print methods for functions
+
+Symbolic functions employ a print method dispatch mechanism similar to the
+one used for classes. The methods are specified with @code{print_func<C>()}
+function options. If you don't specify any special print methods, the function
+will be printed with its name (or LaTeX name, if supplied), followed by a
+comma-separated list of arguments enclosed in parentheses.
+
+For example, this is what GiNaC's @samp{abs()} function is defined like:
+
+@example
+static ex abs_eval(const ex & arg) @{ ... @}
+static ex abs_evalf(const ex & arg) @{ ... @}
+                                                                                
+static void abs_print_latex(const ex & arg, const print_context & c)
+@{
+    c.s << "@{|"; arg.print(c); c.s << "|@}";
+@}
+                                                                                
+static void abs_print_csrc_float(const ex & arg, const print_context & c)
+@{
+    c.s << "fabs("; arg.print(c); c.s << ")";
+@}
+                                                                                
+REGISTER_FUNCTION(abs, eval_func(abs_eval).
+                       evalf_func(abs_evalf).
+                       print_func<print_latex>(abs_print_latex).
+                       print_func<print_csrc_float>(abs_print_csrc_float).
+                       print_func<print_csrc_double>(abs_print_csrc_float));
+@end example
+
+This will display @samp{abs(x)} as @samp{|x|} in LaTeX mode and @code{fabs(x)}
+in non-CLN C source output, but as @code{abs(x)} in all other formats.
+
+There is currently no equivalent of @code{set_print_func()} for functions.
+
+@subsection Adding new output formats
+
+Creating a new output format involves subclassing @code{print_context},
+which is somewhat similar to adding a new algebraic class
+(@pxref{Adding classes}). There is a macro @code{GINAC_DECLARE_PRINT_CONTEXT}
+that needs to go into the class definition, and a corresponding macro
+@code{GINAC_IMPLEMENT_PRINT_CONTEXT} that has to appear at global scope.
+Every @code{print_context} class needs to provide a default constructor
+and a constructor from an @code{std::ostream} and an @code{unsigned}
+options value.
+
+Here is an example for a user-defined @code{print_context} class:
+
+@example
+class print_myformat : public print_dflt
+@{
+    GINAC_DECLARE_PRINT_CONTEXT(print_myformat, print_dflt)
+public:
+    print_myformat(std::ostream & os, unsigned opt = 0)
+     : print_dflt(os, opt) @{@}
+@};
+
+print_myformat::print_myformat() : print_dflt(std::cout) @{@}
+
+GINAC_IMPLEMENT_PRINT_CONTEXT(print_myformat, print_dflt)
+@end example
+
+That's all there is to it. None of the actual expression output logic is
+implemented in this class. It merely serves as a selector for choosing
+a particular format. The algorithms for printing expressions in the new
+format are implemented as print methods, as described above.
+
+@code{print_myformat} is a subclass of @code{print_dflt}, so it behaves
+exactly like GiNaC's default output format:
+
+@example
+@{
+    symbol x("x");
+    ex e = pow(x, 2) + 1;
+
+    // this prints "1+x^2"
+    cout << e << endl;
+    
+    // this also prints "1+x^2"
+    e.print(print_myformat()); cout << endl;
+
+    ...
+@}
+@end example
+
+To fill @code{print_myformat} with life, we need to supply appropriate
+print methods with @code{set_print_func()}, like this:
+
+@example
+// This prints powers with '**' instead of '^'. See the LaTeX output
+// example above for explanations.
+void print_power_as_myformat(const power & p,
+                             const print_myformat & c,
+                             unsigned level)
+@{
+    unsigned power_prec = p.precedence();
+    if (level >= power_prec)
+        c.s << '(';
+    p.op(0).print(c, power_prec);
+    c.s << "**";
+    p.op(1).print(c, power_prec);
+    if (level >= power_prec)
+        c.s << ')';
+@}
+
+@{
+    ...
+    // install a new print method for power objects
+    set_print_func<power, print_myformat>(print_power_as_myformat);
+
+    // now this prints "1+x**2"
+    e.print(print_myformat()); cout << endl;
+
+    // but the default format is still "1+x^2"
+    cout << e << endl;
+@}
+@end example
+
 
-@node Structures, Adding classes, Symbolic functions, Extending GiNaC
+@node Structures, Adding classes, Printing, Extending GiNaC
 @c    node-name, next, previous, up
 @section Structures
 
@@ -5600,10 +7507,15 @@ desired, most notably proper output:
 @end example
 
 By default, any structure types you define will be printed as
-@samp{[structure object]}. To override this, you can specialize the
-template's @code{print()} member function. The member functions of
-GiNaC classes are described in more detail in the next section, but
-it shouldn't be hard to figure out what's going on here:
+@samp{[structure object]}. To override this you can either specialize the
+template's @code{print()} member function, or specify print methods with
+@code{set_print_func<>()}, as described in @ref{Printing}. Unfortunately,
+it's not possible to supply class options like @code{print_func<>()} to
+structures, so for a self-contained structure type you need to resort to
+overriding the @code{print()} function, which is also what we will do here.
+
+The member functions of GiNaC classes are described in more detail in the
+next section, but it shouldn't be hard to figure out what's going on here:
 
 @example
 void sprod::print(const print_context & c, unsigned level) const
@@ -5658,7 +7570,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
 
@@ -5824,7 +7737,7 @@ Note that the unarchiving constructor is @code{sprod::structure} and not
 @code{sprod::unarchive()} function.
 
 
-@node Adding classes, A Comparison With Other CAS, Structures, Extending GiNaC
+@node Adding classes, A comparison with other CAS, Structures, Extending GiNaC
 @c    node-name, next, previous, up
 @section Adding classes
 
@@ -5866,9 +7779,8 @@ The RTTI in GiNaC is based on two mechanisms:
 
 @item
 The @code{basic} class declares a member variable @code{tinfo_key} which
-holds an unsigned integer that identifies the object's class. These numbers
-are defined in the @file{tinfos.h} header file for the built-in GiNaC
-classes. They all start with @code{TINFO_}.
+holds a variable of @code{tinfo_t} type (which is actually just
+@code{const void*}) that identifies the object's class.
 
 @item
 By means of some clever tricks with static members, GiNaC maintains a list
@@ -5881,7 +7793,7 @@ to the unarchiving functions. This class registry is defined in the
 
 The disadvantage of this proprietary RTTI implementation is that there's
 a little more to do when implementing new classes (C++'s RTTI works more
-or less automatic) but don't worry, most of the work is simplified by
+or less automatically) but don't worry, most of the work is simplified by
 macros.
 
 @subsection A minimalistic example
@@ -5905,17 +7817,6 @@ using namespace std;
 using namespace GiNaC;
 @end example
 
-The first thing we have to do is to define a @code{tinfo_key} for our new
-class. This can be any arbitrary unsigned number that is not already taken
-by one of the existing classes but it's better to come up with something
-that is unlikely to clash with keys that might be added in the future. The
-numbers in @file{tinfos.h} are modeled somewhat after the class hierarchy
-which is not a requirement but we are going to stick with this scheme:
-
-@example
-const unsigned TINFO_mystring = 0x42420001U;
-@end example
-
 Now we can write down the class declaration. The class stores a C++
 @code{string} and the user shall be able to construct a @code{mystring}
 object from a C or C++ string:
@@ -5926,8 +7827,8 @@ class mystring : public basic
     GINAC_DECLARE_REGISTERED_CLASS(mystring, basic)
   
 public:
-    mystring(const string &s);
-    mystring(const char *s);
+    mystring(const string & s);
+    mystring(const char * s);
 
 private:
     string str;
@@ -5937,7 +7838,7 @@ GINAC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)
 @end example
 
 The @code{GINAC_DECLARE_REGISTERED_CLASS} and @code{GINAC_IMPLEMENT_REGISTERED_CLASS}
-macros are defined in @file{registrar.h}.  They take the name of the class
+macros are defined in @file{registrar.h}. They take the name of the class
 and its direct superclass as arguments and insert all required declarations
 for the RTTI system. The @code{GINAC_DECLARE_REGISTERED_CLASS} should be
 the first line after the opening brace of the class definition. The
@@ -5946,10 +7847,12 @@ source (at global scope, of course, not inside a function).
 
 @code{GINAC_DECLARE_REGISTERED_CLASS} contains, among other things the
 declarations of the default constructor and a couple of other functions that
-are required.  It also defines a type @code{inherited} which refers to the
+are required. It also defines a type @code{inherited} which refers to the
 superclass so you don't have to modify your code every time you shuffle around
-the class hierarchy.  @code{GINAC_IMPLEMENT_REGISTERED_CLASS} registers the
-class with the GiNaC RTTI.
+the class hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} registers the
+class with the GiNaC RTTI (there is also a
+@code{GINAC_IMPLEMENT_REGISTERED_CLASS_OPT} which allows specifying additional
+options for the class, and which we will be using instead in a few minutes).
 
 Now there are seven member functions we have to implement to get a working
 class:
@@ -5960,27 +7863,27 @@ class:
 @code{mystring()}, the default constructor.
 
 @item
-@code{void archive(archive_node &n)}, the archiving function. This stores all
+@code{void archive(archive_node & n)}, the archiving function. This stores all
 information needed to reconstruct an object of this class inside an
 @code{archive_node}.
 
 @item
-@code{mystring(const archive_node &n, lst &sym_lst)}, the unarchiving
+@code{mystring(const archive_node & n, lst & sym_lst)}, the unarchiving
 constructor. This constructs an instance of the class from the information
 found in an @code{archive_node}.
 
 @item
-@code{ex unarchive(const archive_node &n, lst &sym_lst)}, the static
+@code{ex unarchive(const archive_node & n, lst & sym_lst)}, the static
 unarchiving function. It constructs a new instance by calling the unarchiving
 constructor.
 
 @item
 @cindex @code{compare_same_type()}
-@code{int compare_same_type(const basic &other)}, which is used internally
+@code{int compare_same_type(const basic & other)}, which is used internally
 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
@@ -5989,7 +7892,7 @@ objects for which no reasonable algebraic ordering relationship can be
 defined.
 
 @item
-And, of course, @code{mystring(const string &s)} and @code{mystring(const char *s)}
+And, of course, @code{mystring(const string & s)} and @code{mystring(const char * s)}
 which are the two constructors we declared.
 
 @end itemize
@@ -5997,11 +7900,14 @@ which are the two constructors we declared.
 Let's proceed step-by-step. The default constructor looks like this:
 
 @example
-mystring::mystring() : inherited(TINFO_mystring) @{@}
+mystring::mystring() : inherited(&mystring::tinfo_static) @{@}
 @end example
 
 The golden rule is that in all constructors you have to set the
-@code{tinfo_key} member to the @code{TINFO_*} value of your class. Otherwise
+@code{tinfo_key} member to the @code{&your_class_name::tinfo_static}
+@footnote{Each GiNaC class has a static member called tinfo_static.
+This member is declared by the GINAC_DECLARE_REGISTERED_CLASS macros
+and defined by the GINAC_IMPLEMENT_REGISTERED_CLASS macros.}. Otherwise
 it will be set by the constructor of the superclass and all hell will break
 loose in the RTTI. For your convenience, the @code{basic} class provides
 a constructor that takes a @code{tinfo_key} value, which we are using here
@@ -6018,7 +7924,7 @@ if you don't plan to use archives, but the minimum required implementation
 is really simple.  First, the archiving function:
 
 @example
-void mystring::archive(archive_node &n) const
+void mystring::archive(archive_node & n) const
 @{
     inherited::archive(n);
     n.add_string("string", str);
@@ -6036,7 +7942,7 @@ The unarchiving constructor is basically the inverse of the archiving
 function:
 
 @example
-mystring::mystring(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
+mystring::mystring(const archive_node & n, lst & sym_lst) : inherited(n, sym_lst)
 @{
     n.find_string("string", str);
 @}
@@ -6050,7 +7956,7 @@ by the unarchiving constructor of the @code{basic} class.
 Finally, the unarchiving function:
 
 @example
-ex mystring::unarchive(const archive_node &n, lst &sym_lst)
+ex mystring::unarchive(const archive_node & n, lst & sym_lst)
 @{
     return (new mystring(n, sym_lst))->setflag(status_flags::dynallocated);
 @}
@@ -6069,7 +7975,7 @@ Our @code{compare_same_type()} function uses a provided function to compare
 the string members:
 
 @example
-int mystring::compare_same_type(const basic &other) const
+int mystring::compare_same_type(const basic & other) const
 @{
     const mystring &o = static_cast<const mystring &>(other);
     int cmpval = str.compare(o.str);
@@ -6091,8 +7997,10 @@ all relevant member variables.
 Now the only thing missing is our two new constructors:
 
 @example
-mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s) @{@}
-mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s) @{@}
+mystring::mystring(const string & s)
+    : inherited(&mystring::tinfo_static), str(s) @{@}
+mystring::mystring(const char * s)
+    : inherited(&mystring::tinfo_static), str(s) @{@}
 @end example
 
 No surprises here. We set the @code{str} member from the argument and
@@ -6106,7 +8014,7 @@ ex e = mystring("Hello, world!");
 cout << is_a<mystring>(e) << endl;
  // -> 1 (true)
 
-cout << e.bp->class_name() << endl;
+cout << ex_to<basic>(e).class_name() << endl;
  // -> mystring
 @end example
 
@@ -6118,20 +8026,21 @@ cout << e << endl;
 @end example
 
 Hm, not exactly what we expect, but of course the @code{mystring} class
-doesn't yet know how to print itself. This is done in the @code{print()}
-member function. Let's say that we wanted to print the string surrounded
-by double quotes:
+doesn't yet know how to print itself. This can be done either by implementing
+the @code{print()} member function, or, preferably, by specifying a
+@code{print_func<>()} class option. Let's say that we want to print the string
+surrounded by double quotes:
 
 @example
 class mystring : public basic
 @{
     ...
-public:
-    void print(const print_context &c, unsigned level = 0) const;
+protected:
+    void do_print(const print_context & c, unsigned level = 0) const;
     ...
 @};
 
-void mystring::print(const print_context &c, unsigned level) const
+void mystring::do_print(const print_context & c, unsigned level) const
 @{
     // print_context::s is a reference to an ostream
     c.s << '\"' << str << '\"';
@@ -6139,14 +8048,39 @@ void mystring::print(const print_context &c, unsigned level) const
 @end example
 
 The @code{level} argument is only required for container classes to
-correctly parenthesize the output. Let's try again to print the expression:
+correctly parenthesize the output.
+
+Now we need to tell GiNaC that @code{mystring} objects should use the
+@code{do_print()} member function for printing themselves. For this, we
+replace the line
+
+@example
+GINAC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)
+@end example
+
+with
+
+@example
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mystring, basic,
+  print_func<print_context>(&mystring::do_print))
+@end example
+
+Let's try again to print the expression:
 
 @example
 cout << e << endl;
  // -> "Hello, world!"
 @end example
 
-Much better. The @code{mystring} class can be used in arbitrary expressions:
+Much better. If we wanted to have @code{mystring} objects displayed in a
+different way depending on the output format (default, LaTeX, etc.), we
+would have supplied multiple @code{print_func<>()} options with different
+template parameters (@code{print_dflt}, @code{print_latex}, etc.),
+separated by dots. This is similar to the way options are specified for
+symbolic functions. @xref{Printing}, for a more in-depth description of the
+way expression output is implemented in GiNaC.
+
+The @code{mystring} class can be used in arbitrary expressions:
 
 @example
 e += mystring("GiNaC rulez"); 
@@ -6205,7 +8139,7 @@ public:
 ex mystring::eval(int level) const
 @{
     string new_str;
-    for (int i=0; i<str.length(); i++) @{
+    for (size_t i=0; i<str.length(); i++) @{
         char c = str[i];
         if (c >= 'A' && c <= 'Z') 
             new_str += tolower(c);
@@ -6251,7 +8185,7 @@ required but will make operations with objects of the class more efficient:
 @cindex @code{is_equal_same_type()}
 @example
 unsigned calchash() const;
-bool is_equal_same_type(const basic &other) const;
+bool is_equal_same_type(const basic & other) const;
 @end example
 
 The @code{calchash()} method returns an @code{unsigned} hash value for the
@@ -6303,8 +8237,30 @@ should become a need.
 
 That's it. May the source be with you!
 
+@subsection Upgrading extension classes from older version of GiNaC
+
+If you got some extension classes for GiNaC 1.3.X some changes are
+necessary in order to make your code work with GiNaC 1.4.
+
+@itemize @bullet
+@item constructors which set @code{tinfo_key} such as
+
+@example
+myclass::myclass() : inherited(TINFO_myclass) @{@}
+@end example
+
+need to be rewritten as
+
+@example
+myclass::myclass() : inherited(&myclass::tinfo_static) @{@}
+@end example
+
+@item TINO_myclass is not necessary any more and can be removed.
+
+@end itemize
+
 
-@node A Comparison With Other CAS, Advantages, Adding classes, Top
+@node A comparison with other CAS, Advantages, Adding classes, Top
 @c    node-name, next, previous, up
 @chapter A Comparison With Other CAS
 @cindex advocacy
@@ -6320,7 +8276,7 @@ disadvantages over these systems.
 * Why C++?::                         Attractiveness of C++.
 @end menu
 
-@node Advantages, Disadvantages, A Comparison With Other CAS, A Comparison With Other CAS
+@node Advantages, Disadvantages, A comparison with other CAS, A comparison with other CAS
 @c    node-name, next, previous, up
 @section Advantages
 
@@ -6400,7 +8356,7 @@ speed with other CAS.
 @end itemize
 
 
-@node Disadvantages, Why C++?, Advantages, A Comparison With Other CAS
+@node Disadvantages, Why C++?, Advantages, A comparison with other CAS
 @c    node-name, next, previous, up
 @section Disadvantages
 
@@ -6435,7 +8391,7 @@ yet ANSI compliant, support all needed features.
 @end itemize
 
 
-@node Why C++?, Internal Structures, Disadvantages, A Comparison With Other CAS
+@node Why C++?, Internal structures, Disadvantages, A comparison with other CAS
 @c    node-name, next, previous, up
 @section Why C++?
 
@@ -6452,16 +8408,16 @@ Furthermore, the main developers are more familiar with C++ than with
 any other programming language.
 
 
-@node Internal Structures, Expressions are reference counted, Why C++? , Top
+@node Internal structures, Expressions are reference counted, Why C++? , Top
 @c    node-name, next, previous, up
-@appendix Internal Structures
+@appendix Internal structures
 
 @menu
 * Expressions are reference counted::
 * Internal representation of products and sums::
 @end menu
 
-@node Expressions are reference counted, Internal representation of products and sums, Internal Structures, Internal Structures
+@node Expressions are reference counted, Internal representation of products and sums, Internal structures, Internal structures
 @c    node-name, next, previous, up
 @appendixsection Expressions are reference counted
 
@@ -6545,12 +8501,12 @@ inserted.  But it may be useful to remember that this is not what
 happens.  Knowing this will enable you to write much more efficient
 code.  If you still have an uncertain feeling with copy-on-write
 semantics, we recommend you have a look at the
-@uref{http://www.cerfnet.com/~mpcline/c++-faq-lite/, C++-FAQ lite} by
+@uref{http://www.parashift.com/c++-faq-lite/, C++-FAQ lite} by
 Marshall Cline.  Chapter 16 covers this issue and presents an
 implementation which is pretty close to the one in GiNaC.
 
 
-@node Internal representation of products and sums, Package Tools, Expressions are reference counted, Internal Structures
+@node Internal representation of products and sums, Package tools, Expressions are reference counted, Internal structures
 @c    node-name, next, previous, up
 @appendixsection Internal representation of products and sums
 
@@ -6618,77 +8574,38 @@ expansion and the like are reimplemented for @code{add} and @code{mul},
 but the data structure is inherited from @code{expairseq}.
 
 
-@node Package Tools, ginac-config, Internal representation of products and sums, Top
+@node Package tools, Configure script options, Internal representation of products and sums, Top
 @c    node-name, next, previous, up
-@appendix Package Tools
+@appendix Package tools
 
 If you are creating a software package that uses the GiNaC library,
-setting the correct command line options for the compiler and linker
-can be difficult. GiNaC includes two tools to make this process easier.
-
-@menu
-* ginac-config::   A shell script to detect compiler and linker flags.
-* AM_PATH_GINAC::  Macro for GNU automake.
-@end menu
-
-
-@node ginac-config, AM_PATH_GINAC, Package Tools, Package Tools
-@c    node-name, next, previous, up
-@section @command{ginac-config}
-@cindex ginac-config
-
-@command{ginac-config} is a shell script that you can use to determine
-the compiler and linker command line options required to compile and
-link a program with the GiNaC library.
-
-@command{ginac-config} takes the following flags:
-
-@table @samp
-@item --version
-Prints out the version of GiNaC installed.
-@item --cppflags
-Prints '-I' flags pointing to the installed header files.
-@item --libs
-Prints out the linker flags necessary to link a program against GiNaC.
-@item --prefix[=@var{PREFIX}]
-If @var{PREFIX} is specified, overrides the configured value of @env{$prefix}.
-(And of exec-prefix, unless @code{--exec-prefix} is also specified)
-Otherwise, prints out the configured value of @env{$prefix}.
-@item --exec-prefix[=@var{PREFIX}]
-If @var{PREFIX} is specified, overrides the configured value of @env{$exec_prefix}.
-Otherwise, prints out the configured value of @env{$exec_prefix}.
-@end table
-
-Typically, @command{ginac-config} will be used within a configure
-script, as described below. It, however, can also be used directly from
-the command line using backquotes to compile a simple program. For
-example:
-
+setting the correct command line options for the compiler and linker can
+be difficult.  The @command{pkg-config} utility makes this process
+easier.  GiNaC supplies all necessary data in @file{ginac.pc} (installed
+into @code{/usr/local/lib/pkgconfig} by default). To compile a simple
+program use @footnote{If GiNaC is installed into some non-standard
+directory @var{prefix} one should set the @var{PKG_CONFIG_PATH}
+environment variable to @var{prefix}/lib/pkgconfig for this to work.}
 @example
-c++ -o simple `ginac-config --cppflags` simple.cpp `ginac-config --libs`
+g++ -o simple `pkg-config --cflags --libs ginac` simple.cpp
 @end example
 
 This command line might expand to (for example):
-
 @example
-cc -o simple -I/usr/local/include simple.cpp -L/usr/local/lib \
-  -lginac -lcln -lstdc++
+g++ -o simple -lginac -lcln simple.cpp
 @end example
 
-Not only is the form using @command{ginac-config} easier to type, it will
+Not only is the form using @command{pkg-config} easier to type, it will
 work on any system, no matter how GiNaC was configured.
 
-
-@node AM_PATH_GINAC, Configure script options, ginac-config, Package Tools
-@c    node-name, next, previous, up
-@section @samp{AM_PATH_GINAC}
-@cindex AM_PATH_GINAC
-
-For packages configured using GNU automake, GiNaC also provides
-a macro to automate the process of checking for GiNaC.
+For packages configured using GNU automake, @command{pkg-config} also
+provides the @code{PKG_CHECK_MODULES} macro to automate the process of
+checking for libraries
 
 @example
-AM_PATH_GINAC([@var{MINIMUM-VERSION}, [@var{ACTION-IF-FOUND} [, @var{ACTION-IF-NOT-FOUND}]]])
+PKG_CHECK_MODULES(MYAPP, ginac >= MINIMUM_VERSION, 
+                  [@var{ACTION-IF-FOUND}],
+                  [@var{ACTION-IF-NOT-FOUND}])
 @end example
 
 This macro:
@@ -6696,115 +8613,94 @@ This macro:
 @itemize @bullet
 
 @item
-Determines the location of GiNaC using @command{ginac-config}, which is
-either found in the user's path, or from the environment variable
-@env{GINACLIB_CONFIG}.
+Determines the location of GiNaC using data from @file{ginac.pc}, which is
+either found in the default @command{pkg-config} search path, or from 
+the environment variable @env{PKG_CONFIG_PATH}.
 
 @item
 Tests the installed libraries to make sure that their version
-is later than @var{MINIMUM-VERSION}. (A default version will be used
-if not specified)
+is later than @var{MINIMUM-VERSION}.
 
 @item
-If the required version was found, sets the @env{GINACLIB_CPPFLAGS} variable
-to the output of @command{ginac-config --cppflags} and the @env{GINACLIB_LIBS}
-variable to the output of @command{ginac-config --libs}, and calls
+If the required version was found, sets the @env{MYAPP_CFLAGS} variable
+to the output of @command{pkg-config --cflags ginac} and the @env{MYAPP_LIBS}
+variable to the output of @command{pkg-config --libs ginac}, and calls
 @samp{AC_SUBST()} for these variables so they can be used in generated
 makefiles, and then executes @var{ACTION-IF-FOUND}.
 
 @item
-If the required version was not found, sets @env{GINACLIB_CPPFLAGS} and
-@env{GINACLIB_LIBS} to empty strings, and executes @var{ACTION-IF-NOT-FOUND}.
+If the required version was not found, executes @var{ACTION-IF-NOT-FOUND}.
 
 @end itemize
 
-This macro is in file @file{ginac.m4} which is installed in
-@file{$datadir/aclocal}. Note that if automake was installed with a
-different @samp{--prefix} than GiNaC, you will either have to manually
-move @file{ginac.m4} to automake's @file{$datadir/aclocal}, or give
-aclocal the @samp{-I} option when running it.
-
 @menu
-* Configure script options::  Configuring a package that uses AM_PATH_GINAC.
-* Example package::           Example of a package using AM_PATH_GINAC.
+* Configure script options::  Configuring a package that uses GiNaC
+* Example package::           Example of a package using GiNaC
 @end menu
 
 
-@node Configure script options, Example package, AM_PATH_GINAC, AM_PATH_GINAC
+@node Configure script options, Example package, Package tools, Package tools 
 @c    node-name, next, previous, up
-@subsection Configuring a package that uses @samp{AM_PATH_GINAC}
-
-Simply make sure that @command{ginac-config} is in your path, and run
-the configure script.
-
-Notes:
+@subsection Configuring a package that uses GiNaC
 
-@itemize @bullet
-
-@item
 The directory where the GiNaC libraries are installed needs
-to be found by your system's dynamic linker.
-  
-This is generally done by
+to be found by your system's dynamic linkers (both compile- and run-time
+ones).  See the documentation of your system linker for details.  Also
+make sure that @file{ginac.pc} is in @command{pkg-config}'s search path,
+@xref{pkg-config, ,pkg-config, *manpages*}.
 
-@display
-editing @file{/etc/ld.so.conf} and running @command{ldconfig}
-@end display
+The short summary below describes how to do this on a GNU/Linux
+system.
 
-or by
-   
-@display
-setting the environment variable @env{LD_LIBRARY_PATH},
-@end display
-
-or, as a last resort, 
-@display
-giving a @samp{-R} or @samp{-rpath} flag (depending on your linker) when
-running configure, for instance:
+Suppose GiNaC is installed into the directory @samp{PREFIX}. To tell
+the linkers where to find the library one should
 
+@itemize @bullet
+@item
+edit @file{/etc/ld.so.conf} and run @command{ldconfig}. For example,
 @example
-LDFLAGS=-R/home/cbauer/lib ./configure
+# echo PREFIX/lib >> /etc/ld.so.conf
+# ldconfig
 @end example
-@end display
 
 @item
-You can also specify a @command{ginac-config} not in your path by
-setting the @env{GINACLIB_CONFIG} environment variable to the
-name of the executable
+or set the environment variables @env{LD_LIBRARY_PATH} and @env{LD_RUN_PATH}
+@example
+$ export LD_LIBRARY_PATH=PREFIX/lib
+$ export LD_RUN_PATH=PREFIX/lib
+@end example
 
 @item
-If you move the GiNaC package from its installed location,
-you will either need to modify @command{ginac-config} script
-manually to point to the new location or rebuild GiNaC.
+or give a @samp{-L} and @samp{--rpath} flags when running configure,
+for instance:
 
+@example
+$ LDFLAGS='-Wl,-LPREFIX/lib -Wl,--rpath=PREFIX/lib' ./configure
+@end example
 @end itemize
 
-Advanced note:
-
-@itemize @bullet
-@item
-configure flags
-  
+To tell @command{pkg-config} where the @file{ginac.pc} file is,
+set the @env{PKG_CONFIG_PATH} environment variable:
 @example
---with-ginac-prefix=@var{PREFIX}
---with-ginac-exec-prefix=@var{PREFIX}
+$ export PKG_CONFIG_PATH=PREFIX/lib/pkgconfig
 @end example
 
-are provided to override the prefix and exec-prefix that were stored
-in the @command{ginac-config} shell script by GiNaC's configure. You are
-generally better off configuring GiNaC with the right path to begin with.
-@end itemize
+Finally, run the @command{configure} script
+@example
+$ ./configure 
+@end example
 
+@c There are many other ways to do the same, @xref{Options, ,Command Line Options, ld, GNU ld manual}.
 
-@node Example package, Bibliography, Configure script options, AM_PATH_GINAC
+@node Example package, Bibliography, Configure script options, Package tools
 @c    node-name, next, previous, up
-@subsection Example of a package using @samp{AM_PATH_GINAC}
+@subsection Example of a package using GiNaC
 
 The following shows how to build a simple package using automake
-and the @samp{AM_PATH_GINAC} macro. The program used here is @file{simple.cpp}:
+and the @samp{PKG_CHECK_MODULES} macro. The program used here is @file{simple.cpp}:
 
 @example
+#include <iostream>
 #include <ginac/ginac.h>
 
 int main()
@@ -6820,33 +8716,39 @@ int main()
 You should first read the introductory portions of the automake
 Manual, if you are not already familiar with it.
 
-Two files are needed, @file{configure.in}, which is used to build the
+Two files are needed, @file{configure.ac}, which is used to build the
 configure script:
 
 @example
-dnl Process this file with autoconf to produce a configure script.
-AC_INIT(simple.cpp)
-AM_INIT_AUTOMAKE(simple.cpp, 1.0.0)
+dnl Process this file with autoreconf to produce a configure script.
+AC_INIT([simple], 1.0.0, bogus@@example.net)
+AC_CONFIG_SRCDIR(simple.cpp)
+AM_INIT_AUTOMAKE([foreign 1.8])
 
 AC_PROG_CXX
 AC_PROG_INSTALL
-AC_LANG_CPLUSPLUS
+AC_LANG([C++])
 
-AM_PATH_GINAC(0.9.0, [
-  LIBS="$LIBS $GINACLIB_LIBS"
-  CPPFLAGS="$CPPFLAGS $GINACLIB_CPPFLAGS"  
-], AC_MSG_ERROR([need to have GiNaC installed]))
+PKG_CHECK_MODULES(SIMPLE, ginac >= 1.3.7)
 
 AC_OUTPUT(Makefile)
 @end example
 
-The only command in this which is not standard for automake
-is the @samp{AM_PATH_GINAC} macro.
+The @samp{PKG_CHECK_MODULES} macro does the following: If a GiNaC version
+greater or equal than 1.3.7 is found, then it defines @var{SIMPLE_CFLAGS}
+and @var{SIMPLE_LIBS}. Otherwise, it dies with the error message like
+@example
+configure: error: Package requirements (ginac >= 1.3.7) were not met:
+
+Requested 'ginac >= 1.3.7' but version of GiNaC is 1.3.5
 
-That command does the following: If a GiNaC version greater or equal
-than 0.7.0 is found, then it adds @env{$GINACLIB_LIBS} to @env{$LIBS}
-and @env{$GINACLIB_CPPFLAGS} to @env{$CPPFLAGS}. Otherwise, it dies with
-the error message `need to have GiNaC installed'
+Consider adjusting the PKG_CONFIG_PATH environment variable if you
+installed software in a non-standard prefix.
+
+Alternatively, you may set the environment variables SIMPLE_CFLAGS
+and SIMPLE_LIBS to avoid the need to call pkg-config.
+See the pkg-config man page for more details.
+@end example
 
 And the @file{Makefile.am}, which will be used to build the Makefile.
 
@@ -6854,31 +8756,24 @@ And the @file{Makefile.am}, which will be used to build the Makefile.
 ## Process this file with automake to produce Makefile.in
 bin_PROGRAMS = simple
 simple_SOURCES = simple.cpp
+simple_CPPFLAGS = $(SIMPLE_CFLAGS)
+simple_LDADD = $(SIMPLE_LIBS)
 @end example
 
 This @file{Makefile.am}, says that we are building a single executable,
 from a single source file @file{simple.cpp}. Since every program
-we are building uses GiNaC we simply added the GiNaC options
-to @env{$LIBS} and @env{$CPPFLAGS}, but in other circumstances, we might
-want to specify them on a per-program basis: for instance by
-adding the lines:
-
-@example
-simple_LDADD = $(GINACLIB_LIBS)
-INCLUDES = $(GINACLIB_CPPFLAGS)
-@end example
-
-to the @file{Makefile.am}.
+we are building uses GiNaC we could have simply added @var{SIMPLE_CFLAGS}
+to @var{CPPFLAGS} and @var{SIMPLE_LIBS} to @var{LIBS}. However, it is
+more flexible to specify libraries and complier options on a per-program
+basis.
 
 To try this example out, create a new directory and add the three
 files above to it.
 
-Now execute the following commands:
+Now execute the following command:
 
 @example
-$ automake --add-missing
-$ aclocal
-$ autoconf
+$ autoreconf -i
 @end example
 
 You now have a package that can be built in the normal fashion
@@ -6890,7 +8785,7 @@ $ make install
 @end example
 
 
-@node Bibliography, Concept Index, Example package, Top
+@node Bibliography, Concept index, Example package, Top
 @c    node-name, next, previous, up
 @appendix Bibliography
 
@@ -6935,9 +8830,9 @@ ISBN 3-540-66572-2, 2001, Springer, Heidelberg
 @end itemize
 
 
-@node Concept Index, , Bibliography, Top
+@node Concept index, , Bibliography, Top
 @c    node-name, next, previous, up
-@unnumbered Concept Index
+@unnumbered Concept index
 
 @printindex cp