+In this chapter the most important algorithms provided by GiNaC will be
+described. Some of them are implemented as functions on expressions,
+others are implemented as methods provided by expression objects. If
+they are methods, there exists a wrapper function around it, so you can
+alternatively call it in a functional way as shown in the simple
+example:
+
+@example
+ ...
+ cout << "As method: " << sin(1).evalf() << endl;
+ cout << "As function: " << evalf(sin(1)) << endl;
+ ...
+@end example
+
+@cindex @code{subs()}
+The general rule is that wherever methods accept one or more parameters
+(@var{arg1}, @var{arg2}, @dots{}) the order of arguments the function
+wrapper accepts is the same but preceded by the object to act on
+(@var{object}, @var{arg1}, @var{arg2}, @dots{}). This approach is the
+most natural one in an OO model but it may lead to confusion for MapleV
+users because where they would type @code{A:=x+1; subs(x=2,A);} GiNaC
+would require @code{A=x+1; subs(A,x==2);} (after proper declaration of
+@code{A} and @code{x}). On the other hand, since MapleV returns 3 on
+@code{A:=x^2+3; coeff(A,x,0);} (GiNaC: @code{A=pow(x,2)+3;
+coeff(A,x,0);}) it is clear that MapleV is not trying to be consistent
+here. Also, users of MuPAD will in most cases feel more comfortable
+with GiNaC's convention. All function wrappers are implemented
+as simple inline functions which just call the corresponding method and
+are only provided for users uncomfortable with OO who are dead set to
+avoid method invocations. Generally, nested function wrappers are much
+harder to read than a sequence of methods and should therefore be
+avoided if possible. On the other hand, not everything in GiNaC is a
+method on class @code{ex} and sometimes calling a function cannot be
+avoided.
+
+@menu
+* Information about expressions::
+* Numerical evaluation::
+* Substituting expressions::
+* Pattern matching and advanced substitutions::
+* Applying a function on subexpressions::
+* Visitors and tree traversal::
+* Polynomial arithmetic:: Working with polynomials.
+* Rational expressions:: Working with rational functions.
+* Symbolic differentiation::
+* Series expansion:: Taylor and Laurent expansion.
+* Symmetrization::
+* Built-in functions:: List of predefined mathematical functions.
+* Multiple polylogarithms::
+* Complex expressions::
+* Solving linear systems of equations::
+* Input/output:: Input and output of expressions.
+@end menu
+
+
+@node Information about expressions, Numerical evaluation, Methods and functions, Methods and functions
+@c node-name, next, previous, up
+@section Getting information about expressions
+
+@subsection Checking expression types
+@cindex @code{is_a<@dots{}>()}
+@cindex @code{is_exactly_a<@dots{}>()}
+@cindex @code{ex_to<@dots{}>()}
+@cindex Converting @code{ex} to other classes
+@cindex @code{info()}
+@cindex @code{return_type()}
+@cindex @code{return_type_tinfo()}
+
+Sometimes it's useful to check whether a given expression is a plain number,
+a sum, a polynomial with integer coefficients, or of some other specific type.
+GiNaC provides a couple of functions for this:
+
+@example
+bool is_a<T>(const ex & e);
+bool is_exactly_a<T>(const ex & e);
+bool ex::info(unsigned flag);
+unsigned ex::return_type() const;
+unsigned ex::return_type_tinfo() const;
+@end example
+
+When the test made by @code{is_a<T>()} returns true, it is safe to call
+one of the functions @code{ex_to<T>()}, where @code{T} is one of the
+class names (@xref{The class hierarchy}, for a list of all classes). For
+example, assuming @code{e} is an @code{ex}:
+
+@example
+@{
+ @dots{}
+ if (is_a<numeric>(e))
+ numeric n = ex_to<numeric>(e);
+ @dots{}
+@}
+@end example
+
+@code{is_a<T>(e)} allows you to check whether the top-level object of
+an expression @samp{e} is an instance of the GiNaC class @samp{T}
+(@xref{The class hierarchy}, for a list of all classes). This is most useful,
+e.g., for checking whether an expression is a number, a sum, or a product:
+
+@example
+@{
+ symbol x("x");
+ ex e1 = 42;
+ ex e2 = 4*x - 3;
+ is_a<numeric>(e1); // true
+ is_a<numeric>(e2); // false
+ is_a<add>(e1); // false
+ is_a<add>(e2); // true
+ is_a<mul>(e1); // false
+ is_a<mul>(e2); // false
+@}
+@end example
+
+In contrast, @code{is_exactly_a<T>(e)} allows you to check whether the
+top-level object of an expression @samp{e} is an instance of the GiNaC
+class @samp{T}, not including parent classes.
+
+The @code{info()} method is used for checking certain attributes of
+expressions. The possible values for the @code{flag} argument are defined
+in @file{ginac/flags.h}, the most important being explained in the following
+table:
+
+@cartouche
+@multitable @columnfractions .30 .70
+@item @strong{Flag} @tab @strong{Returns true if the object is@dots{}}
+@item @code{numeric}
+@tab @dots{}a number (same as @code{is_a<numeric>(...)})
+@item @code{real}
+@tab @dots{}a real number, symbol or constant (i.e. is not complex)
+@item @code{rational}
+@tab @dots{}an exact rational number (integers are rational, too)
+@item @code{integer}
+@tab @dots{}a (non-complex) integer
+@item @code{crational}
+@tab @dots{}an exact (complex) rational number (such as @math{2/3+7/2*I})
+@item @code{cinteger}
+@tab @dots{}a (complex) integer (such as @math{2-3*I})
+@item @code{positive}
+@tab @dots{}not complex and greater than 0
+@item @code{negative}
+@tab @dots{}not complex and less than 0
+@item @code{nonnegative}
+@tab @dots{}not complex and greater than or equal to 0
+@item @code{posint}
+@tab @dots{}an integer greater than 0
+@item @code{negint}
+@tab @dots{}an integer less than 0
+@item @code{nonnegint}
+@tab @dots{}an integer greater than or equal to 0
+@item @code{even}
+@tab @dots{}an even integer
+@item @code{odd}
+@tab @dots{}an odd integer
+@item @code{prime}
+@tab @dots{}a prime integer (probabilistic primality test)
+@item @code{relation}
+@tab @dots{}a relation (same as @code{is_a<relational>(...)})
+@item @code{relation_equal}
+@tab @dots{}a @code{==} relation
+@item @code{relation_not_equal}
+@tab @dots{}a @code{!=} relation
+@item @code{relation_less}
+@tab @dots{}a @code{<} relation
+@item @code{relation_less_or_equal}
+@tab @dots{}a @code{<=} relation
+@item @code{relation_greater}
+@tab @dots{}a @code{>} relation
+@item @code{relation_greater_or_equal}
+@tab @dots{}a @code{>=} relation
+@item @code{symbol}
+@tab @dots{}a symbol (same as @code{is_a<symbol>(...)})
+@item @code{list}
+@tab @dots{}a list (same as @code{is_a<lst>(...)})
+@item @code{polynomial}
+@tab @dots{}a polynomial (i.e. only consists of sums and products of numbers and symbols with positive integer powers)
+@item @code{integer_polynomial}
+@tab @dots{}a polynomial with (non-complex) integer coefficients
+@item @code{cinteger_polynomial}
+@tab @dots{}a polynomial with (possibly complex) integer coefficients (such as @math{2-3*I})
+@item @code{rational_polynomial}
+@tab @dots{}a polynomial with (non-complex) rational coefficients
+@item @code{crational_polynomial}
+@tab @dots{}a polynomial with (possibly complex) rational coefficients (such as @math{2/3+7/2*I})
+@item @code{rational_function}
+@tab @dots{}a rational function (@math{x+y}, @math{z/(x+y)})
+@item @code{algebraic}
+@tab @dots{}an algebraic object (@math{sqrt(2)}, @math{sqrt(x)-1})
+@end multitable
+@end cartouche
+
+To determine whether an expression is commutative or non-commutative and if
+so, with which other expressions it would commutate, you use the methods
+@code{return_type()} and @code{return_type_tinfo()}. @xref{Non-commutative objects},
+for an explanation of these.
+
+
+@subsection Accessing subexpressions
+@cindex container
+
+Many GiNaC classes, like @code{add}, @code{mul}, @code{lst}, and
+@code{function}, act as containers for subexpressions. For example, the
+subexpressions of a sum (an @code{add} object) are the individual terms,
+and the subexpressions of a @code{function} are the function's arguments.
+
+@cindex @code{nops()}
+@cindex @code{op()}
+GiNaC provides several ways of accessing subexpressions. The first way is to
+use the two methods
+
+@example
+size_t ex::nops();
+ex ex::op(size_t i);
+@end example
+
+@code{nops()} determines the number of subexpressions (operands) contained
+in the expression, while @code{op(i)} returns the @code{i}-th
+(0..@code{nops()-1}) subexpression. In the case of a @code{power} object,
+@code{op(0)} will return the basis and @code{op(1)} the exponent. For
+@code{indexed} objects, @code{op(0)} is the base expression and @code{op(i)},
+@math{i>0} are the indices.
+
+@cindex iterators
+@cindex @code{const_iterator}
+The second way to access subexpressions is via the STL-style random-access
+iterator class @code{const_iterator} and the methods
+
+@example
+const_iterator ex::begin();
+const_iterator ex::end();
+@end example
+
+@code{begin()} returns an iterator referring to the first subexpression;
+@code{end()} returns an iterator which is one-past the last subexpression.
+If the expression has no subexpressions, then @code{begin() == end()}. These
+iterators can also be used in conjunction with non-modifying STL algorithms.
+
+Here is an example that (non-recursively) prints the subexpressions of a
+given expression in three different ways:
+
+@example
+@{
+ ex e = ...
+
+ // with nops()/op()
+ for (size_t i = 0; i != e.nops(); ++i)
+ cout << e.op(i) << endl;
+
+ // with iterators
+ for (const_iterator i = e.begin(); i != e.end(); ++i)
+ cout << *i << endl;
+
+ // with iterators and STL copy()
+ std::copy(e.begin(), e.end(), std::ostream_iterator<ex>(cout, "\n"));
+@}
+@end example
+
+@cindex @code{const_preorder_iterator}
+@cindex @code{const_postorder_iterator}
+@code{op()}/@code{nops()} and @code{const_iterator} only access an
+expression's immediate children. GiNaC provides two additional iterator
+classes, @code{const_preorder_iterator} and @code{const_postorder_iterator},
+that iterate over all objects in an expression tree, in preorder or postorder,
+respectively. They are STL-style forward iterators, and are created with the
+methods
+
+@example
+const_preorder_iterator ex::preorder_begin();
+const_preorder_iterator ex::preorder_end();
+const_postorder_iterator ex::postorder_begin();
+const_postorder_iterator ex::postorder_end();
+@end example
+
+The following example illustrates the differences between
+@code{const_iterator}, @code{const_preorder_iterator}, and
+@code{const_postorder_iterator}:
+
+@example
+@{
+ symbol A("A"), B("B"), C("C");
+ ex e = lst(lst(A, B), C);
+
+ std::copy(e.begin(), e.end(),
+ std::ostream_iterator<ex>(cout, "\n"));
+ // @{A,B@}
+ // C
+
+ std::copy(e.preorder_begin(), e.preorder_end(),
+ std::ostream_iterator<ex>(cout, "\n"));
+ // @{@{A,B@},C@}
+ // @{A,B@}
+ // A
+ // B
+ // C
+
+ std::copy(e.postorder_begin(), e.postorder_end(),
+ std::ostream_iterator<ex>(cout, "\n"));
+ // A
+ // B
+ // @{A,B@}
+ // C
+ // @{@{A,B@},C@}
+@}
+@end example
+
+@cindex @code{relational} (class)
+Finally, the left-hand side and right-hand side expressions of objects of
+class @code{relational} (and only of these) can also be accessed with the
+methods
+
+@example
+ex ex::lhs();
+ex ex::rhs();
+@end example
+
+
+@subsection Comparing expressions
+@cindex @code{is_equal()}
+@cindex @code{is_zero()}
+
+Expressions can be compared with the usual C++ relational operators like
+@code{==}, @code{>}, and @code{<} but if the expressions contain symbols,
+the result is usually not determinable and the result will be @code{false},
+except in the case of the @code{!=} operator. You should also be aware that
+GiNaC will only do the most trivial test for equality (subtracting both
+expressions), so something like @code{(pow(x,2)+x)/x==x+1} will return
+@code{false}.
+
+Actually, if you construct an expression like @code{a == b}, this will be
+represented by an object of the @code{relational} class (@pxref{Relations})
+which is not evaluated until (explicitly or implicitly) cast to a @code{bool}.
+
+There are also two methods
+
+@example
+bool ex::is_equal(const ex & other);
+bool ex::is_zero();
+@end example
+
+for checking whether one expression is equal to another, or equal to zero,
+respectively. See also the method @code{ex::is_zero_matrix()},
+@pxref{Matrices}.
+
+
+@subsection Ordering expressions
+@cindex @code{ex_is_less} (class)
+@cindex @code{ex_is_equal} (class)
+@cindex @code{compare()}
+
+Sometimes it is necessary to establish a mathematically well-defined ordering
+on a set of arbitrary expressions, for example to use expressions as keys
+in a @code{std::map<>} container, or to bring a vector of expressions into
+a canonical order (which is done internally by GiNaC for sums and products).
+
+The operators @code{<}, @code{>} etc. described in the last section cannot
+be used for this, as they don't implement an ordering relation in the
+mathematical sense. In particular, they are not guaranteed to be
+antisymmetric: if @samp{a} and @samp{b} are different expressions, and
+@code{a < b} yields @code{false}, then @code{b < a} doesn't necessarily
+yield @code{true}.
+
+By default, STL classes and algorithms use the @code{<} and @code{==}
+operators to compare objects, which are unsuitable for expressions, but GiNaC
+provides two functors that can be supplied as proper binary comparison
+predicates to the STL:
+
+@example
+class ex_is_less : public std::binary_function<ex, ex, bool> @{
+public:
+ bool operator()(const ex &lh, const ex &rh) const;
+@};
+
+class ex_is_equal : public std::binary_function<ex, ex, bool> @{
+public:
+ bool operator()(const ex &lh, const ex &rh) const;
+@};
+@end example
+
+For example, to define a @code{map} that maps expressions to strings you
+have to use
+
+@example
+std::map<ex, std::string, ex_is_less> myMap;
+@end example
+
+Omitting the @code{ex_is_less} template parameter will introduce spurious
+bugs because the map operates improperly.
+
+Other examples for the use of the functors:
+
+@example
+std::vector<ex> v;
+// fill vector
+...
+
+// sort vector
+std::sort(v.begin(), v.end(), ex_is_less());
+
+// count the number of expressions equal to '1'
+unsigned num_ones = std::count_if(v.begin(), v.end(),
+ std::bind2nd(ex_is_equal(), 1));
+@end example
+
+The implementation of @code{ex_is_less} uses the member function
+
+@example
+int ex::compare(const ex & other) const;
+@end example
+
+which returns @math{0} if @code{*this} and @code{other} are equal, @math{-1}
+if @code{*this} sorts before @code{other}, and @math{1} if @code{*this} sorts
+after @code{other}.
+
+
+@node Numerical evaluation, Substituting expressions, Information about expressions, Methods and functions
+@c node-name, next, previous, up
+@section Numerical evaluation
+@cindex @code{evalf()}
+
+GiNaC keeps algebraic expressions, numbers and constants in their exact form.
+To evaluate them using floating-point arithmetic you need to call
+
+@example
+ex ex::evalf(int level = 0) const;
+@end example
+
+@cindex @code{Digits}
+The accuracy of the evaluation is controlled by the global object @code{Digits}
+which can be assigned an integer value. The default value of @code{Digits}
+is 17. @xref{Numbers}, for more information and examples.
+
+To evaluate an expression to a @code{double} floating-point number you can
+call @code{evalf()} followed by @code{numeric::to_double()}, like this:
+
+@example
+@{
+ // Approximate sin(x/Pi)
+ symbol x("x");
+ ex e = series(sin(x/Pi), x == 0, 6);
+
+ // Evaluate numerically at x=0.1
+ ex f = evalf(e.subs(x == 0.1));
+
+ // ex_to<numeric> is an unsafe cast, so check the type first
+ if (is_a<numeric>(f)) @{
+ double d = ex_to<numeric>(f).to_double();
+ cout << d << endl;
+ // -> 0.0318256
+ @} else
+ // error
+@}
+@end example
+
+
+@node Substituting expressions, Pattern matching and advanced substitutions, Numerical evaluation, Methods and functions
+@c node-name, next, previous, up
+@section Substituting expressions
+@cindex @code{subs()}
+
+Algebraic objects inside expressions can be replaced with arbitrary
+expressions via the @code{.subs()} method:
+
+@example
+ex ex::subs(const ex & e, unsigned options = 0);
+ex ex::subs(const exmap & m, unsigned options = 0);
+ex ex::subs(const lst & syms, const lst & repls, unsigned options = 0);
+@end example
+
+In the first form, @code{subs()} accepts a relational of the form
+@samp{object == expression} or a @code{lst} of such relationals:
+
+@example
+@{
+ symbol x("x"), y("y");
+
+ ex e1 = 2*x^2-4*x+3;
+ cout << "e1(7) = " << e1.subs(x == 7) << endl;
+ // -> 73
+
+ ex e2 = x*y + x;
+ cout << "e2(-2, 4) = " << e2.subs(lst(x == -2, y == 4)) << endl;
+ // -> -10
+@}
+@end example
+
+If you specify multiple substitutions, they are performed in parallel, so e.g.
+@code{subs(lst(x == y, y == x))} exchanges @samp{x} and @samp{y}.
+
+The second form of @code{subs()} takes an @code{exmap} object which is a
+pair associative container that maps expressions to expressions (currently
+implemented as a @code{std::map}). This is the most efficient one of the
+three @code{subs()} forms and should be used when the number of objects to
+be substituted is large or unknown.
+
+Using this form, the second example from above would look like this:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex e2 = x*y + x;
+
+ exmap m;
+ m[x] = -2;
+ m[y] = 4;
+ cout << "e2(-2, 4) = " << e2.subs(m) << endl;
+@}
+@end example
+
+The third form of @code{subs()} takes two lists, one for the objects to be
+replaced and one for the expressions to be substituted (both lists must
+contain the same number of elements). Using this form, you would write
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex e2 = x*y + x;
+
+ cout << "e2(-2, 4) = " << e2.subs(lst(x, y), lst(-2, 4)) << endl;
+@}
+@end example
+
+The optional last argument to @code{subs()} is a combination of
+@code{subs_options} flags. There are three options available:
+@code{subs_options::no_pattern} disables pattern matching, which makes
+large @code{subs()} operations significantly faster if you are not using
+patterns. The second option, @code{subs_options::algebraic} enables
+algebraic substitutions in products and powers.
+@xref{Pattern matching and advanced substitutions}, for more information
+about patterns and algebraic substitutions. The third option,
+@code{subs_options::no_index_renaming} disables the feature that dummy
+indices are renamed if the substitution could give a result in which a
+dummy index occurs more than two times. This is sometimes necessary if
+you want to use @code{subs()} to rename your dummy indices.
+
+@code{subs()} performs syntactic substitution of any complete algebraic
+object; it does not try to match sub-expressions as is demonstrated by the
+following example:
+
+@example
+@{
+ symbol x("x"), y("y"), z("z");
+
+ ex e1 = pow(x+y, 2);
+ cout << e1.subs(x+y == 4) << endl;
+ // -> 16
+
+ ex e2 = sin(x)*sin(y)*cos(x);
+ cout << e2.subs(sin(x) == cos(x)) << endl;
+ // -> cos(x)^2*sin(y)
+
+ ex e3 = x+y+z;
+ cout << e3.subs(x+y == 4) << endl;
+ // -> x+y+z
+ // (and not 4+z as one might expect)
+@}
+@end example
+
+A more powerful form of substitution using wildcards is described in the
+next section.
+
+
+@node Pattern matching and advanced substitutions, Applying a function on subexpressions, Substituting expressions, Methods and functions
+@c node-name, next, previous, up
+@section Pattern matching and advanced substitutions
+@cindex @code{wildcard} (class)
+@cindex Pattern matching
+
+GiNaC allows the use of patterns for checking whether an expression is of a
+certain form or contains subexpressions of a certain form, and for
+substituting expressions in a more general way.
+
+A @dfn{pattern} is an algebraic expression that optionally contains wildcards.
+A @dfn{wildcard} is a special kind of object (of class @code{wildcard}) that
+represents an arbitrary expression. Every wildcard has a @dfn{label} which is
+an unsigned integer number to allow having multiple different wildcards in a
+pattern. Wildcards are printed as @samp{$label} (this is also the way they
+are specified in @command{ginsh}). In C++ code, wildcard objects are created
+with the call
+
+@example
+ex wild(unsigned label = 0);
+@end example
+
+which is simply a wrapper for the @code{wildcard()} constructor with a shorter
+name.
+
+Some examples for patterns:
+
+@multitable @columnfractions .5 .5
+@item @strong{Constructed as} @tab @strong{Output as}
+@item @code{wild()} @tab @samp{$0}
+@item @code{pow(x,wild())} @tab @samp{x^$0}
+@item @code{atan2(wild(1),wild(2))} @tab @samp{atan2($1,$2)}
+@item @code{indexed(A,idx(wild(),3))} @tab @samp{A.$0}
+@end multitable
+
+Notes:
+
+@itemize @bullet
+@item Wildcards behave like symbols and are subject to the same algebraic
+ rules. E.g., @samp{$0+2*$0} is automatically transformed to @samp{3*$0}.
+@item As shown in the last example, to use wildcards for indices you have to
+ use them as the value of an @code{idx} object. This is because indices must
+ always be of class @code{idx} (or a subclass).
+@item Wildcards only represent expressions or subexpressions. It is not
+ possible to use them as placeholders for other properties like index
+ dimension or variance, representation labels, symmetry of indexed objects
+ etc.
+@item Because wildcards are commutative, it is not possible to use wildcards
+ as part of noncommutative products.
+@item A pattern does not have to contain wildcards. @samp{x} and @samp{x+y}
+ are also valid patterns.
+@end itemize
+
+@subsection Matching expressions
+@cindex @code{match()}
+The most basic application of patterns is to check whether an expression
+matches a given pattern. This is done by the function
+
+@example
+bool ex::match(const ex & pattern);
+bool ex::match(const ex & pattern, lst & repls);
+@end example
+
+This function returns @code{true} when the expression matches the pattern
+and @code{false} if it doesn't. If used in the second form, the actual
+subexpressions matched by the wildcards get returned in the @code{repls}
+object as a list of relations of the form @samp{wildcard == expression}.
+If @code{match()} returns false, the state of @code{repls} is undefined.
+For reproducible results, the list should be empty when passed to
+@code{match()}, but it is also possible to find similarities in multiple
+expressions by passing in the result of a previous match.
+
+The matching algorithm works as follows:
+
+@itemize
+@item A single wildcard matches any expression. If one wildcard appears
+ multiple times in a pattern, it must match the same expression in all
+ places (e.g. @samp{$0} matches anything, and @samp{$0*($0+1)} matches
+ @samp{x*(x+1)} but not @samp{x*(y+1)}).
+@item If the expression is not of the same class as the pattern, the match
+ fails (i.e. a sum only matches a sum, a function only matches a function,
+ etc.).
+@item If the pattern is a function, it only matches the same function
+ (i.e. @samp{sin($0)} matches @samp{sin(x)} but doesn't match @samp{exp(x)}).
+@item Except for sums and products, the match fails if the number of
+ subexpressions (@code{nops()}) is not equal to the number of subexpressions
+ of the pattern.
+@item If there are no subexpressions, the expressions and the pattern must
+ be equal (in the sense of @code{is_equal()}).
+@item Except for sums and products, each subexpression (@code{op()}) must
+ match the corresponding subexpression of the pattern.
+@end itemize
+
+Sums (@code{add}) and products (@code{mul}) are treated in a special way to
+account for their commutativity and associativity:
+
+@itemize
+@item If the pattern contains a term or factor that is a single wildcard,
+ this one is used as the @dfn{global wildcard}. If there is more than one
+ such wildcard, one of them is chosen as the global wildcard in a random
+ way.
+@item Every term/factor of the pattern, except the global wildcard, is
+ matched against every term of the expression in sequence. If no match is
+ found, the whole match fails. Terms that did match are not considered in
+ further matches.
+@item If there are no unmatched terms left, the match succeeds. Otherwise
+ the match fails unless there is a global wildcard in the pattern, in
+ which case this wildcard matches the remaining terms.
+@end itemize
+
+In general, having more than one single wildcard as a term of a sum or a
+factor of a product (such as @samp{a+$0+$1}) will lead to unpredictable or
+ambiguous results.
+
+Here are some examples in @command{ginsh} to demonstrate how it works (the
+@code{match()} function in @command{ginsh} returns @samp{FAIL} if the
+match fails, and the list of wildcard replacements otherwise):
+
+@example
+> match((x+y)^a,(x+y)^a);
+@{@}
+> match((x+y)^a,(x+y)^b);
+FAIL
+> match((x+y)^a,$1^$2);
+@{$1==x+y,$2==a@}
+> match((x+y)^a,$1^$1);
+FAIL
+> match((x+y)^(x+y),$1^$1);
+@{$1==x+y@}
+> match((x+y)^(x+y),$1^$2);
+@{$1==x+y,$2==x+y@}
+> match((a+b)*(a+c),($1+b)*($1+c));
+@{$1==a@}
+> match((a+b)*(a+c),(a+$1)*(a+$2));
+@{$1==b,$2==c@}
+ (Unpredictable. The result might also be [$1==c,$2==b].)
+> match((a+b)*(a+c),($1+$2)*($1+$3));
+ (The result is undefined. Due to the sequential nature of the algorithm
+ and the re-ordering of terms in GiNaC, the match for the first factor
+ may be @{$1==a,$2==b@} in which case the match for the second factor
+ succeeds, or it may be @{$1==b,$2==a@} which causes the second match to
+ fail.)
+> match(a*(x+y)+a*z+b,a*$1+$2);
+ (This is also ambiguous and may return either @{$1==z,$2==a*(x+y)+b@} or
+ @{$1=x+y,$2=a*z+b@}.)
+> match(a+b+c+d+e+f,c);
+FAIL
+> match(a+b+c+d+e+f,c+$0);
+@{$0==a+e+b+f+d@}
+> match(a+b+c+d+e+f,c+e+$0);
+@{$0==a+b+f+d@}
+> match(a+b,a+b+$0);
+@{$0==0@}
+> match(a*b^2,a^$1*b^$2);
+FAIL
+ (The matching is syntactic, not algebraic, and "a" doesn't match "a^$1"
+ even though a==a^1.)
+> match(x*atan2(x,x^2),$0*atan2($0,$0^2));
+@{$0==x@}
+> match(atan2(y,x^2),atan2(y,$0));
+@{$0==x^2@}
+@end example
+
+@subsection Matching parts of expressions
+@cindex @code{has()}
+A more general way to look for patterns in expressions is provided by the
+member function
+
+@example
+bool ex::has(const ex & pattern);
+@end example
+
+This function checks whether a pattern is matched by an expression itself or
+by any of its subexpressions.
+
+Again some examples in @command{ginsh} for illustration (in @command{ginsh},
+@code{has()} returns @samp{1} for @code{true} and @samp{0} for @code{false}):
+
+@example
+> has(x*sin(x+y+2*a),y);
+1
+> has(x*sin(x+y+2*a),x+y);
+0
+ (This is because in GiNaC, "x+y" is not a subexpression of "x+y+2*a" (which
+ has the subexpressions "x", "y" and "2*a".)
+> has(x*sin(x+y+2*a),x+y+$1);
+1
+ (But this is possible.)
+> has(x*sin(2*(x+y)+2*a),x+y);
+0
+ (This fails because "2*(x+y)" automatically gets converted to "2*x+2*y" of
+ which "x+y" is not a subexpression.)
+> has(x+1,x^$1);
+0
+ (Although x^1==x and x^0==1, neither "x" nor "1" are actually of the form
+ "x^something".)
+> has(4*x^2-x+3,$1*x);
+1
+> has(4*x^2+x+3,$1*x);
+0
+ (Another possible pitfall. The first expression matches because the term
+ "-x" has the form "(-1)*x" in GiNaC. To check whether a polynomial
+ contains a linear term you should use the coeff() function instead.)
+@end example
+
+@cindex @code{find()}
+The method
+
+@example
+bool ex::find(const ex & pattern, lst & found);
+@end example
+
+works a bit like @code{has()} but it doesn't stop upon finding the first
+match. Instead, it appends all found matches to the specified list. If there
+are multiple occurrences of the same expression, it is entered only once to
+the list. @code{find()} returns false if no matches were found (in
+@command{ginsh}, it returns an empty list):
+
+@example
+> find(1+x+x^2+x^3,x);
+@{x@}
+> find(1+x+x^2+x^3,y);
+@{@}
+> find(1+x+x^2+x^3,x^$1);
+@{x^3,x^2@}
+ (Note the absence of "x".)
+> expand((sin(x)+sin(y))*(a+b));
+sin(y)*a+sin(x)*b+sin(x)*a+sin(y)*b
+> find(%,sin($1));
+@{sin(y),sin(x)@}
+@end example
+
+@subsection Substituting expressions
+@cindex @code{subs()}
+Probably the most useful application of patterns is to use them for
+substituting expressions with the @code{subs()} method. Wildcards can be
+used in the search patterns as well as in the replacement expressions, where
+they get replaced by the expressions matched by them. @code{subs()} doesn't
+know anything about algebra; it performs purely syntactic substitutions.
+
+Some examples:
+
+@example
+> subs(a^2+b^2+(x+y)^2,$1^2==$1^3);
+b^3+a^3+(x+y)^3
+> subs(a^4+b^4+(x+y)^4,$1^2==$1^3);
+b^4+a^4+(x+y)^4
+> subs((a+b+c)^2,a+b==x);
+(a+b+c)^2
+> subs((a+b+c)^2,a+b+$1==x+$1);
+(x+c)^2
+> subs(a+2*b,a+b==x);
+a+2*b
+> subs(4*x^3-2*x^2+5*x-1,x==a);
+-1+5*a-2*a^2+4*a^3
+> subs(4*x^3-2*x^2+5*x-1,x^$0==a^$0);
+-1+5*x-2*a^2+4*a^3
+> subs(sin(1+sin(x)),sin($1)==cos($1));
+cos(1+cos(x))
+> expand(subs(a*sin(x+y)^2+a*cos(x+y)^2+b,cos($1)^2==1-sin($1)^2));
+a+b
+@end example
+
+The last example would be written in C++ in this way:
+
+@example
+@{
+ symbol a("a"), b("b"), x("x"), y("y");
+ e = a*pow(sin(x+y), 2) + a*pow(cos(x+y), 2) + b;
+ e = e.subs(pow(cos(wild()), 2) == 1-pow(sin(wild()), 2));
+ cout << e.expand() << endl;
+ // -> a+b
+@}
+@end example
+
+@subsection The option algebraic
+Both @code{has()} and @code{subs()} take an optional argument to pass them
+extra options. This section describes what happens if you give the former
+the option @code{has_options::algebraic} or the latter
+@code{subs_options::algebraic}. In that case the matching condition for
+powers and multiplications is changed in such a way that they become
+more intuitive. Intuition says that @code{x*y} is a part of @code{x*y*z}.
+If you use these options you will find that
+@code{(x*y*z).has(x*y, has_options::algebraic)} indeed returns true.
+Besides matching some of the factors of a product also powers match as
+often as is possible without getting negative exponents. For example
+@code{(x^5*y^2*z).subs(x^2*y^2==c, subs_options::algebraic)} will return
+@code{x*c^2*z}. This also works with negative powers:
+@code{(x^(-3)*y^(-2)*z).subs(1/(x*y)==c, subs_options::algebraic)} will
+return @code{x^(-1)*c^2*z}.
+
+@strong{Note:} this only works for multiplications
+and not for locating @code{x+y} within @code{x+y+z}.
+
+
+@node Applying a function on subexpressions, Visitors and tree traversal, Pattern matching and advanced substitutions, Methods and functions
+@c node-name, next, previous, up
+@section Applying a function on subexpressions
+@cindex tree traversal
+@cindex @code{map()}
+
+Sometimes you may want to perform an operation on specific parts of an
+expression while leaving the general structure of it intact. An example
+of this would be a matrix trace operation: the trace of a sum is the sum
+of the traces of the individual terms. That is, the trace should @dfn{map}
+on the sum, by applying itself to each of the sum's operands. It is possible
+to do this manually which usually results in code like this:
+
+@example
+ex calc_trace(ex e)
+@{
+ if (is_a<matrix>(e))
+ return ex_to<matrix>(e).trace();
+ else if (is_a<add>(e)) @{
+ ex sum = 0;
+ for (size_t i=0; i<e.nops(); i++)
+ sum += calc_trace(e.op(i));
+ return sum;
+ @} else if (is_a<mul>)(e)) @{
+ ...
+ @} else @{
+ ...
+ @}
+@}
+@end example
+
+This is, however, slightly inefficient (if the sum is very large it can take
+a long time to add the terms one-by-one), and its applicability is limited to
+a rather small class of expressions. If @code{calc_trace()} is called with
+a relation or a list as its argument, you will probably want the trace to
+be taken on both sides of the relation or of all elements of the list.
+
+GiNaC offers the @code{map()} method to aid in the implementation of such
+operations:
+
+@example
+ex ex::map(map_function & f) const;
+ex ex::map(ex (*f)(const ex & e)) const;
+@end example
+
+In the first (preferred) form, @code{map()} takes a function object that
+is subclassed from the @code{map_function} class. In the second form, it
+takes a pointer to a function that accepts and returns an expression.
+@code{map()} constructs a new expression of the same type, applying the
+specified function on all subexpressions (in the sense of @code{op()}),
+non-recursively.
+
+The use of a function object makes it possible to supply more arguments to
+the function that is being mapped, or to keep local state information.
+The @code{map_function} class declares a virtual function call operator
+that you can overload. Here is a sample implementation of @code{calc_trace()}
+that uses @code{map()} in a recursive fashion:
+
+@example
+struct calc_trace : public map_function @{
+ ex operator()(const ex &e)
+ @{
+ if (is_a<matrix>(e))
+ return ex_to<matrix>(e).trace();
+ else if (is_a<mul>(e)) @{
+ ...
+ @} else
+ return e.map(*this);
+ @}
+@};
+@end example
+
+This function object could then be used like this:
+
+@example
+@{
+ ex M = ... // expression with matrices
+ calc_trace do_trace;
+ ex tr = do_trace(M);
+@}
+@end example
+
+Here is another example for you to meditate over. It removes quadratic
+terms in a variable from an expanded polynomial:
+
+@example
+struct map_rem_quad : public map_function @{
+ ex var;
+ map_rem_quad(const ex & var_) : var(var_) @{@}
+
+ ex operator()(const ex & e)
+ @{
+ if (is_a<add>(e) || is_a<mul>(e))
+ return e.map(*this);
+ else if (is_a<power>(e) &&
+ e.op(0).is_equal(var) && e.op(1).info(info_flags::even))
+ return 0;
+ else
+ return e;
+ @}
+@};
+
+...
+
+@{
+ symbol x("x"), y("y");
+
+ ex e;
+ for (int i=0; i<8; i++)
+ e += pow(x, i) * pow(y, 8-i) * (i+1);
+ cout << e << endl;
+ // -> 4*y^5*x^3+5*y^4*x^4+8*y*x^7+7*y^2*x^6+2*y^7*x+6*y^3*x^5+3*y^6*x^2+y^8
+
+ map_rem_quad rem_quad(x);
+ cout << rem_quad(e) << endl;
+ // -> 4*y^5*x^3+8*y*x^7+2*y^7*x+6*y^3*x^5+y^8
+@}
+@end example
+
+@command{ginsh} offers a slightly different implementation of @code{map()}
+that allows applying algebraic functions to operands. The second argument
+to @code{map()} is an expression containing the wildcard @samp{$0} which
+acts as the placeholder for the operands:
+
+@example
+> map(a*b,sin($0));
+sin(a)*sin(b)
+> map(a+2*b,sin($0));
+sin(a)+sin(2*b)
+> map(@{a,b,c@},$0^2+$0);
+@{a^2+a,b^2+b,c^2+c@}
+@end example
+
+Note that it is only possible to use algebraic functions in the second
+argument. You can not use functions like @samp{diff()}, @samp{op()},
+@samp{subs()} etc. because these are evaluated immediately:
+
+@example
+> map(@{a,b,c@},diff($0,a));
+@{0,0,0@}
+ This is because "diff($0,a)" evaluates to "0", so the command is equivalent
+ to "map(@{a,b,c@},0)".
+@end example
+
+
+@node Visitors and tree traversal, Polynomial arithmetic, Applying a function on subexpressions, Methods and functions
+@c node-name, next, previous, up
+@section Visitors and tree traversal
+@cindex tree traversal
+@cindex @code{visitor} (class)
+@cindex @code{accept()}
+@cindex @code{visit()}
+@cindex @code{traverse()}
+@cindex @code{traverse_preorder()}
+@cindex @code{traverse_postorder()}
+
+Suppose that you need a function that returns a list of all indices appearing
+in an arbitrary expression. The indices can have any dimension, and for
+indices with variance you always want the covariant version returned.
+
+You can't use @code{get_free_indices()} because you also want to include
+dummy indices in the list, and you can't use @code{find()} as it needs
+specific index dimensions (and it would require two passes: one for indices
+with variance, one for plain ones).
+
+The obvious solution to this problem is a tree traversal with a type switch,
+such as the following:
+
+@example
+void gather_indices_helper(const ex & e, lst & l)
+@{
+ if (is_a<varidx>(e)) @{
+ const varidx & vi = ex_to<varidx>(e);
+ l.append(vi.is_covariant() ? vi : vi.toggle_variance());
+ @} else if (is_a<idx>(e)) @{
+ l.append(e);
+ @} else @{
+ size_t n = e.nops();
+ for (size_t i = 0; i < n; ++i)
+ gather_indices_helper(e.op(i), l);
+ @}
+@}
+
+lst gather_indices(const ex & e)
+@{
+ lst l;
+ gather_indices_helper(e, l);
+ l.sort();
+ l.unique();
+ return l;
+@}
+@end example
+
+This works fine but fans of object-oriented programming will feel
+uncomfortable with the type switch. One reason is that there is a possibility
+for subtle bugs regarding derived classes. If we had, for example, written
+
+@example
+ if (is_a<idx>(e)) @{
+ ...
+ @} else if (is_a<varidx>(e)) @{
+ ...
+@end example
+
+in @code{gather_indices_helper}, the code wouldn't have worked because the
+first line "absorbs" all classes derived from @code{idx}, including
+@code{varidx}, so the special case for @code{varidx} would never have been
+executed.
+
+Also, for a large number of classes, a type switch like the above can get
+unwieldy and inefficient (it's a linear search, after all).
+@code{gather_indices_helper} only checks for two classes, but if you had to
+write a function that required a different implementation for nearly
+every GiNaC class, the result would be very hard to maintain and extend.
+
+The cleanest approach to the problem would be to add a new virtual function
+to GiNaC's class hierarchy. In our example, there would be specializations
+for @code{idx} and @code{varidx} while the default implementation in
+@code{basic} performed the tree traversal. Unfortunately, in C++ it's
+impossible to add virtual member functions to existing classes without
+changing their source and recompiling everything. GiNaC comes with source,
+so you could actually do this, but for a small algorithm like the one
+presented this would be impractical.
+
+One solution to this dilemma is the @dfn{Visitor} design pattern,
+which is implemented in GiNaC (actually, Robert Martin's Acyclic Visitor
+variation, described in detail in
+@uref{http://objectmentor.com/publications/acv.pdf}). Instead of adding
+virtual functions to the class hierarchy to implement operations, GiNaC
+provides a single "bouncing" method @code{accept()} that takes an instance
+of a special @code{visitor} class and redirects execution to the one
+@code{visit()} virtual function of the visitor that matches the type of
+object that @code{accept()} was being invoked on.
+
+Visitors in GiNaC must derive from the global @code{visitor} class as well
+as from the class @code{T::visitor} of each class @code{T} they want to
+visit, and implement the member functions @code{void visit(const T &)} for
+each class.
+
+A call of
+
+@example
+void ex::accept(visitor & v) const;
+@end example
+
+will then dispatch to the correct @code{visit()} member function of the
+specified visitor @code{v} for the type of GiNaC object at the root of the
+expression tree (e.g. a @code{symbol}, an @code{idx} or a @code{mul}).
+
+Here is an example of a visitor:
+
+@example
+class my_visitor
+ : public visitor, // this is required
+ public add::visitor, // visit add objects
+ public numeric::visitor, // visit numeric objects
+ public basic::visitor // visit basic objects
+@{
+ void visit(const add & x)
+ @{ cout << "called with an add object" << endl; @}
+
+ void visit(const numeric & x)
+ @{ cout << "called with a numeric object" << endl; @}
+
+ void visit(const basic & x)
+ @{ cout << "called with a basic object" << endl; @}
+@};
+@end example
+
+which can be used as follows:
+
+@example
+...
+ symbol x("x");
+ ex e1 = 42;
+ ex e2 = 4*x-3;
+ ex e3 = 8*x;
+
+ my_visitor v;
+ e1.accept(v);
+ // prints "called with a numeric object"
+ e2.accept(v);
+ // prints "called with an add object"
+ e3.accept(v);
+ // prints "called with a basic object"
+...
+@end example
+
+The @code{visit(const basic &)} method gets called for all objects that are
+not @code{numeric} or @code{add} and acts as an (optional) default.
+
+From a conceptual point of view, the @code{visit()} methods of the visitor
+behave like a newly added virtual function of the visited hierarchy.
+In addition, visitors can store state in member variables, and they can
+be extended by deriving a new visitor from an existing one, thus building
+hierarchies of visitors.
+
+We can now rewrite our index example from above with a visitor:
+
+@example
+class gather_indices_visitor
+ : public visitor, public idx::visitor, public varidx::visitor
+@{
+ lst l;
+
+ void visit(const idx & i)
+ @{
+ l.append(i);
+ @}
+
+ void visit(const varidx & vi)
+ @{
+ l.append(vi.is_covariant() ? vi : vi.toggle_variance());
+ @}
+
+public:
+ const lst & get_result() // utility function
+ @{
+ l.sort();
+ l.unique();
+ return l;
+ @}
+@};
+@end example
+
+What's missing is the tree traversal. We could implement it in
+@code{visit(const basic &)}, but GiNaC has predefined methods for this:
+
+@example
+void ex::traverse_preorder(visitor & v) const;
+void ex::traverse_postorder(visitor & v) const;
+void ex::traverse(visitor & v) const;
+@end example
+
+@code{traverse_preorder()} visits a node @emph{before} visiting its
+subexpressions, while @code{traverse_postorder()} visits a node @emph{after}
+visiting its subexpressions. @code{traverse()} is a synonym for
+@code{traverse_preorder()}.
+
+Here is a new implementation of @code{gather_indices()} that uses the visitor
+and @code{traverse()}:
+
+@example
+lst gather_indices(const ex & e)
+@{
+ gather_indices_visitor v;
+ e.traverse(v);
+ return v.get_result();
+@}
+@end example
+
+Alternatively, you could use pre- or postorder iterators for the tree
+traversal:
+
+@example
+lst gather_indices(const ex & e)
+@{
+ gather_indices_visitor v;
+ for (const_preorder_iterator i = e.preorder_begin();
+ i != e.preorder_end(); ++i) @{
+ i->accept(v);
+ @}
+ return v.get_result();
+@}
+@end example
+
+
+@node Polynomial arithmetic, Rational expressions, Visitors and tree traversal, Methods and functions
+@c node-name, next, previous, up
+@section Polynomial arithmetic
+
+@subsection Testing whether an expression is a polynomial
+@cindex @code{is_polynomial()}
+
+Testing whether an expression is a polynomial in one or more variables
+can be done with the method
+@example
+bool ex::is_polynomial(const ex & vars) const;
+@end example
+In the case of more than
+one variable, the variables are given as a list.
+
+@example
+(x*y*sin(y)).is_polynomial(x) // Returns true.
+(x*y*sin(y)).is_polynomial(lst(x,y)) // Returns false.
+@end example
+
+@subsection Expanding and collecting
+@cindex @code{expand()}
+@cindex @code{collect()}
+@cindex @code{collect_common_factors()}
+
+A polynomial in one or more variables has many equivalent
+representations. Some useful ones serve a specific purpose. Consider
+for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
+21*y*z + 4*z^2} (written down here in output-style). It is equivalent
+to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}. Other
+representations are the recursive ones where one collects for exponents
+in one of the three variable. Since the factors are themselves
+polynomials in the remaining two variables the procedure can be
+repeated. In our example, two possibilities would be @math{(4*y + z)*x
++ 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
+x*z}.
+
+To bring an expression into expanded form, its method
+
+@example
+ex ex::expand(unsigned options = 0);
+@end example
+
+may be called. In our example above, this corresponds to @math{4*x*y +
+x*z + 20*y^2 + 21*y*z + 4*z^2}. Again, since the canonical form in
+GiNaC is not easy to guess you should be prepared to see different
+orderings of terms in such sums!
+
+Another useful representation of multivariate polynomials is as a
+univariate polynomial in one of the variables with the coefficients
+being polynomials in the remaining variables. The method
+@code{collect()} accomplishes this task:
+
+@example
+ex ex::collect(const ex & s, bool distributed = false);
+@end example
+
+The first argument to @code{collect()} can also be a list of objects in which
+case the result is either a recursively collected polynomial, or a polynomial
+in a distributed form with terms like @math{c*x1^e1*...*xn^en}, as specified
+by the @code{distributed} flag.
+
+Note that the original polynomial needs to be in expanded form (for the
+variables concerned) in order for @code{collect()} to be able to find the
+coefficients properly.
+
+The following @command{ginsh} transcript shows an application of @code{collect()}
+together with @code{find()}:
+
+@example
+> a=expand((sin(x)+sin(y))*(1+p+q)*(1+d));
+d*p*sin(x)+p*sin(x)+q*d*sin(x)+q*sin(y)+d*sin(x)+q*d*sin(y)+sin(y)+d*sin(y)
++q*sin(x)+d*sin(y)*p+sin(x)+sin(y)*p
+> collect(a,@{p,q@});
+d*sin(x)+(d*sin(x)+sin(y)+d*sin(y)+sin(x))*p
++(d*sin(x)+sin(y)+d*sin(y)+sin(x))*q+sin(y)+d*sin(y)+sin(x)
+> collect(a,find(a,sin($1)));
+(1+q+d+q*d+d*p+p)*sin(y)+(1+q+d+q*d+d*p+p)*sin(x)
+> collect(a,@{find(a,sin($1)),p,q@});
+(1+(1+d)*p+d+q*(1+d))*sin(x)+(1+(1+d)*p+d+q*(1+d))*sin(y)
+> collect(a,@{find(a,sin($1)),d@});
+(1+q+d*(1+q+p)+p)*sin(y)+(1+q+d*(1+q+p)+p)*sin(x)
+@end example
+
+Polynomials can often be brought into a more compact form by collecting
+common factors from the terms of sums. This is accomplished by the function
+
+@example
+ex collect_common_factors(const ex & e);
+@end example
+
+This function doesn't perform a full factorization but only looks for
+factors which are already explicitly present:
+
+@example
+> collect_common_factors(a*x+a*y);
+(x+y)*a
+> collect_common_factors(a*x^2+2*a*x*y+a*y^2);
+a*(2*x*y+y^2+x^2)
+> collect_common_factors(a*(b*(a+c)*x+b*((a+c)*x+(a+c)*y)*y));
+(c+a)*a*(x*y+y^2+x)*b
+@end example
+
+@subsection Degree and coefficients
+@cindex @code{degree()}
+@cindex @code{ldegree()}
+@cindex @code{coeff()}
+
+The degree and low degree of a polynomial can be obtained using the two
+methods
+
+@example
+int ex::degree(const ex & s);
+int ex::ldegree(const ex & s);
+@end example
+
+which also work reliably on non-expanded input polynomials (they even work
+on rational functions, returning the asymptotic degree). By definition, the
+degree of zero is zero. To extract a coefficient with a certain power from
+an expanded polynomial you use
+
+@example
+ex ex::coeff(const ex & s, int n);
+@end example
+
+You can also obtain the leading and trailing coefficients with the methods
+
+@example
+ex ex::lcoeff(const ex & s);
+ex ex::tcoeff(const ex & s);
+@end example
+
+which are equivalent to @code{coeff(s, degree(s))} and @code{coeff(s, ldegree(s))},
+respectively.
+
+An application is illustrated in the next example, where a multivariate
+polynomial is analyzed:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex PolyInp = 4*pow(x,3)*y + 5*x*pow(y,2) + 3*y
+ - pow(x+y,2) + 2*pow(y+2,2) - 8;
+ ex Poly = PolyInp.expand();
+
+ for (int i=Poly.ldegree(x); i<=Poly.degree(x); ++i) @{
+ cout << "The x^" << i << "-coefficient is "
+ << Poly.coeff(x,i) << endl;
+ @}
+ cout << "As polynomial in y: "
+ << Poly.collect(y) << endl;
+@}
+@end example
+
+When run, it returns an output in the following fashion:
+
+@example
+The x^0-coefficient is y^2+11*y
+The x^1-coefficient is 5*y^2-2*y
+The x^2-coefficient is -1
+The x^3-coefficient is 4*y
+As polynomial in y: -x^2+(5*x+1)*y^2+(-2*x+4*x^3+11)*y
+@end example
+
+As always, the exact output may vary between different versions of GiNaC
+or even from run to run since the internal canonical ordering is not
+within the user's sphere of influence.
+
+@code{degree()}, @code{ldegree()}, @code{coeff()}, @code{lcoeff()},
+@code{tcoeff()} and @code{collect()} can also be used to a certain degree
+with non-polynomial expressions as they not only work with symbols but with
+constants, functions and indexed objects as well:
+
+@example
+@{
+ symbol a("a"), b("b"), c("c"), x("x");
+ idx i(symbol("i"), 3);
+
+ ex e = pow(sin(x) - cos(x), 4);
+ cout << e.degree(cos(x)) << endl;
+ // -> 4
+ cout << e.expand().coeff(sin(x), 3) << endl;
+ // -> -4*cos(x)
+
+ e = indexed(a+b, i) * indexed(b+c, i);
+ e = e.expand(expand_options::expand_indexed);
+ cout << e.collect(indexed(b, i)) << endl;
+ // -> a.i*c.i+(a.i+c.i)*b.i+b.i^2
+@}
+@end example
+
+
+@subsection Polynomial division
+@cindex polynomial division
+@cindex quotient
+@cindex remainder
+@cindex pseudo-remainder
+@cindex @code{quo()}
+@cindex @code{rem()}
+@cindex @code{prem()}
+@cindex @code{divide()}
+
+The two functions
+
+@example
+ex quo(const ex & a, const ex & b, const ex & x);
+ex rem(const ex & a, const ex & b, const ex & x);
+@end example
+
+compute the quotient and remainder of univariate polynomials in the variable
+@samp{x}. The results satisfy @math{a = b*quo(a, b, x) + rem(a, b, x)}.
+
+The additional function
+
+@example
+ex prem(const ex & a, const ex & b, const ex & x);
+@end example
+
+computes the pseudo-remainder of @samp{a} and @samp{b} which satisfies
+@math{c*a = b*q + prem(a, b, x)}, where @math{c = b.lcoeff(x) ^ (a.degree(x) - b.degree(x) + 1)}.
+
+Exact division of multivariate polynomials is performed by the function
+
+@example
+bool divide(const ex & a, const ex & b, ex & q);
+@end example
+
+If @samp{b} divides @samp{a} over the rationals, this function returns @code{true}
+and returns the quotient in the variable @code{q}. Otherwise it returns @code{false}
+in which case the value of @code{q} is undefined.
+
+
+@subsection Unit, content and primitive part
+@cindex @code{unit()}
+@cindex @code{content()}
+@cindex @code{primpart()}
+@cindex @code{unitcontprim()}
+
+The methods
+
+@example
+ex ex::unit(const ex & x);
+ex ex::content(const ex & x);
+ex ex::primpart(const ex & x);
+ex ex::primpart(const ex & x, const ex & c);
+@end example
+
+return the unit part, content part, and primitive polynomial of a multivariate
+polynomial with respect to the variable @samp{x} (the unit part being the sign
+of the leading coefficient, the content part being the GCD of the coefficients,
+and the primitive polynomial being the input polynomial divided by the unit and
+content parts). The second variant of @code{primpart()} expects the previously
+calculated content part of the polynomial in @code{c}, which enables it to
+work faster in the case where the content part has already been computed. The
+product of unit, content, and primitive part is the original polynomial.
+
+Additionally, the method
+
+@example
+void ex::unitcontprim(const ex & x, ex & u, ex & c, ex & p);
+@end example
+
+computes the unit, content, and primitive parts in one go, returning them
+in @code{u}, @code{c}, and @code{p}, respectively.
+
+
+@subsection GCD, LCM and resultant
+@cindex GCD
+@cindex LCM
+@cindex @code{gcd()}
+@cindex @code{lcm()}
+
+The functions for polynomial greatest common divisor and least common
+multiple have the synopsis
+
+@example
+ex gcd(const ex & a, const ex & b);
+ex lcm(const ex & a, const ex & b);
+@end example
+
+The functions @code{gcd()} and @code{lcm()} accept two expressions
+@code{a} and @code{b} as arguments and return a new expression, their
+greatest common divisor or least common multiple, respectively. If the
+polynomials @code{a} and @code{b} are coprime @code{gcd(a,b)} returns 1
+and @code{lcm(a,b)} returns the product of @code{a} and @code{b}. Note that all
+the coefficients must be rationals.
+
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+int main()
+@{
+ symbol x("x"), y("y"), z("z");
+ ex P_a = 4*x*y + x*z + 20*pow(y, 2) + 21*y*z + 4*pow(z, 2);
+ ex P_b = x*y + 3*x*z + 5*pow(y, 2) + 19*y*z + 12*pow(z, 2);
+
+ ex P_gcd = gcd(P_a, P_b);
+ // x + 5*y + 4*z
+ ex P_lcm = lcm(P_a, P_b);
+ // 4*x*y^2 + 13*y*x*z + 20*y^3 + 81*y^2*z + 67*y*z^2 + 3*x*z^2 + 12*z^3
+@}
+@end example
+
+@cindex resultant
+@cindex @code{resultant()}
+
+The resultant of two expressions only makes sense with polynomials.
+It is always computed with respect to a specific symbol within the
+expressions. The function has the interface
+
+@example
+ex resultant(const ex & a, const ex & b, const ex & s);
+@end example
+
+Resultants are symmetric in @code{a} and @code{b}. The following example
+computes the resultant of two expressions with respect to @code{x} and
+@code{y}, respectively:
+
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+int main()
+@{
+ symbol x("x"), y("y");
+
+ ex e1 = x+pow(y,2), e2 = 2*pow(x,3)-1; // x+y^2, 2*x^3-1
+ ex r;
+
+ r = resultant(e1, e2, x);
+ // -> 1+2*y^6
+ r = resultant(e1, e2, y);
+ // -> 1-4*x^3+4*x^6
+@}
+@end example
+
+@subsection Square-free decomposition
+@cindex square-free decomposition
+@cindex factorization
+@cindex @code{sqrfree()}
+
+GiNaC still lacks proper factorization support. Some form of
+factorization is, however, easily implemented by noting that factors
+appearing in a polynomial with power two or more also appear in the
+derivative and hence can easily be found by computing the GCD of the
+original polynomial and its derivatives. Any decent system has an
+interface for this so called square-free factorization. So we provide
+one, too:
+@example
+ex sqrfree(const ex & a, const lst & l = lst());
+@end example
+Here is an example that by the way illustrates how the exact form of the
+result may slightly depend on the order of differentiation, calling for
+some care with subsequent processing of the result:
+@example
+ ...
+ symbol x("x"), y("y");
+ ex BiVarPol = expand(pow(2-2*y,3) * pow(1+x*y,2) * pow(x-2*y,2) * (x+y));
+
+ cout << sqrfree(BiVarPol, lst(x,y)) << endl;
+ // -> 8*(1-y)^3*(y*x^2-2*y+x*(1-2*y^2))^2*(y+x)
+
+ cout << sqrfree(BiVarPol, lst(y,x)) << endl;
+ // -> 8*(1-y)^3*(-y*x^2+2*y+x*(-1+2*y^2))^2*(y+x)
+
+ cout << sqrfree(BiVarPol) << endl;
+ // -> depending on luck, any of the above
+ ...
+@end example
+Note also, how factors with the same exponents are not fully factorized
+with this method.
+
+
+@node Rational expressions, Symbolic differentiation, Polynomial arithmetic, Methods and functions
+@c node-name, next, previous, up
+@section Rational expressions
+
+@subsection The @code{normal} method
+@cindex @code{normal()}
+@cindex simplification
+@cindex temporary replacement
+
+Some basic form of simplification of expressions is called for frequently.
+GiNaC provides the method @code{.normal()}, which converts a rational function
+into an equivalent rational function of the form @samp{numerator/denominator}
+where numerator and denominator are coprime. If the input expression is already
+a fraction, it just finds the GCD of numerator and denominator and cancels it,
+otherwise it performs fraction addition and multiplication.
+
+@code{.normal()} can also be used on expressions which are not rational functions
+as it will replace all non-rational objects (like functions or non-integer
+powers) by temporary symbols to bring the expression to the domain of rational
+functions before performing the normalization, and re-substituting these
+symbols afterwards. This algorithm is also available as a separate method
+@code{.to_rational()}, described below.
+
+This means that both expressions @code{t1} and @code{t2} are indeed
+simplified in this little code snippet:
+
+@example
+@{
+ symbol x("x");
+ ex t1 = (pow(x,2) + 2*x + 1)/(x + 1);
+ ex t2 = (pow(sin(x),2) + 2*sin(x) + 1)/(sin(x) + 1);
+ std::cout << "t1 is " << t1.normal() << std::endl;
+ std::cout << "t2 is " << t2.normal() << std::endl;
+@}
+@end example
+
+Of course this works for multivariate polynomials too, so the ratio of
+the sample-polynomials from the section about GCD and LCM above would be
+normalized to @code{P_a/P_b} = @code{(4*y+z)/(y+3*z)}.
+
+
+@subsection Numerator and denominator
+@cindex numerator
+@cindex denominator
+@cindex @code{numer()}
+@cindex @code{denom()}
+@cindex @code{numer_denom()}
+
+The numerator and denominator of an expression can be obtained with
+
+@example
+ex ex::numer();
+ex ex::denom();
+ex ex::numer_denom();
+@end example
+
+These functions will first normalize the expression as described above and
+then return the numerator, denominator, or both as a list, respectively.
+If you need both numerator and denominator, calling @code{numer_denom()} is
+faster than using @code{numer()} and @code{denom()} separately.
+
+
+@subsection Converting to a polynomial or rational expression
+@cindex @code{to_polynomial()}
+@cindex @code{to_rational()}
+
+Some of the methods described so far only work on polynomials or rational
+functions. GiNaC provides a way to extend the domain of these functions to
+general expressions by using the temporary replacement algorithm described
+above. You do this by calling
+
+@example
+ex ex::to_polynomial(exmap & m);
+ex ex::to_polynomial(lst & l);
+@end example
+or
+@example
+ex ex::to_rational(exmap & m);
+ex ex::to_rational(lst & l);
+@end example
+
+on the expression to be converted. The supplied @code{exmap} or @code{lst}
+will be filled with the generated temporary symbols and their replacement
+expressions in a format that can be used directly for the @code{subs()}
+method. It can also already contain a list of replacements from an earlier
+application of @code{.to_polynomial()} or @code{.to_rational()}, so it's
+possible to use it on multiple expressions and get consistent results.
+
+The difference between @code{.to_polynomial()} and @code{.to_rational()}
+is probably best illustrated with an example:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex a = 2*x/sin(x) - y/(3*sin(x));
+ cout << a << endl;
+
+ lst lp;
+ ex p = a.to_polynomial(lp);
+ cout << " = " << p << "\n with " << lp << endl;
+ // = symbol3*symbol2*y+2*symbol2*x
+ // with @{symbol2==sin(x)^(-1),symbol3==-1/3@}
+
+ lst lr;
+ ex r = a.to_rational(lr);
+ cout << " = " << r << "\n with " << lr << endl;
+ // = -1/3*symbol4^(-1)*y+2*symbol4^(-1)*x
+ // with @{symbol4==sin(x)@}
+@}
+@end example
+
+The following more useful example will print @samp{sin(x)-cos(x)}:
+
+@example
+@{
+ symbol x("x");
+ ex a = pow(sin(x), 2) - pow(cos(x), 2);
+ ex b = sin(x) + cos(x);
+ ex q;
+ exmap m;
+ divide(a.to_polynomial(m), b.to_polynomial(m), q);
+ cout << q.subs(m) << endl;
+@}
+@end example
+
+
+@node Symbolic differentiation, Series expansion, Rational expressions, Methods and functions
+@c node-name, next, previous, up
+@section Symbolic differentiation
+@cindex differentiation
+@cindex @code{diff()}
+@cindex chain rule
+@cindex product rule
+
+GiNaC's objects know how to differentiate themselves. Thus, a
+polynomial (class @code{add}) knows that its derivative is the sum of
+the derivatives of all the monomials:
+
+@example
+@{
+ symbol x("x"), y("y"), z("z");
+ ex P = pow(x, 5) + pow(x, 2) + y;
+
+ cout << P.diff(x,2) << endl;
+ // -> 20*x^3 + 2
+ cout << P.diff(y) << endl; // 1
+ // -> 1
+ cout << P.diff(z) << endl; // 0
+ // -> 0
+@}
+@end example
+
+If a second integer parameter @var{n} is given, the @code{diff} method
+returns the @var{n}th derivative.
+
+If @emph{every} object and every function is told what its derivative
+is, all derivatives of composed objects can be calculated using the
+chain rule and the product rule. Consider, for instance the expression