+> 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");
+ 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()}
+
+The methods
+
+@example
+ex ex::unit(const ex & x);
+ex ex::content(const ex & x);
+ex ex::primpart(const ex & 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");
+ 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
+
+
+@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
+@code{1/cosh(x)}. Since the derivative of @code{cosh(x)} is
+@code{sinh(x)} and the derivative of @code{pow(x,-1)} is
+@code{-pow(x,-2)}, GiNaC can readily compute the composition. It turns
+out that the composition is the generating function for Euler Numbers,
+i.e. the so called @var{n}th Euler number is the coefficient of
+@code{x^n/n!} in the expansion of @code{1/cosh(x)}. We may use this
+identity to code a function that generates Euler numbers in just three
+lines:
+
+@cindex Euler numbers
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+ex EulerNumber(unsigned n)
+@{
+ symbol x;
+ const ex generator = pow(cosh(x),-1);
+ return generator.diff(x,n).subs(x==0);
+@}
+
+int main()
+@{
+ for (unsigned i=0; i<11; i+=2)
+ std::cout << EulerNumber(i) << std::endl;
+ return 0;
+@}
+@end example
+
+When you run it, it produces the sequence @code{1}, @code{-1}, @code{5},
+@code{-61}, @code{1385}, @code{-50521}. We increment the loop variable
+@code{i} by two since all odd Euler numbers vanish anyways.
+
+
+@node Series Expansion, Symmetrization, Symbolic Differentiation, Methods and Functions
+@c node-name, next, previous, up
+@section Series expansion
+@cindex @code{series()}
+@cindex Taylor expansion
+@cindex Laurent expansion
+@cindex @code{pseries} (class)
+@cindex @code{Order()}
+
+Expressions know how to expand themselves as a Taylor series or (more
+generally) a Laurent series. As in most conventional Computer Algebra
+Systems, no distinction is made between those two. There is a class of
+its own for storing such series (@code{class pseries}) and a built-in
+function (called @code{Order}) for storing the order term of the series.
+As a consequence, if you want to work with series, i.e. multiply two
+series, you need to call the method @code{ex::series} again to convert
+it to a series object with the usual structure (expansion plus order
+term). A sample application from special relativity could read:
+
+@example
+#include <ginac/ginac.h>
+using namespace std;
+using namespace GiNaC;
+
+int main()
+@{
+ symbol v("v"), c("c");
+
+ ex gamma = 1/sqrt(1 - pow(v/c,2));
+ ex mass_nonrel = gamma.series(v==0, 10);
+
+ cout << "the relativistic mass increase with v is " << endl
+ << mass_nonrel << endl;
+
+ cout << "the inverse square of this series is " << endl
+ << pow(mass_nonrel,-2).series(v==0, 10) << endl;
+@}
+@end example
+
+Only calling the series method makes the last output simplify to
+@math{1-v^2/c^2+O(v^10)}, without that call we would just have a long
+series raised to the power @math{-2}.
+
+@cindex Machin's formula
+As another instructive application, let us calculate the numerical
+value of Archimedes' constant
+@tex
+$\pi$
+@end tex
+(for which there already exists the built-in constant @code{Pi})
+using John Machin's amazing formula
+@tex
+$\pi=16$~atan~$\!\left(1 \over 5 \right)-4$~atan~$\!\left(1 \over 239 \right)$.
+@end tex
+@ifnottex
+@math{Pi==16*atan(1/5)-4*atan(1/239)}.
+@end ifnottex
+This equation (and similar ones) were used for over 200 years for
+computing digits of pi (see @cite{Pi Unleashed}). We may expand the
+arcus tangent around @code{0} and insert the fractions @code{1/5} and
+@code{1/239}. However, as we have seen, a series in GiNaC carries an
+order term with it and the question arises what the system is supposed
+to do when the fractions are plugged into that order term. The solution
+is to use the function @code{series_to_poly()} to simply strip the order
+term off:
+
+@example
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+ex machin_pi(int degr)
+@{
+ symbol x;
+ ex pi_expansion = series_to_poly(atan(x).series(x,degr));
+ ex pi_approx = 16*pi_expansion.subs(x==numeric(1,5))
+ -4*pi_expansion.subs(x==numeric(1,239));
+ return pi_approx;
+@}
+
+int main()
+@{
+ using std::cout; // just for fun, another way of...
+ using std::endl; // ...dealing with this namespace std.
+ ex pi_frac;
+ for (int i=2; i<12; i+=2) @{
+ pi_frac = machin_pi(i);
+ cout << i << ":\t" << pi_frac << endl
+ << "\t" << pi_frac.evalf() << endl;
+ @}
+ return 0;
+@}
+@end example
+
+Note how we just called @code{.series(x,degr)} instead of
+@code{.series(x==0,degr)}. This is a simple shortcut for @code{ex}'s
+method @code{series()}: if the first argument is a symbol the expression
+is expanded in that symbol around point @code{0}. When you run this
+program, it will type out:
+
+@example
+2: 3804/1195
+ 3.1832635983263598326
+4: 5359397032/1706489875
+ 3.1405970293260603143
+6: 38279241713339684/12184551018734375
+ 3.141621029325034425
+8: 76528487109180192540976/24359780855939418203125
+ 3.141591772182177295
+10: 327853873402258685803048818236/104359128170408663038552734375
+ 3.1415926824043995174
+@end example
+
+
+@node Symmetrization, Built-in Functions, Series Expansion, Methods and Functions
+@c node-name, next, previous, up
+@section Symmetrization
+@cindex @code{symmetrize()}
+@cindex @code{antisymmetrize()}
+@cindex @code{symmetrize_cyclic()}
+
+The three methods
+
+@example
+ex ex::symmetrize(const lst & l);
+ex ex::antisymmetrize(const lst & l);
+ex ex::symmetrize_cyclic(const lst & l);
+@end example
+
+symmetrize an expression by returning the sum over all symmetric,
+antisymmetric or cyclic permutations of the specified list of objects,
+weighted by the number of permutations.
+
+The three additional methods
+
+@example
+ex ex::symmetrize();
+ex ex::antisymmetrize();
+ex ex::symmetrize_cyclic();
+@end example
+
+symmetrize or antisymmetrize an expression over its free indices.
+
+Symmetrization is most useful with indexed expressions but can be used with
+almost any kind of object (anything that is @code{subs()}able):
+
+@example
+@{
+ idx i(symbol("i"), 3), j(symbol("j"), 3), k(symbol("k"), 3);
+ symbol A("A"), B("B"), a("a"), b("b"), c("c");
+
+ cout << indexed(A, i, j).symmetrize() << endl;
+ // -> 1/2*A.j.i+1/2*A.i.j
+ cout << indexed(A, i, j, k).antisymmetrize(lst(i, j)) << endl;
+ // -> -1/2*A.j.i.k+1/2*A.i.j.k
+ cout << lst(a, b, c).symmetrize_cyclic(lst(a, b, c)) << endl;
+ // -> 1/3*@{a,b,c@}+1/3*@{b,c,a@}+1/3*@{c,a,b@}
+@}
+@end example
+
+@node Built-in Functions, Complex Conjugation, Symmetrization, Methods and Functions
+@c node-name, next, previous, up
+@section Predefined mathematical functions
+@c
+@subsection Overview
+
+GiNaC contains the following predefined mathematical functions:
+
+@cartouche
+@multitable @columnfractions .30 .70
+@item @strong{Name} @tab @strong{Function}
+@item @code{abs(x)}
+@tab absolute value
+@cindex @code{abs()}
+@item @code{csgn(x)}
+@tab complex sign
+@cindex @code{conjugate()}
+@item @code{conjugate(x)}
+@tab complex conjugation
+@cindex @code{csgn()}
+@item @code{sqrt(x)}
+@tab square root (not a GiNaC function, rather an alias for @code{pow(x, numeric(1, 2))})
+@cindex @code{sqrt()}
+@item @code{sin(x)}
+@tab sine
+@cindex @code{sin()}
+@item @code{cos(x)}
+@tab cosine
+@cindex @code{cos()}
+@item @code{tan(x)}
+@tab tangent
+@cindex @code{tan()}
+@item @code{asin(x)}
+@tab inverse sine
+@cindex @code{asin()}
+@item @code{acos(x)}
+@tab inverse cosine
+@cindex @code{acos()}
+@item @code{atan(x)}
+@tab inverse tangent
+@cindex @code{atan()}
+@item @code{atan2(y, x)}
+@tab inverse tangent with two arguments
+@item @code{sinh(x)}
+@tab hyperbolic sine
+@cindex @code{sinh()}
+@item @code{cosh(x)}
+@tab hyperbolic cosine
+@cindex @code{cosh()}
+@item @code{tanh(x)}
+@tab hyperbolic tangent
+@cindex @code{tanh()}
+@item @code{asinh(x)}
+@tab inverse hyperbolic sine
+@cindex @code{asinh()}
+@item @code{acosh(x)}
+@tab inverse hyperbolic cosine
+@cindex @code{acosh()}
+@item @code{atanh(x)}
+@tab inverse hyperbolic tangent
+@cindex @code{atanh()}
+@item @code{exp(x)}
+@tab exponential function
+@cindex @code{exp()}
+@item @code{log(x)}
+@tab natural logarithm
+@cindex @code{log()}
+@item @code{Li2(x)}
+@tab dilogarithm
+@cindex @code{Li2()}
+@item @code{Li(m, x)}
+@tab classical polylogarithm as well as multiple polylogarithm
+@cindex @code{Li()}
+@item @code{S(n, p, x)}
+@tab Nielsen's generalized polylogarithm
+@cindex @code{S()}
+@item @code{H(m, x)}
+@tab harmonic polylogarithm
+@cindex @code{H()}
+@item @code{zeta(m)}
+@tab Riemann's zeta function as well as multiple zeta value
+@cindex @code{zeta()}
+@item @code{zeta(m, s)}
+@tab alternating Euler sum
+@cindex @code{zeta()}
+@item @code{zetaderiv(n, x)}
+@tab derivatives of Riemann's zeta function
+@item @code{tgamma(x)}
+@tab gamma function
+@cindex @code{tgamma()}
+@cindex gamma function
+@item @code{lgamma(x)}
+@tab logarithm of gamma function
+@cindex @code{lgamma()}
+@item @code{beta(x, y)}
+@tab beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
+@cindex @code{beta()}
+@item @code{psi(x)}
+@tab psi (digamma) function
+@cindex @code{psi()}
+@item @code{psi(n, x)}
+@tab derivatives of psi function (polygamma functions)
+@item @code{factorial(n)}
+@tab factorial function
+@cindex @code{factorial()}
+@item @code{binomial(n, m)}
+@tab binomial coefficients
+@cindex @code{binomial()}
+@item @code{Order(x)}
+@tab order term function in truncated power series
+@cindex @code{Order()}
+@end multitable
+@end cartouche
+
+@cindex branch cut
+For functions that have a branch cut in the complex plane GiNaC follows
+the conventions for C++ as defined in the ANSI standard as far as
+possible. In particular: the natural logarithm (@code{log}) and the
+square root (@code{sqrt}) both have their branch cuts running along the
+negative real axis where the points on the axis itself belong to the
+upper part (i.e. continuous with quadrant II). The inverse
+trigonometric and hyperbolic functions are not defined for complex
+arguments by the C++ standard, however. In GiNaC we follow the
+conventions used by CLN, which in turn follow the carefully designed
+definitions in the Common Lisp standard. It should be noted that this
+convention is identical to the one used by the C99 standard and by most
+serious CAS. It is to be expected that future revisions of the C++
+standard incorporate these functions in the complex domain in a manner
+compatible with C99.
+
+@subsection Multiple polylogarithms
+
+@cindex polylogarithm
+@cindex Nielsen's generalized polylogarithm
+@cindex harmonic polylogarithm
+@cindex multiple zeta value
+@cindex alternating Euler sum
+@cindex multiple polylogarithm
+
+The multiple polylogarithm is the most generic member of a family of functions,
+to which others like the harmonic polylogarithm, Nielsen's generalized
+polylogarithm and the multiple zeta value belong.
+Everyone of these functions can also be written as a multiple polylogarithm with specific
+parameters. This whole family of functions is therefore often referred to simply as
+multiple polylogarithms, containing @code{Li}, @code{H}, @code{S} and @code{zeta}.
+
+To facilitate the discussion of these functions we distinguish between indices and
+arguments as parameters. In the table above indices are printed as @code{m}, @code{s},
+@code{n} or @code{p}, whereas arguments are printed as @code{x}.
+
+To define a @code{Li}, @code{H} or @code{zeta} with a depth greater than one, you have to
+pass a GiNaC @code{lst} for the indices @code{m} and @code{s}, and in the case of @code{Li}
+for the argument @code{x} as well.
+Note that @code{Li} and @code{zeta} are polymorphic in this respect. They can stand in for
+the classical polylogarithm and Riemann's zeta function (if depth is one), as well as for
+the multiple polylogarithm and the multiple zeta value, respectively. Note also, that
+GiNaC doesn't check whether the @code{lst}s for two parameters do have the same length.
+It is up to the user to ensure this, otherwise evaluating will result in undefined behavior.
+
+The functions print in LaTeX format as
+@tex
+${\rm Li\;\!}_{m_1,m_2,\ldots,m_k}(x_1,x_2,\ldots,x_k)$,
+@end tex
+@tex
+${\rm S}_{n,p}(x)$,
+@end tex
+@tex
+${\rm H\;\!}_{m_1,m_2,\ldots,m_k}(x)$ and
+@end tex
+@tex
+$\zeta(m_1,m_2,\ldots,m_k)$.
+@end tex
+If @code{zeta} is an alternating zeta sum, i.e. @code{zeta(m,s)}, the indices with negative sign
+are printed with a line above, e.g.
+@tex
+$\zeta(5,\overline{2})$.
+@end tex
+The order of indices and arguments in the GiNaC @code{lst}s and in the output is the same.
+
+Definitions and analytical as well as numerical properties of multiple polylogarithms
+are too numerous to be covered here. Instead, the user is referred to the publications listed at the
+end of this section. The implementation in GiNaC adheres to the definitions and conventions therein,
+except for a few differences which will be explicitly stated in the following.
+
+One difference is about the order of the indices and arguments. For GiNaC we adopt the convention
+that the indices and arguments are understood to be in the same order as in which they appear in
+the series representation. This means
+@tex
+${\rm Li\;\!}_{m_1,m_2,m_3}(x,1,1) = {\rm H\;\!}_{m_1,m_2,m_3}(x)$ and
+@end tex
+@tex
+${\rm Li\;\!}_{2,1}(1,1) = \zeta(2,1) = \zeta(3)$, but
+@end tex
+@tex
+$\zeta(1,2)$ evaluates to infinity.
+@end tex
+So in comparison to the referenced publications the order of indices and arguments for @code{Li}
+is reversed.
+
+The functions only evaluate if the indices are integers greater than zero, except for the indices
+@code{s} in @code{zeta} and @code{m} in @code{H}. Since @code{s} will be interpreted as the sequence
+of signs for the corresponding indices @code{m}, it must contain 1 or -1, e.g.
+@code{zeta(lst(3,4), lst(-1,1))} means
+@tex
+$\zeta(\overline{3},4)$.
+@end tex
+The definition of @code{H} allows indices to be 0, 1 or -1 (in expanded notation) or equally to
+be any integer (in compact notation). With GiNaC expanded and compact notation can be mixed,
+e.g. @code{lst(0,0,-1,0,1,0,0)}, @code{lst(0,0,-1,2,0,0)} and @code{lst(-3,2,0,0)} are equivalent as
+indices. The anonymous evaluator @code{eval()} tries to reduce the functions, if possible, to
+the least-generic multiple polylogarithm. If all arguments are unit, it returns @code{zeta}.
+Arguments equal to zero get considered, too. Riemann's zeta function @code{zeta} (with depth one)
+evaluates also for negative integers and positive even integers. For example:
+
+@example
+> Li(@{3,1@},@{x,1@});
+S(2,2,x)
+> H(@{-3,2@},1);
+-zeta(@{3,2@},@{-1,-1@})
+> S(3,1,1);
+1/90*Pi^4
+@end example
+
+It is easy to tell for a given function into which other function it can be rewritten, may
+it be a less-generic or a more-generic one, except for harmonic polylogarithms @code{H}
+with negative indices or trailing zeros (the example above gives a hint). Signs can
+quickly be messed up, for example. Therefore GiNaC offers a C++ function
+@code{convert_H_to_Li()} to deal with the upgrade of a @code{H} to a multiple polylogarithm
+@code{Li} (@code{eval()} already cares for the possible downgrade):
+
+@example
+> convert_H_to_Li(@{0,-2,-1,3@},x);
+Li(@{3,1,3@},@{-x,1,-1@})
+> convert_H_to_Li(@{2,-1,0@},x);
+-Li(@{2,1@},@{x,-1@})*log(x)+2*Li(@{3,1@},@{x,-1@})+Li(@{2,2@},@{x,-1@})
+@end example
+
+Every function apart from the multiple polylogarithm @code{Li} can be numerically evaluated for
+arbitrary real or complex arguments. @code{Li} only evaluates if for all arguments
+@tex
+$x_i$ the condition
+@end tex
+@tex
+$x_1x_2\cdots x_i < 1$ holds.
+@end tex
+
+@example
+> Digits=100;
+100
+> evalf(zeta(@{3,1,3,1@}));
+0.005229569563530960100930652283899231589890420784634635522547448972148869544...
+@end example
+
+Note that the convention for arguments on the branch cut in GiNaC as stated above is
+different from the one Remiddi and Vermaseren have chosen for the harmonic polylogarithm.
+
+If a function evaluates to infinity, no exceptions are raised, but the function is returned
+unevaluated, e.g.
+@tex
+$\zeta(1)$.
+@end tex
+In long expressions this helps a lot with debugging, because you can easily spot
+the divergencies. But on the other hand, you have to make sure for yourself, that no illegal
+cancellations of divergencies happen.
+
+Useful publications:
+
+@cite{Nested Sums, Expansion of Transcendental Functions and Multi-Scale Multi-Loop Integrals},
+S.Moch, P.Uwer, S.Weinzierl, hep-ph/0110083
+
+@cite{Harmonic Polylogarithms},
+E.Remiddi, J.A.M.Vermaseren, Int.J.Mod.Phys. A15 (2000), pp. 725-754
+
+@cite{Special Values of Multiple Polylogarithms},
+J.Borwein, D.Bradley, D.Broadhurst, P.Lisonek, Trans.Amer.Math.Soc. 353/3 (2001), pp. 907-941
+
+@node Complex Conjugation, Solving Linear Systems of Equations, Built-in Functions, Methods and Functions
+@c node-name, next, previous, up
+@section Complex Conjugation
+@c
+@cindex @code{conjugate()}
+
+The method
+
+@example
+ex ex::conjugate();
+@end example
+
+returns the complex conjugate of the expression. For all built-in functions and objects the
+conjugation gives the expected results:
+
+@example
+@{
+ varidx a(symbol("a"), 4), b(symbol("b"), 4);
+ symbol x("x");
+ realsymbol y("y");
+
+ cout << (3*I*x*y + sin(2*Pi*I*y)).conjugate() << endl;
+ // -> -3*I*conjugate(x)*y+sin(-2*I*Pi*y)
+ cout << (dirac_gamma(a)*dirac_gamma(b)*dirac_gamma5()).conjugate() << endl;
+ // -> -gamma5*gamma~b*gamma~a
+@}
+@end example
+
+For symbols in the complex domain the conjugation can not be evaluated and the GiNaC function
+@code{conjugate} is returned. GiNaC functions conjugate by applying the conjugation to their
+arguments. This is the default strategy. If you want to define your own functions and want to
+change this behavior, you have to supply a specialized conjugation method for your function
+(see @ref{Symbolic functions} and the GiNaC source-code for @code{abs} as an example).
+
+@node Solving Linear Systems of Equations, Input/Output, Complex Conjugation, Methods and Functions
+@c node-name, next, previous, up
+@section Solving Linear Systems of Equations
+@cindex @code{lsolve()}
+
+The function @code{lsolve()} provides a convenient wrapper around some
+matrix operations that comes in handy when a system of linear equations
+needs to be solved:
+
+@example
+ex lsolve(const ex &eqns, const ex &symbols, unsigned options=solve_algo::automatic);
+@end example
+
+Here, @code{eqns} is a @code{lst} of equalities (i.e. class
+@code{relational}) while @code{symbols} is a @code{lst} of
+indeterminates. (@xref{The Class Hierarchy}, for an exposition of class
+@code{lst}).
+
+It returns the @code{lst} of solutions as an expression. As an example,
+let us solve the two equations @code{a*x+b*y==3} and @code{x-y==b}:
+
+@example
+@{
+ symbol a("a"), b("b"), x("x"), y("y");
+ lst eqns, vars;
+ eqns = a*x+b*y==3, x-y==b;
+ vars = x, y;
+ cout << lsolve(eqns, vars) << endl;
+ // -> @{x==(3+b^2)/(b+a),y==(3-b*a)/(b+a)@}
+@end example
+
+When the linear equations @code{eqns} are underdetermined, the solution
+will contain one or more tautological entries like @code{x==x},
+depending on the rank of the system. When they are overdetermined, the
+solution will be an empty @code{lst}. Note the third optional parameter
+to @code{lsolve()}: it accepts the same parameters as
+@code{matrix::solve()}. This is because @code{lsolve} is just a wrapper
+around that method.
+
+
+@node Input/Output, Extending GiNaC, Solving Linear Systems of Equations, Methods and Functions
+@c node-name, next, previous, up
+@section Input and output of expressions
+@cindex I/O
+
+@subsection Expression output
+@cindex printing
+@cindex output of expressions
+
+Expressions can simply be written to any stream:
+
+@example
+@{
+ symbol x("x");
+ ex e = 4.5*I+pow(x,2)*3/2;
+ cout << e << endl; // prints '4.5*I+3/2*x^2'
+ // ...
+@end example
+
+The default output format is identical to the @command{ginsh} input syntax and
+to that used by most computer algebra systems, but not directly pastable
+into a GiNaC C++ program (note that in the above example, @code{pow(x,2)}
+is printed as @samp{x^2}).
+
+It is possible to print expressions in a number of different formats with
+a set of stream manipulators;
+
+@example
+std::ostream & dflt(std::ostream & os);
+std::ostream & latex(std::ostream & os);
+std::ostream & tree(std::ostream & os);
+std::ostream & csrc(std::ostream & os);
+std::ostream & csrc_float(std::ostream & os);
+std::ostream & csrc_double(std::ostream & os);
+std::ostream & csrc_cl_N(std::ostream & os);
+std::ostream & index_dimensions(std::ostream & os);
+std::ostream & no_index_dimensions(std::ostream & os);
+@end example
+
+The @code{tree}, @code{latex} and @code{csrc} formats are also available in
+@command{ginsh} via the @code{print()}, @code{print_latex()} and
+@code{print_csrc()} functions, respectively.
+
+@cindex @code{dflt}
+All manipulators affect the stream state permanently. To reset the output
+format to the default, use the @code{dflt} manipulator:
+
+@example
+ // ...
+ cout << latex; // all output to cout will be in LaTeX format from now on
+ cout << e << endl; // prints '4.5 i+\frac@{3@}@{2@} x^@{2@}'
+ cout << sin(x/2) << endl; // prints '\sin(\frac@{1@}@{2@} x)'
+ cout << dflt; // revert to default output format
+ cout << e << endl; // prints '4.5*I+3/2*x^2'
+ // ...
+@end example
+
+If you don't want to affect the format of the stream you're working with,
+you can output to a temporary @code{ostringstream} like this:
+
+@example
+ // ...
+ ostringstream s;
+ s << latex << e; // format of cout remains unchanged
+ cout << s.str() << endl; // prints '4.5 i+\frac@{3@}@{2@} x^@{2@}'
+ // ...
+@end example
+
+@cindex @code{csrc}
+@cindex @code{csrc_float}
+@cindex @code{csrc_double}
+@cindex @code{csrc_cl_N}
+The @code{csrc} (an alias for @code{csrc_double}), @code{csrc_float},
+@code{csrc_double} and @code{csrc_cl_N} manipulators set the output to a
+format that can be directly used in a C or C++ program. The three possible
+formats select the data types used for numbers (@code{csrc_cl_N} uses the
+classes provided by the CLN library):
+
+@example
+ // ...
+ cout << "f = " << csrc_float << e << ";\n";
+ cout << "d = " << csrc_double << e << ";\n";
+ cout << "n = " << csrc_cl_N << e << ";\n";
+ // ...
+@end example
+
+The above example will produce (note the @code{x^2} being converted to
+@code{x*x}):
+
+@example
+f = (3.0/2.0)*(x*x)+std::complex<float>(0.0,4.5000000e+00);
+d = (3.0/2.0)*(x*x)+std::complex<double>(0.0,4.5000000000000000e+00);
+n = cln::cl_RA("3/2")*(x*x)+cln::complex(cln::cl_I("0"),cln::cl_F("4.5_17"));
+@end example
+
+@cindex @code{tree}
+The @code{tree} manipulator allows dumping the internal structure of an
+expression for debugging purposes:
+
+@example
+ // ...
+ cout << tree << e;
+@}
+@end example
+
+produces
+
+@example
+add, hash=0x0, flags=0x3, nops=2
+ power, hash=0x0, flags=0x3, nops=2
+ x (symbol), serial=0, hash=0xc8d5bcdd, flags=0xf
+ 2 (numeric), hash=0x6526b0fa, flags=0xf
+ 3/2 (numeric), hash=0xf9828fbd, flags=0xf
+ -----
+ overall_coeff
+ 4.5L0i (numeric), hash=0xa40a97e0, flags=0xf
+ =====
+@end example
+
+@cindex @code{latex}
+The @code{latex} output format is for LaTeX parsing in mathematical mode.
+It is rather similar to the default format but provides some braces needed
+by LaTeX for delimiting boxes and also converts some common objects to
+conventional LaTeX names. It is possible to give symbols a special name for
+LaTeX output by supplying it as a second argument to the @code{symbol}
+constructor.
+
+For example, the code snippet
+
+@example
+@{
+ symbol x("x", "\\circ");
+ ex e = lgamma(x).series(x==0,3);
+ cout << latex << e << endl;
+@}
+@end example
+
+will print
+
+@example
+ @{(-\ln(\circ))@}+@{(-\gamma_E)@} \circ+@{(\frac@{1@}@{12@} \pi^@{2@})@} \circ^@{2@}+\mathcal@{O@}(\circ^@{3@})
+@end example
+
+@cindex @code{index_dimensions}
+@cindex @code{no_index_dimensions}
+Index dimensions are normally hidden in the output. To make them visible, use
+the @code{index_dimensions} manipulator. The dimensions will be written in
+square brackets behind each index value in the default and LaTeX output
+formats:
+
+@example
+@{
+ symbol x("x"), y("y");
+ varidx mu(symbol("mu"), 4), nu(symbol("nu"), 4);
+ ex e = indexed(x, mu) * indexed(y, nu);
+
+ cout << e << endl;
+ // prints 'x~mu*y~nu'
+ cout << index_dimensions << e << endl;
+ // prints 'x~mu[4]*y~nu[4]'
+ cout << no_index_dimensions << e << endl;
+ // prints 'x~mu*y~nu'
+@}
+@end example
+
+
+@cindex Tree traversal
+If you need any fancy special output format, e.g. for interfacing GiNaC
+with other algebra systems or for producing code for different
+programming languages, you can always traverse the expression tree yourself:
+
+@example
+static void my_print(const ex & e)
+@{
+ if (is_a<function>(e))
+ cout << ex_to<function>(e).get_name();
+ else
+ cout << ex_to<basic>(e).class_name();
+ cout << "(";
+ size_t n = e.nops();
+ if (n)
+ for (size_t i=0; i<n; i++) @{
+ my_print(e.op(i));
+ if (i != n-1)
+ cout << ",";
+ @}
+ else
+ cout << e;
+ cout << ")";
+@}
+
+int main()
+@{
+ my_print(pow(3, x) - 2 * sin(y / Pi)); cout << endl;
+ return 0;
+@}
+@end example
+
+This will produce
+
+@example
+add(power(numeric(3),symbol(x)),mul(sin(mul(power(constant(Pi),numeric(-1)),
+symbol(y))),numeric(-2)))
+@end example
+
+If you need an output format that makes it possible to accurately
+reconstruct an expression by feeding the output to a suitable parser or
+object factory, you should consider storing the expression in an
+@code{archive} object and reading the object properties from there.
+See the section on archiving for more information.
+
+
+@subsection Expression input
+@cindex input of expressions
+
+GiNaC provides no way to directly read an expression from a stream because
+you will usually want the user to be able to enter something like @samp{2*x+sin(y)}
+and have the @samp{x} and @samp{y} correspond to the symbols @code{x} and
+@code{y} you defined in your program and there is no way to specify the
+desired symbols to the @code{>>} stream input operator.
+
+Instead, GiNaC lets you construct an expression from a string, specifying the
+list of symbols to be used:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex e("2*x+sin(y)", lst(x, y));
+@}
+@end example
+
+The input syntax is the same as that used by @command{ginsh} and the stream
+output operator @code{<<}. The symbols in the string are matched by name to
+the symbols in the list and if GiNaC encounters a symbol not specified in
+the list it will throw an exception.
+
+With this constructor, it's also easy to implement interactive GiNaC programs:
+
+@example
+#include <iostream>
+#include <string>
+#include <stdexcept>
+#include <ginac/ginac.h>
+using namespace std;
+using namespace GiNaC;
+
+int main()
+@{
+ symbol x("x");
+ string s;
+
+ cout << "Enter an expression containing 'x': ";
+ getline(cin, s);
+
+ try @{
+ ex e(s, lst(x));
+ cout << "The derivative of " << e << " with respect to x is ";
+ cout << e.diff(x) << ".\n";
+ @} catch (exception &p) @{
+ cerr << p.what() << endl;
+ @}
+@}
+@end example
+
+
+@subsection Archiving
+@cindex @code{archive} (class)
+@cindex archiving
+
+GiNaC allows creating @dfn{archives} of expressions which can be stored
+to or retrieved from files. To create an archive, you declare an object
+of class @code{archive} and archive expressions in it, giving each
+expression a unique name:
+
+@example
+#include <fstream>
+using namespace std;
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+int main()
+@{
+ symbol x("x"), y("y"), z("z");
+
+ ex foo = sin(x + 2*y) + 3*z + 41;
+ ex bar = foo + 1;
+
+ archive a;
+ a.archive_ex(foo, "foo");
+ a.archive_ex(bar, "the second one");
+ // ...
+@end example
+
+The archive can then be written to a file:
+
+@example
+ // ...
+ ofstream out("foobar.gar");
+ out << a;
+ out.close();
+ // ...
+@end example
+
+The file @file{foobar.gar} contains all information that is needed to
+reconstruct the expressions @code{foo} and @code{bar}.
+
+@cindex @command{viewgar}
+The tool @command{viewgar} that comes with GiNaC can be used to view
+the contents of GiNaC archive files:
+
+@example
+$ viewgar foobar.gar
+foo = 41+sin(x+2*y)+3*z
+the second one = 42+sin(x+2*y)+3*z
+@end example
+
+The point of writing archive files is of course that they can later be
+read in again:
+
+@example
+ // ...
+ archive a2;
+ ifstream in("foobar.gar");
+ in >> a2;
+ // ...
+@end example
+
+And the stored expressions can be retrieved by their name:
+
+@example
+ // ...
+ lst syms;
+ syms = x, y;
+
+ ex ex1 = a2.unarchive_ex(syms, "foo");
+ ex ex2 = a2.unarchive_ex(syms, "the second one");
+
+ cout << ex1 << endl; // prints "41+sin(x+2*y)+3*z"
+ cout << ex2 << endl; // prints "42+sin(x+2*y)+3*z"
+ cout << ex1.subs(x == 2) << endl; // prints "41+sin(2+2*y)+3*z"
+@}
+@end example
+
+Note that you have to supply a list of the symbols which are to be inserted
+in the expressions. Symbols in archives are stored by their name only and
+if you don't specify which symbols you have, unarchiving the expression will
+create new symbols with that name. E.g. if you hadn't included @code{x} in
+the @code{syms} list above, the @code{ex1.subs(x == 2)} statement would
+have had no effect because the @code{x} in @code{ex1} would have been a
+different symbol than the @code{x} which was defined at the beginning of
+the program, although both would appear as @samp{x} when printed.
+
+You can also use the information stored in an @code{archive} object to
+output expressions in a format suitable for exact reconstruction. The
+@code{archive} and @code{archive_node} classes have a couple of member
+functions that let you access the stored properties:
+
+@example
+static void my_print2(const archive_node & n)
+@{
+ string class_name;
+ n.find_string("class", class_name);
+ cout << class_name << "(";
+
+ archive_node::propinfovector p;
+ n.get_properties(p);
+
+ size_t num = p.size();
+ for (size_t i=0; i<num; i++) @{
+ const string &name = p[i].name;
+ if (name == "class")
+ continue;
+ cout << name << "=";
+
+ unsigned count = p[i].count;
+ if (count > 1)
+ cout << "@{";
+
+ for (unsigned j=0; j<count; j++) @{
+ switch (p[i].type) @{
+ case archive_node::PTYPE_BOOL: @{
+ bool x;
+ n.find_bool(name, x, j);
+ cout << (x ? "true" : "false");
+ break;
+ @}
+ case archive_node::PTYPE_UNSIGNED: @{
+ unsigned x;
+ n.find_unsigned(name, x, j);
+ cout << x;
+ break;
+ @}
+ case archive_node::PTYPE_STRING: @{
+ string x;
+ n.find_string(name, x, j);
+ cout << '\"' << x << '\"';
+ break;
+ @}
+ case archive_node::PTYPE_NODE: @{
+ const archive_node &x = n.find_ex_node(name, j);
+ my_print2(x);
+ break;
+ @}
+ @}
+
+ if (j != count-1)
+ cout << ",";
+ @}
+
+ if (count > 1)
+ cout << "@}";
+
+ if (i != num-1)
+ cout << ",";
+ @}
+
+ cout << ")";
+@}
+
+int main()
+@{
+ ex e = pow(2, x) - y;
+ archive ar(e, "e");
+ my_print2(ar.get_top_node(0)); cout << endl;
+ return 0;
+@}
+@end example
+
+This will produce:
+
+@example
+add(rest=@{power(basis=numeric(number="2"),exponent=symbol(name="x")),
+symbol(name="y")@},coeff=@{numeric(number="1"),numeric(number="-1")@},
+overall_coeff=numeric(number="0"))
+@end example
+
+Be warned, however, that the set of properties and their meaning for each
+class may change between GiNaC versions.
+
+
+@node Extending GiNaC, What does not belong into GiNaC, Input/Output, Top
+@c node-name, next, previous, up
+@chapter Extending GiNaC
+
+By reading so far you should have gotten a fairly good understanding of
+GiNaC's design patterns. From here on you should start reading the
+sources. All we can do now is issue some recommendations how to tackle
+GiNaC's many loose ends in order to fulfill everybody's dreams. If you
+develop some useful extension please don't hesitate to contact the GiNaC
+authors---they will happily incorporate them into future versions.
+
+@menu
+* What does not belong into GiNaC:: What to avoid.
+* Symbolic functions:: Implementing symbolic functions.
+* Printing:: Adding new output formats.
+* Structures:: Defining new algebraic classes (the easy way).
+* Adding classes:: Defining new algebraic classes (the hard way).
+@end menu
+
+
+@node What does not belong into GiNaC, Symbolic functions, Extending GiNaC, Extending GiNaC
+@c node-name, next, previous, up
+@section What doesn't belong into GiNaC
+
+@cindex @command{ginsh}
+First of all, GiNaC's name must be read literally. It is designed to be
+a library for use within C++. The tiny @command{ginsh} accompanying
+GiNaC makes this even more clear: it doesn't even attempt to provide a
+language. There are no loops or conditional expressions in
+@command{ginsh}, it is merely a window into the library for the
+programmer to test stuff (or to show off). Still, the design of a
+complete CAS with a language of its own, graphical capabilities and all
+this on top of GiNaC is possible and is without doubt a nice project for
+the future.
+
+There are many built-in functions in GiNaC that do not know how to
+evaluate themselves numerically to a precision declared at runtime
+(using @code{Digits}). Some may be evaluated at certain points, but not
+generally. This ought to be fixed. However, doing numerical
+computations with GiNaC's quite abstract classes is doomed to be
+inefficient. For this purpose, the underlying foundation classes
+provided by CLN are much better suited.
+
+
+@node Symbolic functions, Printing, What does not belong into GiNaC, Extending GiNaC
+@c node-name, next, previous, up
+@section Symbolic functions
+
+The easiest and most instructive way to start extending GiNaC is probably to
+create your own symbolic functions. These are implemented with the help of
+two preprocessor macros:
+
+@cindex @code{DECLARE_FUNCTION}
+@cindex @code{REGISTER_FUNCTION}
+@example
+DECLARE_FUNCTION_<n>P(<name>)
+REGISTER_FUNCTION(<name>, <options>)
+@end example
+
+The @code{DECLARE_FUNCTION} macro will usually appear in a header file. It
+declares a C++ function with the given @samp{name} that takes exactly @samp{n}
+parameters of type @code{ex} and returns a newly constructed GiNaC
+@code{function} object that represents your function.
+
+The @code{REGISTER_FUNCTION} macro implements the function. It must be passed
+the same @samp{name} as the respective @code{DECLARE_FUNCTION} macro, and a
+set of options that associate the symbolic function with C++ functions you
+provide to implement the various methods such as evaluation, derivative,
+series expansion etc. They also describe additional attributes the function
+might have, such as symmetry and commutation properties, and a name for
+LaTeX output. Multiple options are separated by the member access operator
+@samp{.} and can be given in an arbitrary order.
+
+(By the way: in case you are worrying about all the macros above we can
+assure you that functions are GiNaC's most macro-intense classes. We have
+done our best to avoid macros where we can.)
+
+@subsection A minimal example
+
+Here is an example for the implementation of a function with two arguments
+that is not further evaluated:
+
+@example
+DECLARE_FUNCTION_2P(myfcn)
+
+REGISTER_FUNCTION(myfcn, dummy())
+@end example
+
+Any code that has seen the @code{DECLARE_FUNCTION} line can use @code{myfcn()}
+in algebraic expressions:
+
+@example
+@{
+ ...
+ symbol x("x");
+ ex e = 2*myfcn(42, 1+3*x) - x;
+ cout << e << endl;
+ // prints '2*myfcn(42,1+3*x)-x'
+ ...
+@}
+@end example
+
+The @code{dummy()} option in the @code{REGISTER_FUNCTION} line signifies
+"no options". A function with no options specified merely acts as a kind of
+container for its arguments. It is a pure "dummy" function with no associated
+logic (which is, however, sometimes perfectly sufficient).
+
+Let's now have a look at the implementation of GiNaC's cosine function for an
+example of how to make an "intelligent" function.
+
+@subsection The cosine function
+
+The GiNaC header file @file{inifcns.h} contains the line
+
+@example
+DECLARE_FUNCTION_1P(cos)
+@end example
+
+which declares to all programs using GiNaC that there is a function @samp{cos}
+that takes one @code{ex} as an argument. This is all they need to know to use
+this function in expressions.
+
+The implementation of the cosine function is in @file{inifcns_trans.cpp}. Here
+is its @code{REGISTER_FUNCTION} line:
+
+@example
+REGISTER_FUNCTION(cos, eval_func(cos_eval).
+ evalf_func(cos_evalf).
+ derivative_func(cos_deriv).
+ latex_name("\\cos"));
+@end example
+
+There are four options defined for the cosine function. One of them
+(@code{latex_name}) gives the function a proper name for LaTeX output; the
+other three indicate the C++ functions in which the "brains" of the cosine
+function are defined.
+
+@cindex @code{hold()}
+@cindex evaluation
+The @code{eval_func()} option specifies the C++ function that implements
+the @code{eval()} method, GiNaC's anonymous evaluator. This function takes
+the same number of arguments as the associated symbolic function (one in this
+case) and returns the (possibly transformed or in some way simplified)
+symbolically evaluated function (@xref{Automatic evaluation}, for a description
+of the automatic evaluation process). If no (further) evaluation is to take
+place, the @code{eval_func()} function must return the original function
+with @code{.hold()}, to avoid a potential infinite recursion. If your
+symbolic functions produce a segmentation fault or stack overflow when
+using them in expressions, you are probably missing a @code{.hold()}
+somewhere.
+
+The @code{eval_func()} function for the cosine looks something like this
+(actually, it doesn't look like this at all, but it should give you an idea
+what is going on):
+
+@example
+static ex cos_eval(const ex & x)
+@{
+ if ("x is a multiple of 2*Pi")
+ return 1;
+ else if ("x is a multiple of Pi")
+ return -1;
+ else if ("x is a multiple of Pi/2")
+ return 0;
+ // more rules...
+
+ else if ("x has the form 'acos(y)'")
+ return y;
+ else if ("x has the form 'asin(y)'")
+ return sqrt(1-y^2);
+ // more rules...
+
+ else
+ return cos(x).hold();
+@}
+@end example
+
+This function is called every time the cosine is used in a symbolic expression:
+
+@example
+@{
+ ...
+ e = cos(Pi);
+ // this calls cos_eval(Pi), and inserts its return value into
+ // the actual expression
+ cout << e << endl;
+ // prints '-1'
+ ...
+@}
+@end example
+
+In this way, @code{cos(4*Pi)} automatically becomes @math{1},
+@code{cos(asin(a+b))} becomes @code{sqrt(1-(a+b)^2)}, etc. If no reasonable
+symbolic transformation can be done, the unmodified function is returned
+with @code{.hold()}.
+
+GiNaC doesn't automatically transform @code{cos(2)} to @samp{-0.416146...}.
+The user has to call @code{evalf()} for that. This is implemented in a
+different function:
+
+@example
+static ex cos_evalf(const ex & x)
+@{
+ if (is_a<numeric>(x))
+ return cos(ex_to<numeric>(x));
+ else
+ return cos(x).hold();
+@}
+@end example
+
+Since we are lazy we defer the problem of numeric evaluation to somebody else,
+in this case the @code{cos()} function for @code{numeric} objects, which in
+turn hands it over to the @code{cos()} function in CLN. The @code{.hold()}
+isn't really needed here, but reminds us that the corresponding @code{eval()}
+function would require it in this place.
+
+Differentiation will surely turn up and so we need to tell @code{cos}
+what its first derivative is (higher derivatives, @code{.diff(x,3)} for
+instance, are then handled automatically by @code{basic::diff} and
+@code{ex::diff}):
+
+@example
+static ex cos_deriv(const ex & x, unsigned diff_param)
+@{
+ return -sin(x);
+@}
+@end example
+
+@cindex product rule
+The second parameter is obligatory but uninteresting at this point. It
+specifies which parameter to differentiate in a partial derivative in
+case the function has more than one parameter, and its main application
+is for correct handling of the chain rule.
+
+An implementation of the series expansion is not needed for @code{cos()} as
+it doesn't have any poles and GiNaC can do Taylor expansion by itself (as
+long as it knows what the derivative of @code{cos()} is). @code{tan()}, on
+the other hand, does have poles and may need to do Laurent expansion:
+
+@example
+static ex tan_series(const ex & x, const relational & rel,
+ int order, unsigned options)
+@{
+ // Find the actual expansion point
+ const ex x_pt = x.subs(rel);
+
+ if ("x_pt is not an odd multiple of Pi/2")
+ throw do_taylor(); // tell function::series() to do Taylor expansion
+
+ // On a pole, expand sin()/cos()
+ return (sin(x)/cos(x)).series(rel, order+2, options);
+@}
+@end example
+
+The @code{series()} implementation of a function @emph{must} return a
+@code{pseries} object, otherwise your code will crash.
+
+@subsection Function options
+
+GiNaC functions understand several more options which are always
+specified as @code{.option(params)}. None of them are required, but you
+need to specify at least one option to @code{REGISTER_FUNCTION()}. There
+is a do-nothing option called @code{dummy()} which you can use to define
+functions without any special options.
+
+@example
+eval_func(<C++ function>)
+evalf_func(<C++ function>)
+derivative_func(<C++ function>)
+series_func(<C++ function>)
+conjugate_func(<C++ function>)
+@end example
+
+These specify the C++ functions that implement symbolic evaluation,
+numeric evaluation, partial derivatives, and series expansion, respectively.
+They correspond to the GiNaC methods @code{eval()}, @code{evalf()},
+@code{diff()} and @code{series()}.
+
+The @code{eval_func()} function needs to use @code{.hold()} if no further
+automatic evaluation is desired or possible.
+
+If no @code{series_func()} is given, GiNaC defaults to simple Taylor
+expansion, which is correct if there are no poles involved. If the function
+has poles in the complex plane, the @code{series_func()} needs to check
+whether the expansion point is on a pole and fall back to Taylor expansion
+if it isn't. Otherwise, the pole usually needs to be regularized by some
+suitable transformation.
+
+@example
+latex_name(const string & n)
+@end example
+
+specifies the LaTeX code that represents the name of the function in LaTeX
+output. The default is to put the function name in an @code{\mbox@{@}}.
+
+@example
+do_not_evalf_params()
+@end example
+
+This tells @code{evalf()} to not recursively evaluate the parameters of the
+function before calling the @code{evalf_func()}.
+
+@example
+set_return_type(unsigned return_type, unsigned return_type_tinfo)
+@end example
+
+This allows you to explicitly specify the commutation properties of the
+function (@xref{Non-commutative objects}, for an explanation of
+(non)commutativity in GiNaC). For example, you can use
+@code{set_return_type(return_types::noncommutative, TINFO_matrix)} to make
+GiNaC treat your function like a matrix. By default, functions inherit the
+commutation properties of their first argument.
+
+@example
+set_symmetry(const symmetry & s)
+@end example
+
+specifies the symmetry properties of the function with respect to its
+arguments. @xref{Indexed objects}, for an explanation of symmetry
+specifications. GiNaC will automatically rearrange the arguments of
+symmetric functions into a canonical order.
+
+Sometimes you may want to have finer control over how functions are
+displayed in the output. For example, the @code{abs()} function prints
+itself as @samp{abs(x)} in the default output format, but as @samp{|x|}
+in LaTeX mode, and @code{fabs(x)} in C source output. This is achieved
+with the
+
+@example
+print_func<C>(<C++ function>)
+@end example
+
+option which is explained in the next section.
+
+
+@node Printing, Structures, Symbolic functions, Extending GiNaC
+@c node-name, next, previous, up
+@section GiNaC's expression output system
+
+GiNaC allows the output of expressions in a variety of different formats
+(@pxref{Input/Output}). This section will explain how expression output
+is implemented internally, and how to define your own output formats or
+change the output format of built-in algebraic objects. You will also want
+to read this section if you plan to write your own algebraic classes or
+functions.
+
+@cindex @code{print_context} (class)
+@cindex @code{print_dflt} (class)
+@cindex @code{print_latex} (class)
+@cindex @code{print_tree} (class)
+@cindex @code{print_csrc} (class)
+All the different output formats are represented by a hierarchy of classes
+rooted in the @code{print_context} class, defined in the @file{print.h}
+header file:
+
+@table @code
+@item print_dflt
+the default output format
+@item print_latex
+output in LaTeX mathematical mode
+@item print_tree
+a dump of the internal expression structure (for debugging)
+@item print_csrc
+the base class for C source output
+@item print_csrc_float
+C source output using the @code{float} type
+@item print_csrc_double
+C source output using the @code{double} type
+@item print_csrc_cl_N
+C source output using CLN types
+@end table
+
+The @code{print_context} base class provides two public data members:
+
+@example
+class print_context
+@{
+ ...
+public:
+ std::ostream & s;
+ unsigned options;
+@};
+@end example
+
+@code{s} is a reference to the stream to output to, while @code{options}
+holds flags and modifiers. Currently, there is only one flag defined:
+@code{print_options::print_index_dimensions} instructs the @code{idx} class
+to print the index dimension which is normally hidden.
+
+When you write something like @code{std::cout << e}, where @code{e} is
+an object of class @code{ex}, GiNaC will construct an appropriate
+@code{print_context} object (of a class depending on the selected output
+format), fill in the @code{s} and @code{options} members, and call
+
+@cindex @code{print()}
+@example
+void ex::print(const print_context & c, unsigned level = 0) const;
+@end example
+
+which in turn forwards the call to the @code{print()} method of the
+top-level algebraic object contained in the expression.
+
+Unlike other methods, GiNaC classes don't usually override their
+@code{print()} method to implement expression output. Instead, the default
+implementation @code{basic::print(c, level)} performs a run-time double
+dispatch to a function selected by the dynamic type of the object and the
+passed @code{print_context}. To this end, GiNaC maintains a separate method
+table for each class, similar to the virtual function table used for ordinary
+(single) virtual function dispatch.
+
+The method table contains one slot for each possible @code{print_context}
+type, indexed by the (internally assigned) serial number of the type. Slots
+may be empty, in which case GiNaC will retry the method lookup with the
+@code{print_context} object's parent class, possibly repeating the process
+until it reaches the @code{print_context} base class. If there's still no
+method defined, the method table of the algebraic object's parent class
+is consulted, and so on, until a matching method is found (eventually it
+will reach the combination @code{basic/print_context}, which prints the
+object's class name enclosed in square brackets).
+
+You can think of the print methods of all the different classes and output
+formats as being arranged in a two-dimensional matrix with one axis listing
+the algebraic classes and the other axis listing the @code{print_context}
+classes.
+
+Subclasses of @code{basic} can, of course, also overload @code{basic::print()}
+to implement printing, but then they won't get any of the benefits of the
+double dispatch mechanism (such as the ability for derived classes to
+inherit only certain print methods from its parent, or the replacement of
+methods at run-time).
+
+@subsection Print methods for classes
+
+The method table for a class is set up either in the definition of the class,
+by passing the appropriate @code{print_func<C>()} option to
+@code{GINAC_IMPLEMENT_REGISTERED_CLASS_OPT()} (@xref{Adding classes}, for
+an example), or at run-time using @code{set_print_func<T, C>()}. The latter
+can also be used to override existing methods dynamically.
+
+The argument to @code{print_func<C>()} and @code{set_print_func<T, C>()} can
+be a member function of the class (or one of its parent classes), a static
+member function, or an ordinary (global) C++ function. The @code{C} template
+parameter specifies the appropriate @code{print_context} type for which the
+method should be invoked, while, in the case of @code{set_print_func<>()}, the
+@code{T} parameter specifies the algebraic class (for @code{print_func<>()},
+the class is the one being implemented by
+@code{GINAC_IMPLEMENT_REGISTERED_CLASS_OPT}).
+
+For print methods that are member functions, their first argument must be of
+a type convertible to a @code{const C &}, and the second argument must be an
+@code{unsigned}.
+
+For static members and global functions, the first argument must be of a type
+convertible to a @code{const T &}, the second argument must be of a type
+convertible to a @code{const C &}, and the third argument must be an
+@code{unsigned}. A global function will, of course, not have access to
+private and protected members of @code{T}.
+
+The @code{unsigned} argument of the print methods (and of @code{ex::print()}
+and @code{basic::print()}) is used for proper parenthesizing of the output
+(and by @code{print_tree} for proper indentation). It can be used for similar
+purposes if you write your own output formats.
+
+The explanations given above may seem complicated, but in practice it's
+really simple, as shown in the following example. Suppose that we want to
+display exponents in LaTeX output not as superscripts but with little
+upwards-pointing arrows. This can be achieved in the following way:
+
+@example
+void my_print_power_as_latex(const power & p,
+ const print_latex & c,
+ unsigned level)
+@{
+ // get the precedence of the 'power' class
+ unsigned power_prec = p.precedence();
+
+ // if the parent operator has the same or a higher precedence
+ // we need parentheses around the power
+ if (level >= power_prec)
+ c.s << '(';
+
+ // print the basis and exponent, each enclosed in braces, and
+ // separated by an uparrow
+ c.s << '@{';
+ p.op(0).print(c, power_prec);
+ c.s << "@}\\uparrow@{";
+ p.op(1).print(c, power_prec);
+ c.s << '@}';
+
+ // don't forget the closing parenthesis
+ if (level >= power_prec)
+ c.s << ')';
+@}
+
+int main()
+@{
+ // a sample expression
+ symbol x("x"), y("y");
+ ex e = -3*pow(x, 3)*pow(y, -2) + pow(x+y, 2) - 1;
+
+ // switch to LaTeX mode
+ cout << latex;
+
+ // this prints "-1+@{(y+x)@}^@{2@}-3 \frac@{x^@{3@}@}@{y^@{2@}@}"
+ cout << e << endl;
+
+ // now we replace the method for the LaTeX output of powers with
+ // our own one
+ set_print_func<power, print_latex>(my_print_power_as_latex);
+
+ // this prints "-1+@{@{(y+x)@}@}\uparrow@{2@}-3 \frac@{@{x@}\uparrow@{3@}@}@{@{y@}\uparrow@{2@}@}"
+ cout << e << endl;
+@}
+@end example
+
+Some notes:
+
+@itemize
+
+@item
+The first argument of @code{my_print_power_as_latex} could also have been
+a @code{const basic &}, the second one a @code{const print_context &}.
+
+@item
+The above code depends on @code{mul} objects converting their operands to
+@code{power} objects for the purpose of printing.
+
+@item
+The output of products including negative powers as fractions is also
+controlled by the @code{mul} class.
+
+@item
+The @code{power/print_latex} method provided by GiNaC prints square roots
+using @code{\sqrt}, but the above code doesn't.
+
+@end itemize
+
+It's not possible to restore a method table entry to its previous or default
+value. Once you have called @code{set_print_func()}, you can only override
+it with another call to @code{set_print_func()}, but you can't easily go back
+to the default behavior again (you can, of course, dig around in the GiNaC
+sources, find the method that is installed at startup
+(@code{power::do_print_latex} in this case), and @code{set_print_func} that
+one; that is, after you circumvent the C++ member access control@dots{}).
+
+@subsection Print methods for functions
+
+Symbolic functions employ a print method dispatch mechanism similar to the
+one used for classes. The methods are specified with @code{print_func<C>()}
+function options. If you don't specify any special print methods, the function
+will be printed with its name (or LaTeX name, if supplied), followed by a
+comma-separated list of arguments enclosed in parentheses.
+
+For example, this is what GiNaC's @samp{abs()} function is defined like:
+
+@example
+static ex abs_eval(const ex & arg) @{ ... @}
+static ex abs_evalf(const ex & arg) @{ ... @}
+
+static void abs_print_latex(const ex & arg, const print_context & c)
+@{
+ c.s << "@{|"; arg.print(c); c.s << "|@}";
+@}
+
+static void abs_print_csrc_float(const ex & arg, const print_context & c)
+@{
+ c.s << "fabs("; arg.print(c); c.s << ")";
+@}
+
+REGISTER_FUNCTION(abs, eval_func(abs_eval).
+ evalf_func(abs_evalf).
+ print_func<print_latex>(abs_print_latex).
+ print_func<print_csrc_float>(abs_print_csrc_float).
+ print_func<print_csrc_double>(abs_print_csrc_float));
+@end example
+
+This will display @samp{abs(x)} as @samp{|x|} in LaTeX mode and @code{fabs(x)}
+in non-CLN C source output, but as @code{abs(x)} in all other formats.
+
+There is currently no equivalent of @code{set_print_func()} for functions.
+
+@subsection Adding new output formats
+
+Creating a new output format involves subclassing @code{print_context},
+which is somewhat similar to adding a new algebraic class
+(@pxref{Adding classes}). There is a macro @code{GINAC_DECLARE_PRINT_CONTEXT}
+that needs to go into the class definition, and a corresponding macro
+@code{GINAC_IMPLEMENT_PRINT_CONTEXT} that has to appear at global scope.
+Every @code{print_context} class needs to provide a default constructor
+and a constructor from an @code{std::ostream} and an @code{unsigned}
+options value.
+
+Here is an example for a user-defined @code{print_context} class:
+
+@example
+class print_myformat : public print_dflt
+@{
+ GINAC_DECLARE_PRINT_CONTEXT(print_myformat, print_dflt)
+public:
+ print_myformat(std::ostream & os, unsigned opt = 0)
+ : print_dflt(os, opt) @{@}
+@};
+
+print_myformat::print_myformat() : print_dflt(std::cout) @{@}
+
+GINAC_IMPLEMENT_PRINT_CONTEXT(print_myformat, print_dflt)
+@end example
+
+That's all there is to it. None of the actual expression output logic is
+implemented in this class. It merely serves as a selector for choosing
+a particular format. The algorithms for printing expressions in the new
+format are implemented as print methods, as described above.
+
+@code{print_myformat} is a subclass of @code{print_dflt}, so it behaves
+exactly like GiNaC's default output format:
+
+@example
+@{
+ symbol x("x");
+ ex e = pow(x, 2) + 1;
+
+ // this prints "1+x^2"
+ cout << e << endl;
+
+ // this also prints "1+x^2"
+ e.print(print_myformat()); cout << endl;
+
+ ...
+@}
+@end example
+
+To fill @code{print_myformat} with life, we need to supply appropriate
+print methods with @code{set_print_func()}, like this:
+
+@example
+// This prints powers with '**' instead of '^'. See the LaTeX output
+// example above for explanations.
+void print_power_as_myformat(const power & p,
+ const print_myformat & c,
+ unsigned level)
+@{
+ unsigned power_prec = p.precedence();
+ if (level >= power_prec)
+ c.s << '(';
+ p.op(0).print(c, power_prec);
+ c.s << "**";
+ p.op(1).print(c, power_prec);
+ if (level >= power_prec)
+ c.s << ')';
+@}
+
+@{
+ ...
+ // install a new print method for power objects
+ set_print_func<power, print_myformat>(print_power_as_myformat);
+
+ // now this prints "1+x**2"
+ e.print(print_myformat()); cout << endl;
+
+ // but the default format is still "1+x^2"
+ cout << e << endl;
+@}
+@end example
+
+
+@node Structures, Adding classes, Printing, Extending GiNaC
+@c node-name, next, previous, up
+@section Structures
+
+If you are doing some very specialized things with GiNaC, or if you just
+need some more organized way to store data in your expressions instead of
+anonymous lists, you may want to implement your own algebraic classes.
+('algebraic class' means any class directly or indirectly derived from
+@code{basic} that can be used in GiNaC expressions).
+
+GiNaC offers two ways of accomplishing this: either by using the
+@code{structure<T>} template class, or by rolling your own class from
+scratch. This section will discuss the @code{structure<T>} template which
+is easier to use but more limited, while the implementation of custom
+GiNaC classes is the topic of the next section. However, you may want to
+read both sections because many common concepts and member functions are
+shared by both concepts, and it will also allow you to decide which approach
+is most suited to your needs.
+
+The @code{structure<T>} template, defined in the GiNaC header file
+@file{structure.h}, wraps a type that you supply (usually a C++ @code{struct}
+or @code{class}) into a GiNaC object that can be used in expressions.
+
+@subsection Example: scalar products
+
+Let's suppose that we need a way to handle some kind of abstract scalar
+product of the form @samp{<x|y>} in expressions. Objects of the scalar
+product class have to store their left and right operands, which can in turn
+be arbitrary expressions. Here is a possible way to represent such a
+product in a C++ @code{struct}:
+
+@example
+#include <iostream>
+using namespace std;
+
+#include <ginac/ginac.h>
+using namespace GiNaC;
+
+struct sprod_s @{
+ ex left, right;
+
+ sprod_s() @{@}
+ sprod_s(ex l, ex r) : left(l), right(r) @{@}
+@};
+@end example
+
+The default constructor is required. Now, to make a GiNaC class out of this
+data structure, we need only one line:
+
+@example
+typedef structure<sprod_s> sprod;
+@end example
+
+That's it. This line constructs an algebraic class @code{sprod} which
+contains objects of type @code{sprod_s}. We can now use @code{sprod} in
+expressions like any other GiNaC class:
+
+@example
+...
+ symbol a("a"), b("b");
+ ex e = sprod(sprod_s(a, b));
+...