+ 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
+@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==c,$2==b@}
+ (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 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))
+
+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
+
+
+@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
+
+
+@node Polynomial Arithmetic, Rational Expressions, Visitors and Tree Traversal, Methods and Functions
+@c node-name, next, previous, up
+@section Polynomial arithmetic
+
+@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