+@cindex @code{tensor} (class)
+@subsection Predefined tensors
+
+Some frequently used special tensors such as the delta, epsilon and metric
+tensors are predefined in GiNaC. They have special properties when
+contracted with other tensor expressions and some of them have constant
+matrix representations (they will evaluate to a number when numeric
+indices are specified).
+
+@cindex @code{delta_tensor()}
+@subsubsection Delta tensor
+
+The delta tensor takes two indices, is symmetric and has the matrix
+representation @code{diag(1,1,1,...)}. It is constructed by the function
+@code{delta_tensor()}:
+
+@example
+@{
+ symbol A("A"), B("B");
+
+ idx i(symbol("i"), 3), j(symbol("j"), 3),
+ k(symbol("k"), 3), l(symbol("l"), 3);
+
+ ex e = indexed(A, i, j) * indexed(B, k, l)
+ * delta_tensor(i, k) * delta_tensor(j, l) << endl;
+ cout << e.simplify_indexed() << endl;
+ // -> B.i.j*A.i.j
+
+ cout << delta_tensor(i, i) << endl;
+ // -> 3
+@}
+@end example
+
+@cindex @code{metric_tensor()}
+@subsubsection General metric tensor
+
+The function @code{metric_tensor()} creates a general symmetric metric
+tensor with two indices that can be used to raise/lower tensor indices. The
+metric tensor is denoted as @samp{g} in the output and if its indices are of
+mixed variance it is automatically replaced by a delta tensor:
+
+@example
+@{
+ symbol A("A");
+
+ varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
+
+ ex e = metric_tensor(mu, nu) * indexed(A, nu.toggle_variance(), rho);
+ cout << e.simplify_indexed() << endl;
+ // -> A~mu~rho
+
+ e = delta_tensor(mu, nu.toggle_variance()) * metric_tensor(nu, rho);
+ cout << e.simplify_indexed() << endl;
+ // -> g~mu~rho
+
+ e = metric_tensor(mu.toggle_variance(), nu.toggle_variance())
+ * metric_tensor(nu, rho);
+ cout << e.simplify_indexed() << endl;
+ // -> delta.mu~rho
+
+ e = metric_tensor(nu.toggle_variance(), rho.toggle_variance())
+ * metric_tensor(mu, nu) * (delta_tensor(mu.toggle_variance(), rho)
+ + indexed(A, mu.toggle_variance(), rho));
+ cout << e.simplify_indexed() << endl;
+ // -> 4+A.rho~rho
+@}
+@end example
+
+@cindex @code{lorentz_g()}
+@subsubsection Minkowski metric tensor
+
+The Minkowski metric tensor is a special metric tensor with a constant
+matrix representation which is either @code{diag(1, -1, -1, ...)} (negative
+signature, the default) or @code{diag(-1, 1, 1, ...)} (positive signature).
+It is created with the function @code{lorentz_g()} (although it is output as
+@samp{eta}):
+
+@example
+@{
+ varidx mu(symbol("mu"), 4);
+
+ e = delta_tensor(varidx(0, 4), mu.toggle_variance())
+ * lorentz_g(mu, varidx(0, 4)); // negative signature
+ cout << e.simplify_indexed() << endl;
+ // -> 1
+
+ e = delta_tensor(varidx(0, 4), mu.toggle_variance())
+ * lorentz_g(mu, varidx(0, 4), true); // positive signature
+ cout << e.simplify_indexed() << endl;
+ // -> -1
+@}
+@end example
+
+@cindex @code{spinor_metric()}
+@subsubsection Spinor metric tensor
+
+The function @code{spinor_metric()} creates an antisymmetric tensor with
+two indices that is used to raise/lower indices of 2-component spinors.
+It is output as @samp{eps}:
+
+@example
+@{
+ symbol psi("psi");
+
+ spinidx A(symbol("A")), B(symbol("B")), C(symbol("C"));
+ ex A_co = A.toggle_variance(), B_co = B.toggle_variance();
+
+ e = spinor_metric(A, B) * indexed(psi, B_co);
+ cout << e.simplify_indexed() << endl;
+ // -> psi~A
+
+ e = spinor_metric(A, B) * indexed(psi, A_co);
+ cout << e.simplify_indexed() << endl;
+ // -> -psi~B
+
+ e = spinor_metric(A_co, B_co) * indexed(psi, B);
+ cout << e.simplify_indexed() << endl;
+ // -> -psi.A
+
+ e = spinor_metric(A_co, B_co) * indexed(psi, A);
+ cout << e.simplify_indexed() << endl;
+ // -> psi.B
+
+ e = spinor_metric(A_co, B_co) * spinor_metric(A, B);
+ cout << e.simplify_indexed() << endl;
+ // -> 2
+
+ e = spinor_metric(A_co, B_co) * spinor_metric(B, C);
+ cout << e.simplify_indexed() << endl;
+ // -> -delta.A~C
+@}
+@end example
+
+The matrix representation of the spinor metric is @code{[[ [[ 0, 1 ]], [[ -1, 0 ]] ]]}.
+
+@cindex @code{epsilon_tensor()}
+@cindex @code{lorentz_eps()}
+@subsubsection Epsilon tensor
+
+The epsilon tensor is totally antisymmetric, its number of indices is equal
+to the dimension of the index space (the indices must all be of the same
+numeric dimension), and @samp{eps.1.2.3...} (resp. @samp{eps~0~1~2...}) is
+defined to be 1. Its behaviour with indices that have a variance also
+depends on the signature of the metric. Epsilon tensors are output as
+@samp{eps}.
+
+There are three functions defined to create epsilon tensors in 2, 3 and 4
+dimensions:
+
+@example
+ex epsilon_tensor(const ex & i1, const ex & i2);
+ex epsilon_tensor(const ex & i1, const ex & i2, const ex & i3);
+ex lorentz_eps(const ex & i1, const ex & i2, const ex & i3, const ex & i4, bool pos_sig = false);
+@end example
+
+The first two functions create an epsilon tensor in 2 or 3 Euclidean
+dimensions, the last function creates an epsilon tensor in a 4-dimensional
+Minkowski space (the last @code{bool} argument specifies whether the metric
+has negative or positive signature, as in the case of the Minkowski metric
+tensor).
+
+@subsection Linear algebra
+
+The @code{matrix} class can be used with indices to do some simple linear
+algebra (linear combinations and products of vectors and matrices, traces
+and scalar products):
+
+@example
+@{
+ idx i(symbol("i"), 2), j(symbol("j"), 2);
+ symbol x("x"), y("y");
+
+ matrix A(2, 2, lst(1, 2, 3, 4)), X(2, 1, lst(x, y));
+
+ cout << indexed(A, i, i) << endl;
+ // -> 5
+
+ ex e = indexed(A, i, j) * indexed(X, j);
+ cout << e.simplify_indexed() << endl;
+ // -> [[ [[2*y+x]], [[4*y+3*x]] ]].i
+
+ e = indexed(A, i, j) * indexed(X, i) + indexed(X, j) * 2;
+ cout << e.simplify_indexed() << endl;
+ // -> [[ [[3*y+3*x,6*y+2*x]] ]].j
+@}
+@end example
+
+You can of course obtain the same results with the @code{matrix::add()},
+@code{matrix::mul()} and @code{matrix::trace()} methods but with indices you
+don't have to worry about transposing matrices.
+
+Matrix indices always start at 0 and their dimension must match the number
+of rows/columns of the matrix. Matrices with one row or one column are
+vectors and can have one or two indices (it doesn't matter whether it's a
+row or a column vector). Other matrices must have two indices.
+
+You should be careful when using indices with variance on matrices. GiNaC
+doesn't look at the variance and doesn't know that @samp{F~mu~nu} and
+@samp{F.mu.nu} are different matrices. In this case you should use only
+one form for @samp{F} and explicitly multiply it with a matrix representation
+of the metric tensor.
+
+
+@node Non-commutative objects, Methods and Functions, Indexed objects, Basic Concepts
+@c node-name, next, previous, up
+@section Non-commutative objects
+
+GiNaC is equipped to handle certain non-commutative algebras. Three classes of
+non-commutative objects are built-in which are mostly of use in high energy
+physics:
+
+@itemize
+@item Clifford (Dirac) algebra (class @code{clifford})
+@item su(3) Lie algebra (class @code{color})
+@item Matrices (unindexed) (class @code{matrix})
+@end itemize
+
+The @code{clifford} and @code{color} classes are subclasses of
+@code{indexed} because the elements of these algebras ususally carry
+indices.
+
+Unlike most computer algebra systems, GiNaC does not primarily provide an
+operator (often denoted @samp{&*}) for representing inert products of
+arbitrary objects. Rather, non-commutativity in GiNaC is a property of the
+classes of objects involved, and non-commutative products are formed with
+the usual @samp{*} operator, as are ordinary products. GiNaC is capable of
+figuring out by itself which objects commute and will group the factors
+by their class. Consider this example:
+
+@example
+ ...
+ varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
+ idx a(symbol("a"), 8), b(symbol("b"), 8);
+ ex e = -dirac_gamma(mu) * (2*color_T(a)) * 8 * color_T(b) * dirac_gamma(nu);
+ cout << e << endl;
+ // -> -16*(gamma~mu*gamma~nu)*(T.a*T.b)
+ ...
+@end example
+
+As can be seen, GiNaC pulls out the overall commutative factor @samp{-16} and
+groups the non-commutative factors (the gammas and the su(3) generators)
+together while preserving the order of factors within each class (because
+Clifford objects commute with color objects). The resulting expression is a
+@emph{commutative} product with two factors that are themselves non-commutative
+products (@samp{gamma~mu*gamma~nu} and @samp{T.a*T.b}). For clarification,
+parentheses are placed around the non-commutative products in the output.
+
+@cindex @code{ncmul} (class)
+Non-commutative products are internally represented by objects of the class
+@code{ncmul}, as opposed to commutative products which are handled by the
+@code{mul} class. You will normally not have to worry about this distinction,
+though.
+
+The advantage of this approach is that you never have to worry about using
+(or forgetting to use) a special operator when constructing non-commutative
+expressions. Also, non-commutative products in GiNaC are more intelligent
+than in other computer algebra systems; they can, for example, automatically
+canonicalize themselves according to rules specified in the implementation
+of the non-commutative classes. The drawback is that to work with other than
+the built-in algebras you have to implement new classes yourself. 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.
+
+@cindex @code{return_type()}
+@cindex @code{return_type_tinfo()}
+Information about the commutativity of an object or expression can be
+obtained with the two member functions
+
+@example
+unsigned ex::return_type(void) const;
+unsigned ex::return_type_tinfo(void) const;
+@end example
+
+The @code{return_type()} function returns one of three values (defined in
+the header file @file{flags.h}), corresponding to three categories of
+expressions in GiNaC:
+
+@itemize
+@item @code{return_types::commutative}: Commutes 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
+ 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
+ @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.
+
+Here are a couple of examples:
+
+@cartouche
+@multitable @columnfractions 0.33 0.33 0.34
+@item @strong{Expression} @tab @strong{@code{return_type()}} @tab @strong{@code{return_type_tinfo()}}
+@item @code{42} @tab @code{commutative} @tab -
+@item @code{2*x-y} @tab @code{commutative} @tab -
+@item @code{dirac_ONE()} @tab @code{noncommutative} @tab @code{TINFO_clifford}
+@item @code{dirac_gamma(mu)*dirac_gamma(nu)} @tab @code{noncommutative} @tab @code{TINFO_clifford}
+@item @code{2*color_T(a)} @tab @code{noncommutative} @tab @code{TINFO_color}
+@item @code{dirac_ONE()*color_T(a)} @tab @code{noncommutative_composite} @tab -
+@end multitable
+@end cartouche
+
+Note: the @code{return_type_tinfo()} of Clifford objects is only equal to
+@code{TINFO_clifford} for objects with a representation label of zero.
+Other representation labels yield a different @code{return_type_tinfo()},
+but it's the same for any two objects with the same label. This is also true
+for color objects.
+
+As a last note, powers of non-commutative objects are not allowed in GiNaC.
+You can use the function
+
+@example
+ex ncpow(const ex & basis, unsigned exponent)
+@end example
+
+instead which returns an expanded product (e.g. @code{ncpow(a*b, 2)} becomes
+@samp{a*b*a*b} if @samp{a} and @samp{b} are non-commutative expressions).
+
+
+@cindex @code{clifford} (class)
+@subsection Clifford algebra
+
+@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
+
+@example
+ex dirac_gamma(const ex & mu, unsigned char rl = 0);
+@end example
+
+which takes two arguments: the index and a @dfn{representation label} in the
+range 0 to 255 which is used to distinguish elements of different Clifford
+algebras (this is also called a @dfn{spin line index}). Gammas with different
+labels commute with each other. The dimension of the index can be 4 or (in
+the framework of dimensional regularization) any symbolic value. Spinor
+indices on Dirac gammas are not supported in GiNaC.
+
+@cindex @code{dirac_ONE()}
+The unity element of a Clifford algebra is constructed by
+
+@example
+ex dirac_ONE(unsigned char rl = 0);
+@end example
+
+@cindex @code{dirac_gamma5()}
+and there's a special element @samp{gamma5} that commutes with all other
+gammas and in 4 dimensions equals @samp{gamma~0 gamma~1 gamma~2 gamma~3},
+provided by
+
+@example
+ex dirac_gamma5(unsigned char rl = 0);
+@end example
+
+@cindex @code{dirac_gamma6()}
+@cindex @code{dirac_gamma7()}
+The two additional functions
+
+@example
+ex dirac_gamma6(unsigned char rl = 0);
+ex dirac_gamma7(unsigned char rl = 0);
+@end example
+
+return @code{dirac_ONE(rl) + dirac_gamma5(rl)} and @code{dirac_ONE(rl) - dirac_gamma5(rl)},
+respectively.
+
+@cindex @code{dirac_slash()}
+Finally, the function
+
+@example
+ex dirac_slash(const ex & e, const ex & dim, unsigned char rl = 0);
+@end example
+
+creates a term of the form @samp{e.mu gamma~mu} with a new and unique index
+whose dimension is given by the @code{dim} argument.
+
+In products of dirac gammas, superfluous unity elements are automatically
+removed, squares are replaced by their values and @samp{gamma5} is
+anticommuted to the front. The @code{simplify_indexed()} function performs
+contractions in gamma strings, for example
+
+@example
+@{
+ ...
+ symbol a("a"), b("b"), D("D");
+ varidx mu(symbol("mu"), D);
+ ex e = dirac_gamma(mu) * dirac_slash(a, D)
+ * dirac_gamma(mu.toggle_variance());
+ cout << e << endl;
+ // -> (gamma~mu*gamma~symbol10*gamma.mu)*a.symbol10
+ e = e.simplify_indexed();
+ cout << e << endl;
+ // -> -gamma~symbol10*a.symbol10*D+2*gamma~symbol10*a.symbol10
+ cout << e.subs(D == 4) << endl;
+ // -> -2*gamma~symbol10*a.symbol10
+ // [ == -2 * dirac_slash(a, D) ]
+ ...
+@}
+@end example
+
+@cindex @code{dirac_trace()}
+To calculate the trace of an expression containing strings of Dirac gammas
+you use the function
+
+@example
+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
+@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}.
+
+The value of the trace itself is also usually different in 4 and in
+@math{D != 4} dimensions:
+
+@example
+@{
+ // 4 dimensions
+ varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4), rho(symbol("rho"), 4);
+ ex e = dirac_gamma(mu) * dirac_gamma(nu) *
+ dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
+ cout << dirac_trace(e).simplify_indexed() << endl;
+ // -> -8*eta~rho~nu
+@}
+...
+@{
+ // D dimensions
+ symbol D("D");
+ varidx mu(symbol("mu"), D), nu(symbol("nu"), D), rho(symbol("rho"), D);
+ ex e = dirac_gamma(mu) * dirac_gamma(nu) *
+ dirac_gamma(mu.toggle_variance()) * dirac_gamma(rho);
+ cout << dirac_trace(e).simplify_indexed() << endl;
+ // -> 8*eta~rho~nu-4*eta~rho~nu*D
+@}
+@end example
+
+Here is an example for using @code{dirac_trace()} to compute a value that
+appears in the calculation of the one-loop vacuum polarization amplitude in
+QED:
+
+@example
+@{
+ symbol q("q"), l("l"), m("m"), ldotq("ldotq"), D("D");
+ varidx mu(symbol("mu"), D), nu(symbol("nu"), D);
+
+ scalar_products sp;
+ sp.add(l, l, pow(l, 2));
+ sp.add(l, q, ldotq);
+
+ ex e = dirac_gamma(mu) *
+ (dirac_slash(l, D) + dirac_slash(q, D) + m * dirac_ONE()) *
+ dirac_gamma(mu.toggle_variance()) *
+ (dirac_slash(l, D) + m * dirac_ONE());
+ e = dirac_trace(e).simplify_indexed(sp);
+ e = e.collect(lst(l, ldotq, m), true);
+ cout << e << endl;
+ // -> (8-4*D)*l^2+(8-4*D)*ldotq+4*D*m^2
+@}
+@end example
+
+The @code{canonicalize_clifford()} function reorders all gamma products that
+appear in an expression to a canonical (but not necessarily simple) form.
+You can use this to compare two expressions or for further simplifications:
+
+@example
+@{
+ varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
+ ex e = dirac_gamma(mu) * dirac_gamma(nu) + dirac_gamma(nu) * dirac_gamma(mu);
+ cout << e << endl;
+ // -> gamma~mu*gamma~nu+gamma~nu*gamma~mu
+
+ e = canonicalize_clifford(e);
+ cout << e << endl;
+ // -> 2*eta~mu~nu
+@}
+@end example
+
+
+@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 commute 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
+
+@cindex @code{color_d()}
+@cindex @code{color_f()}
+and 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}.
+
+@cindex @code{color_h()}
+There's an additional function
+
+@example
+ex color_h(const ex & a, const ex & b, const ex & c);
+@end example
+
+which returns the linear combination @samp{color_d(a, b, c)+I*color_f(a, b, c)}.
+
+The function @code{simplify_indexed()} performs some simplifications on
+expressions containing color objects:
+
+@example
+@{
+ ...
+ idx a(symbol("a"), 8), b(symbol("b"), 8), c(symbol("c"), 8),
+ k(symbol("k"), 8), l(symbol("l"), 8);
+
+ e = color_d(a, b, l) * color_f(a, b, k);
+ cout << e.simplify_indexed() << endl;
+ // -> 0
+
+ e = color_d(a, b, l) * color_d(a, b, k);
+ cout << e.simplify_indexed() << endl;
+ // -> 5/3*delta.k.l
+
+ e = color_f(l, a, b) * color_f(a, b, k);
+ cout << e.simplify_indexed() << endl;
+ // -> 3*delta.k.l
+
+ e = color_h(a, b, c) * color_h(a, b, c);
+ cout << e.simplify_indexed() << endl;
+ // -> -32/3
+
+ e = color_h(a, b, c) * color_T(b) * color_T(c);
+ cout << e.simplify_indexed() << endl;
+ // -> -2/3*T.a
+
+ e = color_h(a, b, c) * color_T(a) * color_T(b) * color_T(c);
+ cout << e.simplify_indexed() << endl;
+ // -> -8/9*ONE
+
+ e = color_T(k) * color_T(a) * color_T(b) * color_T(k);
+ cout << e.simplify_indexed() << endl;
+ // -> 1/4*delta.b.a*ONE-1/6*T.a*T.b
+ ...
+@end example
+
+@cindex @code{color_trace()}
+To calculate the trace of an expression containing color objects you use the
+function
+
+@example
+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:
+
+@example
+ ...
+ e = color_trace(4 * color_T(a) * color_T(b) * color_T(c));
+ cout << e << endl;
+ // -> -I*f.a.c.b+d.a.c.b
+@}
+@end example
+
+
+@node Methods and Functions, Information About Expressions, Non-commutative objects, Top
+@c node-name, next, previous, up
+@chapter Methods and Functions
+@cindex polynomial
+
+In this chapter the most important algorithms provided by GiNaC will be
+described. Some of them are implemented as functions on expressions,
+others are implemented as methods provided by expression objects. If
+they are methods, there exists a wrapper function around it, so you can
+alternatively call it in a functional way as shown in the simple
+example:
+
+@example
+ ...
+ cout << "As method: " << sin(1).evalf() << endl;
+ cout << "As function: " << evalf(sin(1)) << endl;
+ ...
+@end example
+
+@cindex @code{subs()}
+The general rule is that wherever methods accept one or more parameters
+(@var{arg1}, @var{arg2}, @dots{}) the order of arguments the function
+wrapper accepts is the same but preceded by the object to act on
+(@var{object}, @var{arg1}, @var{arg2}, @dots{}). This approach is the
+most natural one in an OO model but it may lead to confusion for MapleV
+users because where they would type @code{A:=x+1; subs(x=2,A);} GiNaC
+would require @code{A=x+1; subs(A,x==2);} (after proper declaration of
+@code{A} and @code{x}). On the other hand, since MapleV returns 3 on
+@code{A:=x^2+3; coeff(A,x,0);} (GiNaC: @code{A=pow(x,2)+3;
+coeff(A,x,0);}) it is clear that MapleV is not trying to be consistent
+here. Also, users of MuPAD will in most cases feel more comfortable
+with GiNaC's convention. All function wrappers are implemented
+as simple inline functions which just call the corresponding method and
+are only provided for users uncomfortable with OO who are dead set to
+avoid method invocations. Generally, nested function wrappers are much
+harder to read than a sequence of methods and should therefore be
+avoided if possible. On the other hand, not everything in GiNaC is a
+method on class @code{ex} and sometimes calling a function cannot be
+avoided.
+
+@menu
+* Information About Expressions::
+* Substituting Expressions::
+* Pattern Matching and Advanced Substitutions::
+* 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.
+* Input/Output:: Input and output of expressions.
+@end menu
+
+
+@node Information About Expressions, Substituting Expressions, Methods and Functions, Methods and Functions
+@c node-name, next, previous, up
+@section Getting information about expressions
+
+@subsection Checking expression types
+@cindex @code{is_ex_of_type()}
+@cindex @code{ex_to_numeric()}
+@cindex @code{ex_to_@dots{}}
+@cindex @code{Converting 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 (the first one is actually a macro):
+
+@example
+bool is_ex_of_type(const ex & e, TYPENAME t);
+bool ex::info(unsigned flag);
+unsigned ex::return_type(void) const;
+unsigned ex::return_type_tinfo(void) const;
+@end example
+
+When the test made by @code{is_ex_of_type()} returns true, it is safe to
+call one of the functions @code{ex_to_@dots{}}, where @code{@dots{}} 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_ex_of_type(e, numeric))
+ numeric n = ex_to_numeric(e);
+ @dots{}
+@}
+@end example
+
+@code{is_ex_of_type()} 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_ex_of_type(e1, numeric); // true
+ is_ex_of_type(e2, numeric); // false
+ is_ex_of_type(e1, add); // false
+ is_ex_of_type(e2, add); // true
+ is_ex_of_type(e1, mul); // false
+ is_ex_of_type(e2, mul); // false
+@}
+@end example
+
+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_ex_of_type(..., numeric)})
+@item @code{real}
+@tab @dots{}a real integer, rational or float (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_ex_of_type(..., 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_ex_of_type(..., symbol)})
+@item @code{list}
+@tab @dots{}a list (same as @code{is_ex_of_type(..., 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 commute, 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
+
+@example
+unsigned ex::nops();
+ex ex::op(unsigned 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.
+
+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
+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 (@xref{Relations}.)
+which is not evaluated until (explicitly or implicitely) 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.
+
+@strong{Warning:} You will also find an @code{ex::compare()} method in the
+GiNaC header files. This method is however only to be used internally by
+GiNaC to establish a canonical sort order for terms, and using it to compare
+expressions will give very surprising results.
+
+
+@node Substituting Expressions, Pattern Matching and Advanced Substitutions, Information About Expressions, 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);
+ex ex::subs(const lst & syms, const lst & repls);
+@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 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}.
+
+@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, Polynomial Arithmetic, Substituting Expressions, Methods and Functions
+@c node-name, next, previous, up
+@section Pattern matching and advanced substitutions
+
+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
+
+@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
+amgiguous 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 if 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
+
+@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+y),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+y),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{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
+
+
+@node Polynomial Arithmetic, Rational Expressions, Pattern Matching and Advanced Substitutions, Methods and Functions
+@c node-name, next, previous, up
+@section Polynomial arithmetic
+
+@subsection Expanding and collecting
+@cindex @code{expand()}
+@cindex @code{collect()}
+
+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 expample, 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();
+@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 easily guessable 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 in order
+for @code{collect()} to be able to find the coefficients properly.
+
+@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). 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
+#include <ginac/ginac.h>
+using namespace std;
+using namespace GiNaC;
+
+int main()
+@{
+ 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");
+ 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 symbol & x);
+ex rem(const ex & a, const ex & b, const symbol & 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 symbol & 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()}
+
+The methods
+
+@example
+ex ex::unit(const symbol & x);
+ex ex::content(const symbol & x);
+ex ex::primpart(const symbol & x);
+@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.
+
+
+@subsection GCD and LCM
+@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}.
+
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+int main()
+@{
+ symbol x("x"), y("y"), z("z");