This is a tutorial that documents GiNaC @value{VERSION}, an open
framework for symbolic computation within the C++ programming language.
-Copyright (C) 1999-2003 Johannes Gutenberg University Mainz, Germany
+Copyright (C) 1999-2004 Johannes Gutenberg University Mainz, Germany
Permission is granted to make and distribute verbatim copies of
this manual provided the copyright notice and this permission notice
@subtitle An open framework for symbolic computation within the C++ programming language
@subtitle @value{UPDATED}
@author The GiNaC Group:
-@author Christian Bauer, Alexander Frink, Richard Kreckel
+@author Christian Bauer, Alexander Frink, Richard Kreckel, Jens Vollinga
@page
@vskip 0pt plus 1filll
-Copyright @copyright{} 1999-2003 Johannes Gutenberg University Mainz, Germany
+Copyright @copyright{} 1999-2004 Johannes Gutenberg University Mainz, Germany
@sp 2
Permission is granted to make and distribute verbatim copies of
this manual provided the copyright notice and this permission notice
@section License
The GiNaC framework for symbolic computation within the C++ programming
-language is Copyright @copyright{} 1999-2003 Johannes Gutenberg
+language is Copyright @copyright{} 1999-2004 Johannes Gutenberg
University Mainz, Germany.
This program is free software; you can redistribute it and/or
* Matrices:: Matrices.
* Indexed objects:: Handling indexed quantities.
* Non-commutative objects:: Algebras with non-commutative products.
+* Hash Maps:: A faster alternative to std::map<>.
@end menu
hierarchy and describe the classes of objects that are handled by
@code{ex}.
+@subsection Note: Expressions and STL containers
+
+GiNaC expressions (@code{ex} objects) have value semantics (they can be
+assigned, reassigned and copied like integral types) but the operator
+@code{<} doesn't provide a well-defined ordering on them. In STL-speak,
+expressions are @samp{Assignable} but not @samp{LessThanComparable}.
+
+This implies that in order to use expressions in sorted containers such as
+@code{std::map<>} and @code{std::set<>} you have to supply a suitable
+comparison predicate. GiNaC provides such a predicate, called
+@code{ex_is_less}. For example, a set of expressions should be defined
+as @code{std::set<ex, ex_is_less>}.
+
+Unsorted containers such as @code{std::vector<>} and @code{std::list<>}
+don't pose a problem. A @code{std::vector<ex>} works as expected.
+
+@xref{Information About Expressions}, for more about comparing and ordering
+expressions.
+
@node Automatic evaluation, Error handling, Expressions, Basic Concepts
@c node-name, next, previous, up
@itemize @bullet
@item
-at most of complexity @math{O(n log n)}
+at most of complexity
+@tex
+$O(n\log n)$
+@end tex
+@ifnottex
+@math{O(n log n)}
+@end ifnottex
@item
algebraically correct, possibly except for a set of measure zero (e.g.
@math{x/x} is transformed to @math{1} although this is incorrect for @math{x=0})
The terms of sums and products (and some other things like the arguments of
symmetric functions, the indices of symmetric tensors etc.) are re-ordered
into a canonical form that is deterministic, but not lexicographical or in
-any other way easily guessable (it almost always depends on the number and
+any other way easy to guess (it almost always depends on the number and
order of the symbols you define). However, constructing the same expression
twice, either implicitly or explicitly, will always result in the same
canonical form.
The general rule is that when you construct expressions, GiNaC automatically
creates them in canonical form, which might differ from the form you typed in
your program. This may create some awkward looking output (@samp{-y+x} instead
-of @samp{y-x}) but allows for more efficient operation and usually yields
+of @samp{x-y}) but allows for more efficient operation and usually yields
some immediate simplifications.
@cindex @code{eval()}
The @code{pole_error} class has a member function
@example
-int pole_error::degree(void) const;
+int pole_error::degree() const;
@end example
that returns the order of the singularity (or 0 when the pole is
logarithmic or the order is undefined).
-When using GiNaC it is useful to arrange for exceptions to be catched in
+When using GiNaC it is useful to arrange for exceptions to be caught in
the main program even if you don't want to do any special error handling.
Otherwise whenever an error occurs in GiNaC, it will be delegated to the
default exception handler of your C++ compiler's run-time system which
using namespace std;
using namespace GiNaC;
-int main(void)
+int main()
@{
try @{
...
@end ifnottex
@dots{}
@item @code{pseries} @tab Power Series, e.g. @math{x-1/6*x^3+1/120*x^5+O(x^7)}
-@item @code{function} @tab A symbolic function like @math{sin(2*x)}
+@item @code{function} @tab A symbolic function like
+@tex
+$\sin 2x$
+@end tex
+@ifnottex
+@math{sin(2*x)}
+@end ifnottex
@item @code{lst} @tab Lists of expressions @{@math{x}, @math{2*y}, @math{3+z}@}
@item @code{matrix} @tab @math{m}x@math{n} matrices of expressions
@item @code{relational} @tab A relation like the identity @math{x}@code{==}@math{y}
@item @code{varidx} @tab Index with variance
@item @code{spinidx} @tab Index with variance and dot (used in Weyl-van-der-Waerden spinor formalism)
@item @code{wildcard} @tab Wildcard for pattern matching
+@item @code{structure} @tab Template for user-defined classes
@end multitable
@end cartouche
in the calling function. Again, comparing them (using @code{operator==}
for instance) will always reveal their difference. Watch out, please.
+@cindex @code{realsymbol()}
+Symbols are expected to stand in for complex values by default, i.e. they live
+in the complex domain. As a consequence, operations like complex conjugation,
+for example (see @ref{Complex Conjugation}), do @emph{not} evaluate if applied
+to such symbols. Likewise @code{log(exp(x))} does not evaluate to @code{x},
+because of the unknown imaginary part of @code{x}.
+On the other hand, if you are sure that your symbols will hold only real values, you
+would like to have such functions evaluated. Therefore GiNaC allows you to specify
+the domain of the symbol. Instead of @code{symbol x("x");} you can write
+@code{realsymbol x("x");} to tell GiNaC that @code{x} stands in for real values.
+
@cindex @code{subs()}
Although symbols can be assigned expressions for internal reasons, you
should not do it (and we are not going to tell you how it is done). If
@end multitable
@end cartouche
+@subsection Converting numbers
+
+Sometimes it is desirable to convert a @code{numeric} object back to a
+built-in arithmetic type (@code{int}, @code{double}, etc.). The @code{numeric}
+class provides a couple of methods for this purpose:
+
+@cindex @code{to_int()}
+@cindex @code{to_long()}
+@cindex @code{to_double()}
+@cindex @code{to_cl_N()}
+@example
+int numeric::to_int() const;
+long numeric::to_long() const;
+double numeric::to_double() const;
+cln::cl_N numeric::to_cl_N() const;
+@end example
+
+@code{to_int()} and @code{to_long()} only work when the number they are
+applied on is an exact integer. Otherwise the program will halt with a
+message like @samp{Not a 32-bit integer}. @code{to_double()} applied on a
+rational number will return a floating-point approximation. Both
+@code{to_int()/to_long()} and @code{to_double()} discard the imaginary
+part of complex numbers.
+
@node Constants, Fundamental containers, Numbers, Basic Concepts
@c node-name, next, previous, up
@cindex @code{prepend()}
@cindex @code{remove_first()}
@cindex @code{remove_last()}
+@cindex @code{remove_all()}
The GiNaC class @code{lst} serves for holding a @dfn{list} of arbitrary
expressions. They are not as ubiquitous as in many other computer algebra
packages, but are sometimes used to supply a variable number of arguments of
-the same type to GiNaC methods such as @code{subs()} and @code{to_rational()},
-so you should have a basic understanding of them.
+the same type to GiNaC methods such as @code{subs()} and some @code{matrix}
+constructors, so you should have a basic understanding of them.
-Lists of up to 16 expressions can be directly constructed from single
+Lists can be constructed by assigning a comma-separated sequence of
expressions:
@example
@{
symbol x("x"), y("y");
- lst l(x, 2, y, x+y);
- // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y'
+ lst l;
+ l = x, 2, y, x+y;
+ // now, l is a list holding the expressions 'x', '2', 'y', and 'x+y',
+ // in that order
+ ...
+@end example
+
+There are also constructors that allow direct creation of lists of up to
+16 expressions, which is often more convenient but slightly less efficient:
+
+@example
+ ...
+ // This produces the same list 'l' as above:
+ // lst l(x, 2, y, x+y);
+ // lst l = lst(x, 2, y, x+y);
...
@end example
...
@end example
+As with the standard @code{list<T>} container, accessing random elements of a
+@code{lst} is generally an operation of order @math{O(N)}. Faster read-only
+sequential access to the elements of a list is possible with the
+iterator types provided by the @code{lst} class:
+
+@example
+typedef ... lst::const_iterator;
+typedef ... lst::const_reverse_iterator;
+lst::const_iterator lst::begin() const;
+lst::const_iterator lst::end() const;
+lst::const_reverse_iterator lst::rbegin() const;
+lst::const_reverse_iterator lst::rend() const;
+@end example
+
+For example, to print the elements of a list individually you can use:
+
+@example
+ ...
+ // O(N)
+ for (lst::const_iterator i = l.begin(); i != l.end(); ++i)
+ cout << *i << endl;
+ ...
+@end example
+
+which is one order faster than
+
+@example
+ ...
+ // O(N^2)
+ for (size_t i = 0; i < l.nops(); ++i)
+ cout << l.op(i) << endl;
+ ...
+@end example
+
+These iterators also allow you to use some of the algorithms provided by
+the C++ standard library:
+
+@example
+ ...
+ // print the elements of the list (requires #include <iterator>)
+ std::copy(l.begin(), l.end(), ostream_iterator<ex>(cout, "\n"));
+
+ // sum up the elements of the list (requires #include <numeric>)
+ ex sum = std::accumulate(l.begin(), l.end(), ex(0));
+ cout << sum << endl; // prints '2+2*x+2*y'
+ ...
+@end example
+
@code{lst} is one of the few GiNaC classes that allow in-place modifications
(the only other one is @code{matrix}). You can modify single elements:
@example
...
+ l[1] = 42; // l is now @{x, 42, y, x+y@}
l.let_op(1) = 7; // l is now @{x, 7, y, x+y@}
...
@end example
...
@end example
+You can remove all the elements of a list with @code{remove_all()}:
+
+@example
+ ...
+ l.remove_all(); // l is now empty
+ ...
+@end example
+
You can bring the elements of a list into a canonical order with @code{sort()}:
@example
...
- lst l1(x, 2, y, x+y);
- lst l2(2, x+y, x, y);
+ lst l1, l2;
+ l1 = x, 2, y, x+y;
+ l2 = 2, x+y, x, y;
l1.sort();
l2.sort();
// l1 and l2 are now equal
@example
...
- lst l3(x, 2, 2, 2, y, x+y, y+x);
+ lst l3;
+ l3 = x, 2, 2, 2, y, x+y, y+x;
l3.unique(); // l3 is now @{x, 2, y, x+y@}
@}
@end example
second one in the range 0@dots{}@math{n-1}.
There are a couple of ways to construct matrices, with or without preset
-elements:
+elements. The constructor
+
+@example
+matrix::matrix(unsigned r, unsigned c);
+@end example
+
+creates a matrix with @samp{r} rows and @samp{c} columns with all elements
+set to zero.
+
+The fastest way to create a matrix with preinitialized elements is to assign
+a list of comma-separated expressions to an empty matrix (see below for an
+example). But you can also specify the elements as a (flat) list with
+
+@example
+matrix::matrix(unsigned r, unsigned c, const lst & l);
+@end example
+
+The function
@cindex @code{lst_to_matrix()}
+@example
+ex lst_to_matrix(const lst & l);
+@end example
+
+constructs a matrix from a list of lists, each list representing a matrix row.
+
+There is also a set of functions for creating some special types of
+matrices:
+
@cindex @code{diag_matrix()}
@cindex @code{unit_matrix()}
@cindex @code{symbolic_matrix()}
@example
-matrix::matrix(unsigned r, unsigned c);
-matrix::matrix(unsigned r, unsigned c, const lst & l);
-ex lst_to_matrix(const lst & l);
ex diag_matrix(const lst & l);
ex unit_matrix(unsigned x);
ex unit_matrix(unsigned r, unsigned c);
ex symbolic_matrix(unsigned r, unsigned c, const string & base_name, const string & tex_base_name);
@end example
-The first two functions are @code{matrix} constructors which create a matrix
-with @samp{r} rows and @samp{c} columns. The matrix elements can be
-initialized from a (flat) list of expressions @samp{l}. Otherwise they are
-all set to zero. The @code{lst_to_matrix()} function constructs a matrix
-from a list of lists, each list representing a matrix row. @code{diag_matrix()}
-constructs a diagonal matrix given the list of diagonal elements.
-@code{unit_matrix()} creates an @samp{x} by @samp{x} (or @samp{r} by @samp{c})
-unit matrix. And finally, @code{symbolic_matrix} constructs a matrix filled
-with newly generated symbols made of the specified base name and the
-position of each element in the matrix.
+@code{diag_matrix()} constructs a diagonal matrix given the list of diagonal
+elements. @code{unit_matrix()} creates an @samp{x} by @samp{x} (or @samp{r}
+by @samp{c}) unit matrix. And finally, @code{symbolic_matrix} constructs a
+matrix filled with newly generated symbols made of the specified base name
+and the position of each element in the matrix.
Matrix elements can be accessed and set using the parenthesis (function call)
operator:
the @code{op()} method. But C++-style subscripting with square brackets
@samp{[]} is not available.
-Here are a couple of examples of constructing matrices:
+Here are a couple of examples for constructing matrices:
@example
@{
symbol a("a"), b("b");
matrix M(2, 2);
- M(0, 0) = a;
- M(1, 1) = b;
+ M = a, 0,
+ 0, b;
cout << M << endl;
// -> [[a,0],[0,b]]
+ matrix M2(2, 2);
+ M2(0, 0) = a;
+ M2(1, 1) = b;
+ cout << M2 << endl;
+ // -> [[a,0],[0,b]]
+
cout << matrix(2, 2, lst(a, 0, 0, b)) << endl;
// -> [[a,0],[0,b]]
@end example
@cindex @code{transpose()}
-@cindex @code{inverse()}
There are three ways to do arithmetic with matrices. The first (and most
-efficient one) is to use the methods provided by the @code{matrix} class:
+direct one) is to use the methods provided by the @code{matrix} class:
@example
matrix matrix::add(const matrix & other) const;
matrix matrix::mul(const matrix & other) const;
matrix matrix::mul_scalar(const ex & other) const;
matrix matrix::pow(const ex & expn) const;
-matrix matrix::transpose(void) const;
-matrix matrix::inverse(void) const;
+matrix matrix::transpose() const;
@end example
All of these methods return the result as a new matrix object. Here is an
@example
@{
- matrix A(2, 2, lst(1, 2, 3, 4));
- matrix B(2, 2, lst(-1, 0, 2, 1));
- matrix C(2, 2, lst(8, 4, 2, 1));
+ matrix A(2, 2), B(2, 2), C(2, 2);
+ A = 1, 2,
+ 3, 4;
+ B = -1, 0,
+ 2, 1;
+ C = 8, 4,
+ 2, 1;
matrix result = A.mul(B).sub(C.mul_scalar(2));
cout << result << endl;
@cindex @code{trace()}
@cindex @code{charpoly()}
@example
-ex matrix::determinant(unsigned algo = determinant_algo::automatic) const;
-ex matrix::trace(void) const;
-ex matrix::charpoly(const symbol & lambda) const;
+ex matrix::determinant(unsigned algo=determinant_algo::automatic) const;
+ex matrix::trace() const;
+ex matrix::charpoly(const ex & lambda) const;
@end example
-The @samp{algo} argument of @code{determinant()} allows to select between
-different algorithms for calculating the determinant. The possible values
-are defined in the @file{flags.h} header file. By default, GiNaC uses a
-heuristic to automatically select an algorithm that is likely to give the
-result most quickly.
+The @samp{algo} argument of @code{determinant()} allows to select
+between different algorithms for calculating the determinant. The
+asymptotic speed (as parametrized by the matrix size) can greatly differ
+between those algorithms, depending on the nature of the matrix'
+entries. The possible values are defined in the @file{flags.h} header
+file. By default, GiNaC uses a heuristic to automatically select an
+algorithm that is likely (but not guaranteed) to give the result most
+quickly.
+
+@cindex @code{inverse()}
+@cindex @code{solve()}
+Matrices may also be inverted using the @code{ex matrix::inverse()}
+method and linear systems may be solved with:
+
+@example
+matrix matrix::solve(const matrix & vars, const matrix & rhs, unsigned algo=solve_algo::automatic) const;
+@end example
+
+Assuming the matrix object this method is applied on is an @code{m}
+times @code{n} matrix, then @code{vars} must be a @code{n} times
+@code{p} matrix of symbolic indeterminates and @code{rhs} a @code{m}
+times @code{p} matrix. The returned matrix then has dimension @code{n}
+times @code{p} and in the case of an underdetermined system will still
+contain some of the indeterminates from @code{vars}. If the system is
+overdetermined, an exception is thrown.
@node Indexed objects, Non-commutative objects, Matrices, Basic Concepts
symbol A("A");
cout << indexed(A, i, j) << endl;
// -> A.i.j
+ cout << index_dimensions << indexed(A, i, j) << endl;
+ // -> A.i[3].j[3]
+ cout << dflt; // reset cout to default output format (dimensions hidden)
...
@end example
the symbol @code{A} as its base expression and the two indices @code{i} and
@code{j}.
+The dimensions of indices are normally not visible in the output, but one
+can request them to be printed with the @code{index_dimensions} manipulator,
+as shown above.
+
Note the difference between the indices @code{i} and @code{j} which are of
class @code{idx}, and the index values which are the symbols @code{i_sym}
and @code{j_sym}. The indices of indexed objects cannot directly be symbols
The methods
@example
-ex idx::get_value(void);
-ex idx::get_dimension(void);
+ex idx::get_value();
+ex idx::get_dimension();
@end example
return the value and dimension of an @code{idx} object. If you have an index
There are also the methods
@example
-bool idx::is_numeric(void);
-bool idx::is_symbolic(void);
-bool idx::is_dim_numeric(void);
-bool idx::is_dim_symbolic(void);
+bool idx::is_numeric();
+bool idx::is_symbolic();
+bool idx::is_dim_numeric();
+bool idx::is_dim_symbolic();
@end example
for checking whether the value and dimension are numeric or symbolic
constructor. The two methods
@example
-bool varidx::is_covariant(void);
-bool varidx::is_contravariant(void);
+bool varidx::is_covariant();
+bool varidx::is_contravariant();
@end example
allow you to check the variance of a @code{varidx} object (use @code{ex_to<varidx>()}
method
@example
-ex varidx::toggle_variance(void);
+ex varidx::toggle_variance();
@end example
which makes a new index with the same value and dimension but the opposite
methods
@example
-bool spinidx::is_dotted(void);
-bool spinidx::is_undotted(void);
+bool spinidx::is_dotted();
+bool spinidx::is_undotted();
@end example
allow you to check whether or not a @code{spinidx} object is dotted (use
Finally, the two methods
@example
-ex spinidx::toggle_dot(void);
-ex spinidx::toggle_variance_dot(void);
+ex spinidx::toggle_dot();
+ex spinidx::toggle_variance_dot();
@end example
create a new index with the same value and dimension but opposite dottedness
@end example
@cindex @code{get_free_indices()}
-@cindex Dummy index
+@cindex dummy index
@subsection Dummy indices
GiNaC treats certain symbolic index pairs as @dfn{dummy indices} meaning
there is the method
@example
-ex ex::simplify_indexed(void);
+ex ex::simplify_indexed();
ex ex::simplify_indexed(const scalar_products & sp);
@end example
symbol x("x"), y("y");
// A is a 2x2 matrix, X is a 2x1 vector
- matrix A(2, 2, lst(1, 2, 3, 4)), X(2, 1, lst(x, y));
+ matrix A(2, 2), X(2, 1);
+ A = 1, 2,
+ 3, 4;
+ X = x, y;
cout << indexed(A, i, i) << endl;
// -> 5
of the metric tensor.
-@node Non-commutative objects, Methods and Functions, Indexed objects, Basic Concepts
+@node Non-commutative objects, Hash Maps, Indexed objects, Basic Concepts
@c node-name, next, previous, up
@section Non-commutative objects
obtained with the two member functions
@example
-unsigned ex::return_type(void) const;
-unsigned ex::return_type_tinfo(void) const;
+unsigned ex::return_type() const;
+unsigned ex::return_type_tinfo() const;
@end example
The @code{return_type()} function returns one of three values (defined in
e = canonicalize_clifford(e);
cout << e << endl;
- // -> 2*eta~mu~nu
+ // -> 2*ONE*eta~mu~nu
@}
@end example
@end example
-@node Methods and Functions, Information About Expressions, Non-commutative objects, Top
+@node Hash Maps, Methods and Functions, Non-commutative objects, Basic Concepts
+@c node-name, next, previous, up
+@section Hash Maps
+@cindex hash maps
+@cindex @code{exhashmap} (class)
+
+For your convenience, GiNaC offers the container template @code{exhashmap<T>}
+that can be used as a drop-in replacement for the STL
+@code{std::map<ex, T, ex_is_less>}, using hash tables to provide faster,
+typically constant-time, element look-up than @code{map<>}.
+
+@code{exhashmap<>} supports all @code{map<>} members and operations, with the
+following differences:
+
+@itemize @bullet
+@item
+no @code{lower_bound()} and @code{upper_bound()} methods
+@item
+no reverse iterators, no @code{rbegin()}/@code{rend()}
+@item
+no @code{operator<(exhashmap, exhashmap)}
+@item
+the comparison function object @code{key_compare} is hardcoded to
+@code{ex_is_less}
+@item
+the constructor @code{exhashmap(size_t n)} allows specifying the minimum
+initial hash table size (the actual table size after construction may be
+larger than the specified value)
+@item
+the method @code{size_t bucket_count()} returns the current size of the hash
+table
+@item
+@code{insert()} and @code{erase()} operations invalidate all iterators
+@end itemize
+
+
+@node Methods and Functions, Information About Expressions, Hash Maps, Top
@c node-name, next, previous, up
@chapter Methods and Functions
@cindex polynomial
@menu
* Information About Expressions::
+* Numerical Evaluation::
* Substituting Expressions::
* Pattern Matching and Advanced Substitutions::
* Applying a Function on Subexpressions::
+* Visitors and Tree Traversal::
* Polynomial Arithmetic:: Working with polynomials.
* Rational Expressions:: Working with rational functions.
* Symbolic Differentiation::
* Series Expansion:: Taylor and Laurent expansion.
* Symmetrization::
* Built-in Functions:: List of predefined mathematical functions.
+* Multiple polylogarithms::
+* Complex Conjugation::
+* Built-in Functions:: List of predefined mathematical functions.
+* Solving Linear Systems of Equations::
* Input/Output:: Input and output of expressions.
@end menu
-@node Information About Expressions, Substituting Expressions, Methods and Functions, Methods and Functions
+@node Information About Expressions, Numerical Evaluation, Methods and Functions, Methods and Functions
@c node-name, next, previous, up
@section Getting information about expressions
bool is_a<T>(const ex & e);
bool is_exactly_a<T>(const ex & e);
bool ex::info(unsigned flag);
-unsigned ex::return_type(void) const;
-unsigned ex::return_type_tinfo(void) const;
+unsigned ex::return_type() const;
+unsigned ex::return_type_tinfo() const;
@end example
When the test made by @code{is_a<T>()} returns true, it is safe to call
@subsection Accessing subexpressions
+@cindex container
+
+Many GiNaC classes, like @code{add}, @code{mul}, @code{lst}, and
+@code{function}, act as containers for subexpressions. For example, the
+subexpressions of a sum (an @code{add} object) are the individual terms,
+and the subexpressions of a @code{function} are the function's arguments.
+
@cindex @code{nops()}
@cindex @code{op()}
-@cindex container
-@cindex @code{relational} (class)
+GiNaC provides several ways of accessing subexpressions. The first way is to
+use the two methods
+
+@example
+size_t ex::nops();
+ex ex::op(size_t i);
+@end example
+
+@code{nops()} determines the number of subexpressions (operands) contained
+in the expression, while @code{op(i)} returns the @code{i}-th
+(0..@code{nops()-1}) subexpression. In the case of a @code{power} object,
+@code{op(0)} will return the basis and @code{op(1)} the exponent. For
+@code{indexed} objects, @code{op(0)} is the base expression and @code{op(i)},
+@math{i>0} are the indices.
+
+@cindex iterators
+@cindex @code{const_iterator}
+The second way to access subexpressions is via the STL-style random-access
+iterator class @code{const_iterator} and the methods
+
+@example
+const_iterator ex::begin();
+const_iterator ex::end();
+@end example
+
+@code{begin()} returns an iterator referring to the first subexpression;
+@code{end()} returns an iterator which is one-past the last subexpression.
+If the expression has no subexpressions, then @code{begin() == end()}. These
+iterators can also be used in conjunction with non-modifying STL algorithms.
+
+Here is an example that (non-recursively) prints the subexpressions of a
+given expression in three different ways:
+
+@example
+@{
+ ex e = ...
+
+ // with nops()/op()
+ for (size_t i = 0; i != e.nops(); ++i)
+ cout << e.op(i) << endl;
-GiNaC provides the two methods
+ // with iterators
+ for (const_iterator i = e.begin(); i != e.end(); ++i)
+ cout << *i << endl;
+
+ // with iterators and STL copy()
+ std::copy(e.begin(), e.end(), std::ostream_iterator<ex>(cout, "\n"));
+@}
+@end example
+
+@cindex @code{const_preorder_iterator}
+@cindex @code{const_postorder_iterator}
+@code{op()}/@code{nops()} and @code{const_iterator} only access an
+expression's immediate children. GiNaC provides two additional iterator
+classes, @code{const_preorder_iterator} and @code{const_postorder_iterator},
+that iterate over all objects in an expression tree, in preorder or postorder,
+respectively. They are STL-style forward iterators, and are created with the
+methods
@example
-unsigned ex::nops();
-ex ex::op(unsigned i);
+const_preorder_iterator ex::preorder_begin();
+const_preorder_iterator ex::preorder_end();
+const_postorder_iterator ex::postorder_begin();
+const_postorder_iterator ex::postorder_end();
@end example
-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 following example illustrates the differences between
+@code{const_iterator}, @code{const_preorder_iterator}, and
+@code{const_postorder_iterator}:
+
+@example
+@{
+ symbol A("A"), B("B"), C("C");
+ ex e = lst(lst(A, B), C);
+
+ std::copy(e.begin(), e.end(),
+ std::ostream_iterator<ex>(cout, "\n"));
+ // @{A,B@}
+ // C
+
+ std::copy(e.preorder_begin(), e.preorder_end(),
+ std::ostream_iterator<ex>(cout, "\n"));
+ // @{@{A,B@},C@}
+ // @{A,B@}
+ // A
+ // B
+ // C
+
+ std::copy(e.postorder_begin(), e.postorder_end(),
+ std::ostream_iterator<ex>(cout, "\n"));
+ // A
+ // B
+ // @{A,B@}
+ // C
+ // @{@{A,B@},C@}
+@}
+@end example
-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
+@cindex @code{relational} (class)
+Finally, the left-hand side and right-hand side expressions of objects of
+class @code{relational} (and only of these) can also be accessed with the
+methods
@example
ex ex::lhs();
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.
+
+@subsection Ordering expressions
+@cindex @code{ex_is_less} (class)
+@cindex @code{ex_is_equal} (class)
+@cindex @code{compare()}
+
+Sometimes it is necessary to establish a mathematically well-defined ordering
+on a set of arbitrary expressions, for example to use expressions as keys
+in a @code{std::map<>} container, or to bring a vector of expressions into
+a canonical order (which is done internally by GiNaC for sums and products).
+
+The operators @code{<}, @code{>} etc. described in the last section cannot
+be used for this, as they don't implement an ordering relation in the
+mathematical sense. In particular, they are not guaranteed to be
+antisymmetric: if @samp{a} and @samp{b} are different expressions, and
+@code{a < b} yields @code{false}, then @code{b < a} doesn't necessarily
+yield @code{true}.
+
+By default, STL classes and algorithms use the @code{<} and @code{==}
+operators to compare objects, which are unsuitable for expressions, but GiNaC
+provides two functors that can be supplied as proper binary comparison
+predicates to the STL:
+
+@example
+class ex_is_less : public std::binary_function<ex, ex, bool> @{
+public:
+ bool operator()(const ex &lh, const ex &rh) const;
+@};
+
+class ex_is_equal : public std::binary_function<ex, ex, bool> @{
+public:
+ bool operator()(const ex &lh, const ex &rh) const;
+@};
+@end example
+
+For example, to define a @code{map} that maps expressions to strings you
+have to use
+
+@example
+std::map<ex, std::string, ex_is_less> myMap;
+@end example
+
+Omitting the @code{ex_is_less} template parameter will introduce spurious
+bugs because the map operates improperly.
+
+Other examples for the use of the functors:
+
+@example
+std::vector<ex> v;
+// fill vector
+...
+
+// sort vector
+std::sort(v.begin(), v.end(), ex_is_less());
+
+// count the number of expressions equal to '1'
+unsigned num_ones = std::count_if(v.begin(), v.end(),
+ std::bind2nd(ex_is_equal(), 1));
+@end example
+
+The implementation of @code{ex_is_less} uses the member function
+
+@example
+int ex::compare(const ex & other) const;
+@end example
+
+which returns @math{0} if @code{*this} and @code{other} are equal, @math{-1}
+if @code{*this} sorts before @code{other}, and @math{1} if @code{*this} sorts
+after @code{other}.
+
+
+@node Numerical Evaluation, Substituting Expressions, Information About Expressions, Methods and Functions
+@c node-name, next, previous, up
+@section Numerical Evaluation
+@cindex @code{evalf()}
+
+GiNaC keeps algebraic expressions, numbers and constants in their exact form.
+To evaluate them using floating-point arithmetic you need to call
+
+@example
+ex ex::evalf(int level = 0) const;
+@end example
+
+@cindex @code{Digits}
+The accuracy of the evaluation is controlled by the global object @code{Digits}
+which can be assigned an integer value. The default value of @code{Digits}
+is 17. @xref{Numbers}, for more information and examples.
+
+To evaluate an expression to a @code{double} floating-point number you can
+call @code{evalf()} followed by @code{numeric::to_double()}, like this:
+
+@example
+@{
+ // Approximate sin(x/Pi)
+ symbol x("x");
+ ex e = series(sin(x/Pi), x == 0, 6);
+
+ // Evaluate numerically at x=0.1
+ ex f = evalf(e.subs(x == 0.1));
+
+ // ex_to<numeric> is an unsafe cast, so check the type first
+ if (is_a<numeric>(f)) @{
+ double d = ex_to<numeric>(f).to_double();
+ cout << d << endl;
+ // -> 0.0318256
+ @} else
+ // error
+@}
+@end example
-@node Substituting Expressions, Pattern Matching and Advanced Substitutions, Information About Expressions, Methods and Functions
+@node Substituting Expressions, Pattern Matching and Advanced Substitutions, Numerical Evaluation, Methods and Functions
@c node-name, next, previous, up
@section Substituting expressions
@cindex @code{subs()}
expressions via the @code{.subs()} method:
@example
-ex ex::subs(const ex & e);
-ex ex::subs(const lst & syms, const lst & repls);
+ex ex::subs(const ex & e, unsigned options = 0);
+ex ex::subs(const exmap & m, unsigned options = 0);
+ex ex::subs(const lst & syms, const lst & repls, unsigned options = 0);
@end example
In the first form, @code{subs()} accepts a relational of the form
If you specify multiple substitutions, they are performed in parallel, so e.g.
@code{subs(lst(x == y, y == x))} exchanges @samp{x} and @samp{y}.
-The second form of @code{subs()} takes two lists, one for the objects to be
+The second form of @code{subs()} takes an @code{exmap} object which is a
+pair associative container that maps expressions to expressions (currently
+implemented as a @code{std::map}). This is the most efficient one of the
+three @code{subs()} forms and should be used when the number of objects to
+be substituted is large or unknown.
+
+Using this form, the second example from above would look like this:
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex e2 = x*y + x;
+
+ exmap m;
+ m[x] = -2;
+ m[y] = 4;
+ cout << "e2(-2, 4) = " << e2.subs(m) << endl;
+@}
+@end example
+
+The third form of @code{subs()} takes two lists, one for the objects to be
replaced and one for the expressions to be substituted (both lists must
contain the same number of elements). Using this form, you would write
-@code{subs(lst(x, y), lst(y, x))} to exchange @samp{x} and @samp{y}.
+
+@example
+@{
+ symbol x("x"), y("y");
+ ex e2 = x*y + x;
+
+ cout << "e2(-2, 4) = " << e2.subs(lst(x, y), lst(-2, 4)) << endl;
+@}
+@end example
+
+The optional last argument to @code{subs()} is a combination of
+@code{subs_options} flags. There are two options available:
+@code{subs_options::no_pattern} disables pattern matching, which makes
+large @code{subs()} operations significantly faster if you are not using
+patterns. The second option, @code{subs_options::algebraic} enables
+algebraic substitutions in products and powers.
+@ref{Pattern Matching and Advanced Substitutions}, for more information
+about patterns and algebraic substitutions.
@code{subs()} performs syntactic substitution of any complete algebraic
object; it does not try to match sub-expressions as is demonstrated by the
@end example
works a bit like @code{has()} but it doesn't stop upon finding the first
-match. Instead, it inserts all found matches into the specified list. If
-there are multiple occurrences of the same expression, it is entered only
-once to the list. @code{find()} returns false if no matches were found (in
+match. Instead, it appends all found matches to the specified list. If there
+are multiple occurrences of the same expression, it is entered only once to
+the list. @code{find()} returns false if no matches were found (in
@command{ginsh}, it returns an empty list):
@example
@}
@end example
+@subsection Algebraic substitutions
+Supplying the @code{subs_options::algebraic} option to @code{subs()}
+enables smarter, algebraic substitutions in products and powers. If you want
+to substitute some factors of a product, you only need to list these factors
+in your pattern. Furthermore, if an (integer) power of some expression occurs
+in your pattern and in the expression that you want the substitution to occur
+in, it can be substituted as many times as possible, without getting negative
+powers.
+
+An example clarifies it all (hopefully):
+
+@example
+cout << (a*a*a*a+b*b*b*b+pow(x+y,4)).subs(wild()*wild()==pow(wild(),3),
+ subs_options::algebraic) << endl;
+// --> (y+x)^6+b^6+a^6
+
+cout << ((a+b+c)*(a+b+c)).subs(a+b==x,subs_options::algebraic) << endl;
+// --> (c+b+a)^2
+// Powers and products are smart, but addition is just the same.
+
+cout << ((a+b+c)*(a+b+c)).subs(a+b+wild()==x+wild(), subs_options::algebraic)
+ << endl;
+// --> (x+c)^2
+// As I said: addition is just the same.
+
+cout << (pow(a,5)*pow(b,7)+2*b).subs(b*b*a==x,subs_options::algebraic) << endl;
+// --> x^3*b*a^2+2*b
+
+cout << (pow(a,-5)*pow(b,-7)+2*b).subs(1/(b*b*a)==x,subs_options::algebraic)
+ << endl;
+// --> 2*b+x^3*b^(-1)*a^(-2)
+
+cout << (4*x*x*x-2*x*x+5*x-1).subs(x==a,subs_options::algebraic) << endl;
+// --> -1-2*a^2+4*a^3+5*a
+
+cout << (4*x*x*x-2*x*x+5*x-1).subs(pow(x,wild())==pow(a,wild()),
+ subs_options::algebraic) << endl;
+// --> -1+5*x+4*x^3-2*x^2
+// You should not really need this kind of patterns very often now.
+// But perhaps this it's-not-a-bug-it's-a-feature (c/sh)ould still change.
+
+cout << ex(sin(1+sin(x))).subs(sin(wild())==cos(wild()),
+ subs_options::algebraic) << endl;
+// --> cos(1+cos(x))
+
+cout << expand((a*sin(x+y)*sin(x+y)+a*cos(x+y)*cos(x+y)+b)
+ .subs((pow(cos(wild()),2)==1-pow(sin(wild()),2)),
+ subs_options::algebraic)) << endl;
+// --> b+a
+@end example
+
-@node Applying a Function on Subexpressions, Polynomial Arithmetic, Pattern Matching and Advanced Substitutions, Methods and Functions
+@node Applying a Function on Subexpressions, Visitors and Tree Traversal, Pattern Matching and Advanced Substitutions, Methods and Functions
@c node-name, next, previous, up
@section Applying a Function on Subexpressions
-@cindex Tree traversal
+@cindex tree traversal
@cindex @code{map()}
Sometimes you may want to perform an operation on specific parts of an
return ex_to<matrix>(e).trace();
else if (is_a<add>(e)) @{
ex sum = 0;
- for (unsigned i=0; i<e.nops(); i++)
+ for (size_t i=0; i<e.nops(); i++)
sum += calc_trace(e.op(i));
return sum;
@} else if (is_a<mul>)(e)) @{
@end example
-@node Polynomial Arithmetic, Rational Expressions, Applying a Function on Subexpressions, Methods and Functions
+@node Visitors and Tree Traversal, Polynomial Arithmetic, Applying a Function on Subexpressions, Methods and Functions
@c node-name, next, previous, up
-@section Polynomial arithmetic
+@section Visitors and Tree Traversal
+@cindex tree traversal
+@cindex @code{visitor} (class)
+@cindex @code{accept()}
+@cindex @code{visit()}
+@cindex @code{traverse()}
+@cindex @code{traverse_preorder()}
+@cindex @code{traverse_postorder()}
-@subsection Expanding and collecting
-@cindex @code{expand()}
-@cindex @code{collect()}
-@cindex @code{collect_common_factors()}
+Suppose that you need a function that returns a list of all indices appearing
+in an arbitrary expression. The indices can have any dimension, and for
+indices with variance you always want the covariant version returned.
-A polynomial in one or more variables has many equivalent
-representations. Some useful ones serve a specific purpose. Consider
-for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
-21*y*z + 4*z^2} (written down here in output-style). It is equivalent
-to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}. Other
-representations are the recursive ones where one collects for exponents
-in one of the three variable. Since the factors are themselves
-polynomials in the remaining two variables the procedure can be
-repeated. In our example, two possibilities would be @math{(4*y + z)*x
-+ 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
-x*z}.
+You can't use @code{get_free_indices()} because you also want to include
+dummy indices in the list, and you can't use @code{find()} as it needs
+specific index dimensions (and it would require two passes: one for indices
+with variance, one for plain ones).
-To bring an expression into expanded form, its method
+The obvious solution to this problem is a tree traversal with a type switch,
+such as the following:
@example
-ex ex::expand();
-@end example
+void gather_indices_helper(const ex & e, lst & l)
+@{
+ if (is_a<varidx>(e)) @{
+ const varidx & vi = ex_to<varidx>(e);
+ l.append(vi.is_covariant() ? vi : vi.toggle_variance());
+ @} else if (is_a<idx>(e)) @{
+ l.append(e);
+ @} else @{
+ size_t n = e.nops();
+ for (size_t i = 0; i < n; ++i)
+ gather_indices_helper(e.op(i), l);
+ @}
+@}
-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!
+lst gather_indices(const ex & e)
+@{
+ lst l;
+ gather_indices_helper(e, l);
+ l.sort();
+ l.unique();
+ return l;
+@}
+@end example
-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:
+This works fine but fans of object-oriented programming will feel
+uncomfortable with the type switch. One reason is that there is a possibility
+for subtle bugs regarding derived classes. If we had, for example, written
@example
-ex ex::collect(const ex & s, bool distributed = false);
+ if (is_a<idx>(e)) @{
+ ...
+ @} else if (is_a<varidx>(e)) @{
+ ...
@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.
+in @code{gather_indices_helper}, the code wouldn't have worked because the
+first line "absorbs" all classes derived from @code{idx}, including
+@code{varidx}, so the special case for @code{varidx} would never have been
+executed.
-Note that the original polynomial needs to be in expanded form (for the
-variables concerned) in order for @code{collect()} to be able to find the
-coefficients properly.
+Also, for a large number of classes, a type switch like the above can get
+unwieldy and inefficient (it's a linear search, after all).
+@code{gather_indices_helper} only checks for two classes, but if you had to
+write a function that required a different implementation for nearly
+every GiNaC class, the result would be very hard to maintain and extend.
-The following @command{ginsh} transcript shows an application of @code{collect()}
-together with @code{find()}:
+The cleanest approach to the problem would be to add a new virtual function
+to GiNaC's class hierarchy. In our example, there would be specializations
+for @code{idx} and @code{varidx} while the default implementation in
+@code{basic} performed the tree traversal. Unfortunately, in C++ it's
+impossible to add virtual member functions to existing classes without
+changing their source and recompiling everything. GiNaC comes with source,
+so you could actually do this, but for a small algorithm like the one
+presented this would be impractical.
-@example
+One solution to this dilemma is the @dfn{Visitor} design pattern,
+which is implemented in GiNaC (actually, Robert Martin's Acyclic Visitor
+variation, described in detail in
+@uref{http://objectmentor.com/publications/acv.pdf}). Instead of adding
+virtual functions to the class hierarchy to implement operations, GiNaC
+provides a single "bouncing" method @code{accept()} that takes an instance
+of a special @code{visitor} class and redirects execution to the one
+@code{visit()} virtual function of the visitor that matches the type of
+object that @code{accept()} was being invoked on.
+
+Visitors in GiNaC must derive from the global @code{visitor} class as well
+as from the class @code{T::visitor} of each class @code{T} they want to
+visit, and implement the member functions @code{void visit(const T &)} for
+each class.
+
+A call of
+
+@example
+void ex::accept(visitor & v) const;
+@end example
+
+will then dispatch to the correct @code{visit()} member function of the
+specified visitor @code{v} for the type of GiNaC object at the root of the
+expression tree (e.g. a @code{symbol}, an @code{idx} or a @code{mul}).
+
+Here is an example of a visitor:
+
+@example
+class my_visitor
+ : public visitor, // this is required
+ public add::visitor, // visit add objects
+ public numeric::visitor, // visit numeric objects
+ public basic::visitor // visit basic objects
+@{
+ void visit(const add & x)
+ @{ cout << "called with an add object" << endl; @}
+
+ void visit(const numeric & x)
+ @{ cout << "called with a numeric object" << endl; @}
+
+ void visit(const basic & x)
+ @{ cout << "called with a basic object" << endl; @}
+@};
+@end example
+
+which can be used as follows:
+
+@example
+...
+ symbol x("x");
+ ex e1 = 42;
+ ex e2 = 4*x-3;
+ ex e3 = 8*x;
+
+ my_visitor v;
+ e1.accept(v);
+ // prints "called with a numeric object"
+ e2.accept(v);
+ // prints "called with an add object"
+ e3.accept(v);
+ // prints "called with a basic object"
+...
+@end example
+
+The @code{visit(const basic &)} method gets called for all objects that are
+not @code{numeric} or @code{add} and acts as an (optional) default.
+
+From a conceptual point of view, the @code{visit()} methods of the visitor
+behave like a newly added virtual function of the visited hierarchy.
+In addition, visitors can store state in member variables, and they can
+be extended by deriving a new visitor from an existing one, thus building
+hierarchies of visitors.
+
+We can now rewrite our index example from above with a visitor:
+
+@example
+class gather_indices_visitor
+ : public visitor, public idx::visitor, public varidx::visitor
+@{
+ lst l;
+
+ void visit(const idx & i)
+ @{
+ l.append(i);
+ @}
+
+ void visit(const varidx & vi)
+ @{
+ l.append(vi.is_covariant() ? vi : vi.toggle_variance());
+ @}
+
+public:
+ const lst & get_result() // utility function
+ @{
+ l.sort();
+ l.unique();
+ return l;
+ @}
+@};
+@end example
+
+What's missing is the tree traversal. We could implement it in
+@code{visit(const basic &)}, but GiNaC has predefined methods for this:
+
+@example
+void ex::traverse_preorder(visitor & v) const;
+void ex::traverse_postorder(visitor & v) const;
+void ex::traverse(visitor & v) const;
+@end example
+
+@code{traverse_preorder()} visits a node @emph{before} visiting its
+subexpressions, while @code{traverse_postorder()} visits a node @emph{after}
+visiting its subexpressions. @code{traverse()} is a synonym for
+@code{traverse_preorder()}.
+
+Here is a new implementation of @code{gather_indices()} that uses the visitor
+and @code{traverse()}:
+
+@example
+lst gather_indices(const ex & e)
+@{
+ gather_indices_visitor v;
+ e.traverse(v);
+ return v.get_result();
+@}
+@end example
+
+Alternatively, you could use pre- or postorder iterators for the tree
+traversal:
+
+@example
+lst gather_indices(const ex & e)
+@{
+ gather_indices_visitor v;
+ for (const_preorder_iterator i = e.preorder_begin();
+ i != e.preorder_end(); ++i) @{
+ i->accept(v);
+ @}
+ return v.get_result();
+@}
+@end example
+
+
+@node Polynomial Arithmetic, Rational Expressions, Visitors and Tree Traversal, Methods and Functions
+@c node-name, next, previous, up
+@section Polynomial arithmetic
+
+@subsection Expanding and collecting
+@cindex @code{expand()}
+@cindex @code{collect()}
+@cindex @code{collect_common_factors()}
+
+A polynomial in one or more variables has many equivalent
+representations. Some useful ones serve a specific purpose. Consider
+for example the trivariate polynomial @math{4*x*y + x*z + 20*y^2 +
+21*y*z + 4*z^2} (written down here in output-style). It is equivalent
+to the factorized polynomial @math{(x + 5*y + 4*z)*(4*y + z)}. Other
+representations are the recursive ones where one collects for exponents
+in one of the three variable. Since the factors are themselves
+polynomials in the remaining two variables the procedure can be
+repeated. In our example, two possibilities would be @math{(4*y + z)*x
++ 20*y^2 + 21*y*z + 4*z^2} and @math{20*y^2 + (21*z + 4*x)*y + 4*z^2 +
+x*z}.
+
+To bring an expression into expanded form, its method
+
+@example
+ex ex::expand(unsigned options = 0);
+@end example
+
+may be called. In our example above, this corresponds to @math{4*x*y +
+x*z + 20*y^2 + 21*y*z + 4*z^2}. Again, since the canonical form in
+GiNaC is not easy to guess you should be prepared to see different
+orderings of terms in such sums!
+
+Another useful representation of multivariate polynomials is as a
+univariate polynomial in one of the variables with the coefficients
+being polynomials in the remaining variables. The method
+@code{collect()} accomplishes this task:
+
+@example
+ex ex::collect(const ex & s, bool distributed = false);
+@end example
+
+The first argument to @code{collect()} can also be a list of objects in which
+case the result is either a recursively collected polynomial, or a polynomial
+in a distributed form with terms like @math{c*x1^e1*...*xn^en}, as specified
+by the @code{distributed} flag.
+
+Note that the original polynomial needs to be in expanded form (for the
+variables concerned) in order for @code{collect()} to be able to find the
+coefficients properly.
+
+The following @command{ginsh} transcript shows an application of @code{collect()}
+together with @code{find()}:
+
+@example
> a=expand((sin(x)+sin(y))*(1+p+q)*(1+d));
d*p*sin(x)+p*sin(x)+q*d*sin(x)+q*sin(y)+d*sin(x)+q*d*sin(y)+sin(y)+d*sin(y)+q*sin(x)+d*sin(y)*p+sin(x)+sin(y)*p
> collect(a,@{p,q@});
int ex::ldegree(const ex & s);
@end example
-These functions only work reliably if the input polynomial is collected in
-terms of the object @samp{s}. Otherwise, they are only guaranteed to return
-the upper/lower bounds of the exponents. If you need accurate results, you
-have to call @code{expand()} and/or @code{collect()} on the input polynomial.
-For example
-
-@example
-> a=(x+1)^2-x^2;
-(1+x)^2-x^2;
-> degree(a,x);
-2
-> degree(expand(a),x);
-1
-@end example
-
-@code{degree()} also works on rational functions, returning the asymptotic
-degree:
-
-@example
-> degree((x+1)/(x^3+1),x);
--2
-@end example
-
-If the input is not a polynomial or rational function in the variable @samp{s},
-the behavior of @code{degree()} and @code{ldegree()} is undefined.
-
-To extract a coefficient with a certain power from an expanded
-polynomial you use
+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);
The two functions
@example
-ex quo(const ex & a, const ex & b, const symbol & x);
-ex rem(const ex & a, const ex & b, const symbol & x);
+ex quo(const ex & a, const ex & b, const ex & x);
+ex rem(const ex & a, const ex & b, const ex & x);
@end example
compute the quotient and remainder of univariate polynomials in the variable
The additional function
@example
-ex prem(const ex & a, const ex & b, const symbol & x);
+ex prem(const ex & a, const ex & b, const ex & x);
@end example
computes the pseudo-remainder of @samp{a} and @samp{b} which satisfies
The methods
@example
-ex ex::unit(const symbol & x);
-ex ex::content(const symbol & x);
-ex ex::primpart(const symbol & x);
+ex ex::unit(const ex & x);
+ex ex::content(const ex & x);
+ex ex::primpart(const ex & x);
@end example
return the unit part, content part, and primitive polynomial of a multivariate
faster than using @code{numer()} and @code{denom()} separately.
-@subsection Converting to a rational expression
+@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
above. You do this by calling
@example
-ex ex::to_rational(lst &l);
+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{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_rational()}, so it's possible to use it on multiple expressions
-and get consistent results.
+on the expression to be converted. The supplied @code{exmap} or @code{lst}
+will be filled with the generated temporary symbols and their replacement
+expressions in a format that can be used directly for the @code{subs()}
+method. It can also already contain a list of replacements from an earlier
+application of @code{.to_polynomial()} or @code{.to_rational()}, so it's
+possible to use it on multiple expressions and get consistent results.
+
+The difference 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
-For example,
+The following more useful example will print @samp{sin(x)-cos(x)}:
@example
@{
ex a = pow(sin(x), 2) - pow(cos(x), 2);
ex b = sin(x) + cos(x);
ex q;
- lst l;
- divide(a.to_rational(l), b.to_rational(l), q);
- cout << q.subs(l) << endl;
+ exmap m;
+ divide(a.to_polynomial(m), b.to_polynomial(m), q);
+ cout << q.subs(m) << endl;
@}
@end example
-will print @samp{sin(x)-cos(x)}.
-
@node Symbolic Differentiation, Series Expansion, Rational Expressions, Methods and Functions
@c node-name, next, previous, up
$\pi$
@end tex
(for which there already exists the built-in constant @code{Pi})
-using Machin's amazing formula
+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
-We may expand the arcus tangent around @code{0} and insert the fractions
-@code{1/5} and @code{1/239}. But, 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:
+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>
@}
@end example
-
-@node Built-in Functions, Input/Output, Symmetrization, Methods and Functions
+@node Built-in Functions, Multiple polylogarithms, Symmetrization, Methods and Functions
@c node-name, next, previous, up
@section Predefined mathematical functions
+@c
+@subsection Overview
GiNaC contains the following predefined mathematical functions:
@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))})
@tab natural logarithm
@cindex @code{log()}
@item @code{Li2(x)}
-@tab Dilogarithm
+@tab dilogarithm
@cindex @code{Li2()}
-@item @code{zeta(x)}
-@tab Riemann's zeta function
+@item @code{Li(m, x)}
+@tab classical polylogarithm as well as multiple polylogarithm
+@cindex @code{Li()}
+@item @code{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{zeta(n, x)}
+@item @code{zetaderiv(n, x)}
@tab derivatives of Riemann's zeta function
@item @code{tgamma(x)}
-@tab Gamma function
+@tab gamma function
@cindex @code{tgamma()}
-@cindex Gamma function
+@cindex gamma function
@item @code{lgamma(x)}
-@tab logarithm of Gamma function
+@tab logarithm of gamma function
@cindex @code{lgamma()}
@item @code{beta(x, y)}
-@tab Beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
+@tab beta function (@code{tgamma(x)*tgamma(y)/tgamma(x+y)})
@cindex @code{beta()}
@item @code{psi(x)}
@tab psi (digamma) function
standard incorporate these functions in the complex domain in a manner
compatible with C99.
+@node Multiple polylogarithms, Complex Conjugation, Built-in Functions, Methods and Functions
+@c node-name, next, previous, up
+@subsection Multiple polylogarithms
+
+@cindex polylogarithm
+@cindex Nielsen's generalized polylogarithm
+@cindex harmonic polylogarithm
+@cindex multiple zeta value
+@cindex alternating Euler sum
+@cindex multiple polylogarithm
+
+The multiple polylogarithm is the most generic member of a family of functions,
+to which others like the harmonic polylogarithm, Nielsen's generalized
+polylogarithm and the multiple zeta value belong.
+Everyone of these functions can also be written as a multiple polylogarithm with specific
+parameters. This whole family of functions is therefore often referred to simply as
+multiple polylogarithms, containing @code{Li}, @code{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, Multiple polylogarithms, 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, Built-in Functions, Methods and Functions
+
+@node Input/Output, Extending GiNaC, Solving Linear Systems of Equations, Methods and Functions
@c node-name, next, previous, up
@section Input and output of expressions
@cindex I/O
@cindex printing
@cindex output of expressions
-The easiest way to print an expression is to write it to a stream:
+Expressions can simply be written to any stream:
@example
@{
symbol x("x");
- ex e = 4.5+pow(x,2)*3/2;
- cout << e << endl; // prints '(4.5)+3/2*x^2'
+ ex e = 4.5*I+pow(x,2)*3/2;
+ cout << e << endl; // prints '4.5*I+3/2*x^2'
// ...
@end example
-The output format is identical to the @command{ginsh} input syntax and
+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
-the method
+a set of stream manipulators;
@example
-void ex::print(const print_context & c, unsigned level = 0);
+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
-@cindex @code{print_context} (class)
-The type of @code{print_context} object passed in determines the format
-of the output. The possible types are defined in @file{ginac/print.h}.
-All constructors of @code{print_context} and derived classes take an
-@code{ostream &} as their first argument.
+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.
-To print an expression in a way that can be directly used in a C or C++
-program, you pass a @code{print_csrc} object like this:
+@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 << "float f = ";
- e.print(print_csrc_float(cout));
- cout << ";\n";
+ 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
- cout << "double d = ";
- e.print(print_csrc_double(cout));
- cout << ";\n";
+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:
- cout << "cl_N n = ";
- e.print(print_csrc_cl_N(cout));
- cout << ";\n";
+@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
-The three possible types mostly affect the way in which floating point
-numbers are written.
+@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}):
+The above example will produce (note the @code{x^2} being converted to
+@code{x*x}):
@example
-float f = (3.0/2.0)*(x*x)+4.500000e+00;
-double d = (3.0/2.0)*(x*x)+4.5000000000000000e+00;
-cl_N n = cln::cl_RA("3/2")*(x*x)+cln::cl_F("4.5_17");
+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
-The @code{print_context} type @code{print_tree} provides a dump of the
-internal structure of an expression for debugging purposes:
+@cindex @code{tree}
+The @code{tree} manipulator allows dumping the internal structure of an
+expression for debugging purposes:
@example
// ...
- e.print(print_tree(cout));
+ cout << tree << e;
@}
@end example
@example
add, hash=0x0, flags=0x3, nops=2
- power, hash=0x9, flags=0x3, nops=2
- x (symbol), serial=3, hash=0x44a113a6, flags=0xf
- 2 (numeric), hash=0x80000042, flags=0xf
- 3/2 (numeric), hash=0x80000061, flags=0xf
+ 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.5L0 (numeric), hash=0x8000004b, flags=0xf
+ 4.5L0i (numeric), hash=0xa40a97e0, flags=0xf
=====
@end example
-This kind of output is also available in @command{ginsh} as the @code{print()}
-function.
-
-Another useful output format is for LaTeX parsing in mathematical mode.
-It is rather similar to the default @code{print_context} 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.
+@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");
- ex foo = lgamma(x).series(x==0,3);
- foo.print(print_latex(std::cout));
+@{
+ 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
-will print out:
+@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
- @{(-\ln(x))@}+@{(-\gamma_E)@} x+@{(1/12 \pi^2)@} x^@{2@}+\mathcal@{O@}(x^3)
+@{
+ 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
if (is_a<function>(e))
cout << ex_to<function>(e).get_name();
else
- cout << e.bp->class_name();
+ cout << ex_to<basic>(e).class_name();
cout << "(";
- unsigned n = e.nops();
+ size_t n = e.nops();
if (n)
- for (unsigned i=0; i<n; i++) @{
+ for (size_t i=0; i<n; i++) @{
my_print(e.op(i));
if (i != n-1)
cout << ",";
cout << ")";
@}
-int main(void)
+int main()
@{
my_print(pow(3, x) - 2 * sin(y / Pi)); cout << endl;
return 0;
desired symbols to the @code{>>} stream input operator.
Instead, GiNaC lets you construct an expression from a string, specifying the
-list of symbols and indices to be used:
+list of symbols to be used:
@example
@{
- symbol x("x"), y("y"), p("p");
- idx i(symbol("i"), 3);
- ex e("2*x+sin(y)+p.i", lst(x, y, p, i));
+ 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 and indices in the string are matched
-by name to the symbols and indices in the list and if GiNaC encounters a
-symbol or index not specified in the list it will throw an exception. Only
-indices whose values are single symbols can be used (i.e. numeric indices
-or compound indices as in "A.(2*n+1)" are not allowed).
+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:
// ...
@end example
-And the stored expressions can be retrieved by their name:
+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));
+...
+@end example
+
+Note the difference between @code{sprod} which is the algebraic class, and
+@code{sprod_s} which is the unadorned C++ structure containing the @code{left}
+and @code{right} data members. As shown above, an @code{sprod} can be
+constructed from an @code{sprod_s} object.
+
+If you find the nested @code{sprod(sprod_s())} constructor too unwieldy,
+you could define a little wrapper function like this:
+
+@example
+inline ex make_sprod(ex left, ex right)
+@{
+ return sprod(sprod_s(left, right));
+@}
+@end example
+
+The @code{sprod_s} object contained in @code{sprod} can be accessed with
+the GiNaC @code{ex_to<>()} function followed by the @code{->} operator or
+@code{get_struct()}:
+
+@example
+...
+ cout << ex_to<sprod>(e)->left << endl;
+ // -> a
+ cout << ex_to<sprod>(e).get_struct().right << endl;
+ // -> b
+...
+@end example
+
+You only have read access to the members of @code{sprod_s}.
+
+The type definition of @code{sprod} is enough to write your own algorithms
+that deal with scalar products, for example:
@example
- // ...
- lst syms(x, y);
+ex swap_sprod(ex p)
+@{
+ if (is_a<sprod>(p)) @{
+ const sprod_s & sp = ex_to<sprod>(p).get_struct();
+ return make_sprod(sp.right, sp.left);
+ @} else
+ return p;
+@}
- ex ex1 = a2.unarchive_ex(syms, "foo");
- ex ex2 = a2.unarchive_ex(syms, "the second one");
+...
+ f = swap_sprod(e);
+ // f is now <b|a>
+...
+@end example
- 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"
-@}
+@subsection Structure output
+
+While the @code{sprod} type is useable it still leaves something to be
+desired, most notably proper output:
+
+@example
+...
+ cout << e << endl;
+ // -> [structure object]
+...
@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.
+By default, any structure types you define will be printed as
+@samp{[structure object]}. To override this you can either specialize the
+template's @code{print()} member function, or specify print methods with
+@code{set_print_func<>()}, as described in @ref{Printing}. Unfortunately,
+it's not possible to supply class options like @code{print_func<>()} to
+structures, so for a self-contained structure type you need to resort to
+overriding the @code{print()} function, which is also what we will do here.
-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:
+The member functions of GiNaC classes are described in more detail in the
+next section, but it shouldn't be hard to figure out what's going on here:
@example
-static void my_print2(const archive_node & n)
+void sprod::print(const print_context & c, unsigned level) const
@{
- string class_name;
- n.find_string("class", class_name);
- cout << class_name << "(";
+ // tree debug output handled by superclass
+ if (is_a<print_tree>(c))
+ inherited::print(c, level);
- archive_node::propinfovector p;
- n.get_properties(p);
+ // get the contained sprod_s object
+ const sprod_s & sp = get_struct();
- unsigned num = p.size();
- for (unsigned i=0; i<num; i++) @{
- const string &name = p[i].name;
- if (name == "class")
- continue;
- cout << name << "=";
+ // print_context::s is a reference to an ostream
+ c.s << "<" << sp.left << "|" << sp.right << ">";
+@}
+@end example
- unsigned count = p[i].count;
- if (count > 1)
- cout << "@{";
+Now we can print expressions containing scalar products:
- 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;
- @}
- @}
+@example
+...
+ cout << e << endl;
+ // -> <a|b>
+ cout << swap_sprod(e) << endl;
+ // -> <b|a>
+...
+@end example
- if (j != count-1)
- cout << ",";
- @}
+@subsection Comparing structures
- if (count > 1)
- cout << "@}";
+The @code{sprod} class defined so far still has one important drawback: all
+scalar products are treated as being equal because GiNaC doesn't know how to
+compare objects of type @code{sprod_s}. This can lead to some confusing
+and undesired behavior:
- if (i != num-1)
- cout << ",";
- @}
+@example
+...
+ cout << make_sprod(a, b) - make_sprod(a*a, b*b) << endl;
+ // -> 0
+ cout << make_sprod(a, b) + make_sprod(a*a, b*b) << endl;
+ // -> 2*<a|b> or 2*<a^2|b^2> (which one is undefined)
+...
+@end example
- cout << ")";
+To remedy this, we first need to define the operators @code{==} and @code{<}
+for objects of type @code{sprod_s}:
+
+@example
+inline bool operator==(const sprod_s & lhs, const sprod_s & rhs)
+@{
+ return lhs.left.is_equal(rhs.left) && lhs.right.is_equal(rhs.right);
@}
-int main(void)
+inline bool operator<(const sprod_s & lhs, const sprod_s & rhs)
@{
- ex e = pow(2, x) - y;
- archive ar(e, "e");
- my_print2(ar.get_top_node(0)); cout << endl;
- return 0;
+ return lhs.left.compare(rhs.left) < 0 ? true : lhs.right.compare(rhs.right) < 0;
@}
@end example
-This will produce:
+The ordering established by the @code{<} operator doesn't have to make any
+algebraic sense, but it needs to be well defined. Note that we can't use
+expressions like @code{lhs.left == rhs.left} or @code{lhs.left < rhs.left}
+in the implementation of these operators because they would construct
+GiNaC @code{relational} objects which in the case of @code{<} do not
+establish a well defined ordering (for arbitrary expressions, GiNaC can't
+decide which one is algebraically 'less').
+
+Next, we need to change our definition of the @code{sprod} type to let
+GiNaC know that an ordering relation exists for the embedded objects:
@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"))
+typedef structure<sprod_s, compare_std_less> sprod;
@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.
-* Adding classes:: Defining new algebraic classes.
-@end menu
+@code{sprod} objects then behave as expected:
+@example
+...
+ cout << make_sprod(a, b) - make_sprod(a*a, b*b) << endl;
+ // -> <a|b>-<a^2|b^2>
+ cout << make_sprod(a, b) + make_sprod(a*a, b*b) << endl;
+ // -> <a|b>+<a^2|b^2>
+ cout << make_sprod(a, b) - make_sprod(a, b) << endl;
+ // -> 0
+ cout << make_sprod(a, b) + make_sprod(a, b) << endl;
+ // -> 2*<a|b>
+...
+@end example
-@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.
+The @code{compare_std_less} policy parameter tells GiNaC to use the
+@code{std::less} and @code{std::equal_to} functors to compare objects of
+type @code{sprod_s}. By default, these functors forward their work to the
+standard @code{<} and @code{==} operators, which we have overloaded.
+Alternatively, we could have specialized @code{std::less} and
+@code{std::equal_to} for class @code{sprod_s}.
-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.
+GiNaC provides two other comparison policies for @code{structure<T>}
+objects: the default @code{compare_all_equal}, and @code{compare_bitwise}
+which does a bit-wise comparison of the contained @code{T} objects.
+This should be used with extreme care because it only works reliably with
+built-in integral types, and it also compares any padding (filler bytes of
+undefined value) that the @code{T} class might have.
+@subsection Subexpressions
-@node Symbolic functions, Adding classes, What does not belong into GiNaC, Extending GiNaC
-@c node-name, next, previous, up
-@section Symbolic functions
+Our scalar product class has two subexpressions: the left and right
+operands. It might be a good idea to make them accessible via the standard
+@code{nops()} and @code{op()} methods:
-The easiest and most instructive way to start with is probably to
-implement your own function. GiNaC's functions are objects of class
-@code{function}. The preprocessor is then used to convert the function
-names to objects with a corresponding serial number that is used
-internally to identify them. You usually need not worry about this
-number. New functions may be inserted into the system via a kind of
-`registry'. It is your responsibility to care for some functions that
-are called when the user invokes certain methods. These are usual
-C++-functions accepting a number of @code{ex} as arguments and returning
-one @code{ex}. As an example, if we have a look at a simplified
-implementation of the cosine trigonometric function, we first need a
-function that is called when one wishes to @code{eval} it. It could
-look something like this:
-
-@example
-static ex cos_eval_method(const ex & x)
+@example
+size_t sprod::nops() const
@{
- // if (!x%(2*Pi)) return 1
- // if (!x%Pi) return -1
- // if (!x%Pi/2) return 0
- // care for other cases...
- return cos(x).hold();
+ return 2;
@}
-@end example
-@cindex @code{hold()}
-@cindex evaluation
-The last line returns @code{cos(x)} if we don't know what else to do and
-stops a potential recursive evaluation by saying @code{.hold()}, which
-sets a flag to the expression signaling that it has been evaluated. We
-should also implement a method for numerical evaluation and since we are
-lazy we sweep the problem under the rug by calling someone else's
-function that does so, in this case the one in class @code{numeric}:
-
-@example
-static ex cos_evalf(const ex & x)
+ex sprod::op(size_t i) const
@{
- if (is_a<numeric>(x))
- return cos(ex_to<numeric>(x));
- else
- return cos(x).hold();
+ switch (i) @{
+ case 0:
+ return get_struct().left;
+ case 1:
+ return get_struct().right;
+ default:
+ throw std::range_error("sprod::op(): no such operand");
+ @}
@}
@end example
-Differentiation will surely turn up and so we need to tell @code{cos}
-what the first derivative is (higher derivatives (@code{.diff(x,3)} for
-instance are then handled automatically by @code{basic::diff} and
-@code{ex::diff}):
+Implementing @code{nops()} and @code{op()} for container types such as
+@code{sprod} has two other nice side effects:
+
+@itemize @bullet
+@item
+@code{has()} works as expected
+@item
+GiNaC generates better hash keys for the objects (the default implementation
+of @code{calchash()} takes subexpressions into account)
+@end itemize
+
+@cindex @code{let_op()}
+There is a non-const variant of @code{op()} called @code{let_op()} that
+allows replacing subexpressions:
@example
-static ex cos_deriv(const ex & x, unsigned diff_param)
+ex & sprod::let_op(size_t i)
@{
- return -sin(x);
+ // every non-const member function must call this
+ ensure_if_modifiable();
+
+ switch (i) @{
+ case 0:
+ return get_struct().left;
+ case 1:
+ return get_struct().right;
+ default:
+ throw std::range_error("sprod::let_op(): no such operand");
+ @}
@}
@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. For Taylor expansion, it is
-enough to know how to differentiate. But if the function you want to
-implement does have a pole somewhere in the complex plane, you need to
-write another method for Laurent expansion around that point.
+Once we have provided @code{let_op()} we also get @code{subs()} and
+@code{map()} for free. In fact, every container class that returns a non-null
+@code{nops()} value must either implement @code{let_op()} or provide custom
+implementations of @code{subs()} and @code{map()}.
-Now that all the ingredients for @code{cos} have been set up, we need
-to tell the system about it. This is done by a macro and we are not
-going to describe how it expands, please consult your preprocessor if you
-are curious:
+In turn, the availability of @code{map()} enables the recursive behavior of a
+couple of other default method implementations, in particular @code{evalf()},
+@code{evalm()}, @code{normal()}, @code{diff()} and @code{expand()}. Although
+we probably want to provide our own version of @code{expand()} for scalar
+products that turns expressions like @samp{<a+b|c>} into @samp{<a|c>+<b|c>}.
+This is left as an exercise for the reader.
-@example
-REGISTER_FUNCTION(cos, eval_func(cos_eval).
- evalf_func(cos_evalf).
- derivative_func(cos_deriv));
-@end example
-
-The first argument is the function's name used for calling it and for
-output. The second binds the corresponding methods as options to this
-object. Options are separated by a dot and can be given in an arbitrary
-order. GiNaC functions understand several more options which are always
-specified as @code{.option(params)}, for example a method for series
-expansion @code{.series_func(cos_series)}. Again, if no series
-expansion method is given, GiNaC defaults to simple Taylor expansion,
-which is correct if there are no poles involved as is the case for the
-@code{cos} function. The way GiNaC handles poles in case there are any
-is best understood by studying one of the examples, like the Gamma
-(@code{tgamma}) function for instance. (In essence the function first
-checks if there is a pole at the evaluation point and falls back to
-Taylor expansion if there isn't. Then, the pole is regularized by some
-suitable transformation.) Also, the new function needs to be declared
-somewhere. This may also be done by a convenient preprocessor macro:
+The @code{structure<T>} template defines many more member functions that
+you can override by specialization to customize the behavior of your
+structures. You are referred to the next section for a description of
+some of these (especially @code{eval()}). There is, however, one topic
+that shall be addressed here, as it demonstrates one peculiarity of the
+@code{structure<T>} template: archiving.
+
+@subsection Archiving structures
+
+If you don't know how the archiving of GiNaC objects is implemented, you
+should first read the next section and then come back here. You're back?
+Good.
+
+To implement archiving for structures it is not enough to provide
+specializations for the @code{archive()} member function and the
+unarchiving constructor (the @code{unarchive()} function has a default
+implementation). You also need to provide a unique name (as a string literal)
+for each structure type you define. This is because in GiNaC archives,
+the class of an object is stored as a string, the class name.
+
+By default, this class name (as returned by the @code{class_name()} member
+function) is @samp{structure} for all structure classes. This works as long
+as you have only defined one structure type, but if you use two or more you
+need to provide a different name for each by specializing the
+@code{get_class_name()} member function. Here is a sample implementation
+for enabling archiving of the scalar product type defined above:
@example
-DECLARE_FUNCTION_1P(cos)
+const char *sprod::get_class_name() @{ return "sprod"; @}
+
+void sprod::archive(archive_node & n) const
+@{
+ inherited::archive(n);
+ n.add_ex("left", get_struct().left);
+ n.add_ex("right", get_struct().right);
+@}
+
+sprod::structure(const archive_node & n, lst & sym_lst) : inherited(n, sym_lst)
+@{
+ n.find_ex("left", get_struct().left, sym_lst);
+ n.find_ex("right", get_struct().right, sym_lst);
+@}
@end example
-The suffix @code{_1P} stands for @emph{one parameter}. Of course, this
-implementation of @code{cos} is very incomplete and lacks several safety
-mechanisms. Please, have a look at the real implementation in GiNaC.
-(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.)
+Note that the unarchiving constructor is @code{sprod::structure} and not
+@code{sprod::sprod}, and that we don't need to supply an
+@code{sprod::unarchive()} function.
-@node Adding classes, A Comparison With Other CAS, Symbolic functions, Extending GiNaC
+@node Adding classes, A Comparison With Other CAS, Structures, Extending GiNaC
@c node-name, next, previous, up
@section Adding classes
-If you are doing some very specialized things with GiNaC you may find that
-you have to implement your own algebraic classes to fit your needs. This
-section will explain how to do this by giving the example of a simple
-'string' class. After reading this section you will know how to properly
-declare a GiNaC class and what the minimum required member functions are
-that you have to implement. We only cover the implementation of a 'leaf'
-class here (i.e. one that doesn't contain subexpressions). Creating a
-container class like, for example, a class representing tensor products is
-more involved but this section should give you enough information so you can
-consult the source to GiNaC's predefined classes if you want to implement
-something more complicated.
+The @code{structure<T>} template provides an way to extend GiNaC with custom
+algebraic classes that is easy to use but has its limitations, the most
+severe of which being that you can't add any new member functions to
+structures. To be able to do this, you need to write a new class definition
+from scratch.
+
+This section will explain how to implement new algebraic classes in GiNaC by
+giving the example of a simple 'string' class. After reading this section
+you will know how to properly declare a GiNaC class and what the minimum
+required member functions are that you have to implement. We only cover the
+implementation of a 'leaf' class here (i.e. one that doesn't contain
+subexpressions). Creating a container class like, for example, a class
+representing tensor products is more involved but this section should give
+you enough information so you can consult the source to GiNaC's predefined
+classes if you want to implement something more complicated.
@subsection GiNaC's run-time type information system
The disadvantage of this proprietary RTTI implementation is that there's
a little more to do when implementing new classes (C++'s RTTI works more
-or less automatic) but don't worry, most of the work is simplified by
+or less automatically) but don't worry, most of the work is simplified by
macros.
@subsection A minimalistic example
@end example
The @code{GINAC_DECLARE_REGISTERED_CLASS} and @code{GINAC_IMPLEMENT_REGISTERED_CLASS}
-macros are defined in @file{registrar.h}. They take the name of the class
+macros are defined in @file{registrar.h}. They take the name of the class
and its direct superclass as arguments and insert all required declarations
for the RTTI system. The @code{GINAC_DECLARE_REGISTERED_CLASS} should be
the first line after the opening brace of the class definition. The
source (at global scope, of course, not inside a function).
@code{GINAC_DECLARE_REGISTERED_CLASS} contains, among other things the
-declarations of the default and copy constructor, the destructor, the
-assignment operator and a couple of other functions that are required. It
-also defines a type @code{inherited} which refers to the superclass so you
-don't have to modify your code every time you shuffle around the class
-hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} implements the copy
-constructor, the destructor and the assignment operator.
-
-Now there are nine member functions we have to implement to get a working
+declarations of the default constructor and a couple of other functions that
+are required. It also defines a type @code{inherited} which refers to the
+superclass so you don't have to modify your code every time you shuffle around
+the class hierarchy. @code{GINAC_IMPLEMENT_REGISTERED_CLASS} registers the
+class with the GiNaC RTTI (there is also a
+@code{GINAC_IMPLEMENT_REGISTERED_CLASS_OPT} which allows specifying additional
+options for the class, and which we will be using instead in a few minutes).
+
+Now there are seven member functions we have to implement to get a working
class:
@itemize
@item
@code{mystring()}, the default constructor.
-@item
-@code{void destroy(bool call_parent)}, which is used in the destructor and the
-assignment operator to free dynamically allocated members. The @code{call_parent}
-specifies whether the @code{destroy()} function of the superclass is to be
-called also.
-
-@item
-@code{void copy(const mystring &other)}, which is used in the copy constructor
-and assignment operator to copy the member variables over from another
-object of the same class.
-
@item
@code{void archive(archive_node &n)}, the archiving function. This stores all
information needed to reconstruct an object of this class inside an
@code{archive_node}.
@item
-@code{mystring(const archive_node &n, const lst &sym_lst)}, the unarchiving
+@code{mystring(const archive_node &n, lst &sym_lst)}, the unarchiving
constructor. This constructs an instance of the class from the information
found in an @code{archive_node}.
@item
-@code{ex unarchive(const archive_node &n, const lst &sym_lst)}, the static
+@code{ex unarchive(const archive_node &n, lst &sym_lst)}, the static
unarchiving function. It constructs a new instance by calling the unarchiving
constructor.
@item
+@cindex @code{compare_same_type()}
@code{int compare_same_type(const basic &other)}, which is used internally
by GiNaC to establish a canonical sort order for terms. It returns 0, +1 or
-1, depending on the relative order of this object and the @code{other}
Let's proceed step-by-step. The default constructor looks like this:
@example
-mystring::mystring() : inherited(TINFO_mystring)
-@{
- // dynamically allocate resources here if required
-@}
+mystring::mystring() : inherited(TINFO_mystring) @{@}
@end example
The golden rule is that in all constructors you have to set the
it will be set by the constructor of the superclass and all hell will break
loose in the RTTI. For your convenience, the @code{basic} class provides
a constructor that takes a @code{tinfo_key} value, which we are using here
-(remember that in our case @code{inherited = basic}). If the superclass
+(remember that in our case @code{inherited == basic}). If the superclass
didn't have such a constructor, we would have to set the @code{tinfo_key}
to the right value manually.
In the default constructor you should set all other member variables to
reasonable default values (we don't need that here since our @code{str}
-member gets set to an empty string automatically). The constructor(s) are of
-course also the right place to allocate any dynamic resources you require.
-
-Next, the @code{destroy()} function:
-
-@example
-void mystring::destroy(bool call_parent)
-@{
- // free dynamically allocated resources here if required
- if (call_parent)
- inherited::destroy(call_parent);
-@}
-@end example
-
-This function is where we free all dynamically allocated resources. We
-don't have any so we're not doing anything here, but if we had, for
-example, used a C-style @code{char *} to store our string, this would be
-the place to @code{delete[]} the string storage. If @code{call_parent}
-is true, we have to call the @code{destroy()} function of the superclass
-after we're done (to mimic C++'s automatic invocation of superclass
-destructors where @code{destroy()} is called from outside a destructor).
-
-The @code{copy()} function just copies over the member variables from
-another object:
-
-@example
-void mystring::copy(const mystring &other)
-@{
- inherited::copy(other);
- str = other.str;
-@}
-@end example
-
-We can simply overwrite the member variables here. There's no need to worry
-about dynamically allocated storage. The assignment operator (which is
-automatically defined by @code{GINAC_IMPLEMENT_REGISTERED_CLASS}, as you
-recall) calls @code{destroy()} before it calls @code{copy()}. You have to
-explicitly call the @code{copy()} function of the superclass here so
-all the member variables will get copied.
+member gets set to an empty string automatically).
Next are the three functions for archiving. You have to implement them even
if you don't plan to use archives, but the minimum required implementation
function:
@example
-mystring::mystring(const archive_node &n, const lst &sym_lst) : inherited(n, sym_lst)
+mystring::mystring(const archive_node &n, lst &sym_lst) : inherited(n, sym_lst)
@{
n.find_string("string", str);
@}
Finally, the unarchiving function:
@example
-ex mystring::unarchive(const archive_node &n, const lst &sym_lst)
+ex mystring::unarchive(const archive_node &n, lst &sym_lst)
@{
return (new mystring(n, sym_lst))->setflag(status_flags::dynallocated);
@}
Now the only thing missing is our two new constructors:
@example
-mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s)
-@{
- // dynamically allocate resources here if required
-@}
-
-mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s)
-@{
- // dynamically allocate resources here if required
-@}
+mystring::mystring(const string &s) : inherited(TINFO_mystring), str(s) @{@}
+mystring::mystring(const char *s) : inherited(TINFO_mystring), str(s) @{@}
@end example
No surprises here. We set the @code{str} member from the argument and
@end example
Hm, not exactly what we expect, but of course the @code{mystring} class
-doesn't yet know how to print itself. This is done in the @code{print()}
-member function. Let's say that we wanted to print the string surrounded
-by double quotes:
+doesn't yet know how to print itself. This can be done either by implementing
+the @code{print()} member function, or, preferably, by specifying a
+@code{print_func<>()} class option. Let's say that we want to print the string
+surrounded by double quotes:
@example
class mystring : public basic
@{
...
-public:
- void print(const print_context &c, unsigned level = 0) const;
+protected:
+ void do_print(const print_context &c, unsigned level = 0) const;
...
@};
-void mystring::print(const print_context &c, unsigned level) const
+void mystring::do_print(const print_context &c, unsigned level) const
@{
// print_context::s is a reference to an ostream
c.s << '\"' << str << '\"';
@end example
The @code{level} argument is only required for container classes to
-correctly parenthesize the output. Let's try again to print the expression:
+correctly parenthesize the output.
+
+Now we need to tell GiNaC that @code{mystring} objects should use the
+@code{do_print()} member function for printing themselves. For this, we
+replace the line
+
+@example
+GINAC_IMPLEMENT_REGISTERED_CLASS(mystring, basic)
+@end example
+
+with
+
+@example
+GINAC_IMPLEMENT_REGISTERED_CLASS_OPT(mystring, basic,
+ print_func<print_context>(&mystring::do_print))
+@end example
+
+Let's try again to print the expression:
@example
cout << e << endl;
// -> "Hello, world!"
@end example
-Much better. The @code{mystring} class can be used in arbitrary expressions:
+Much better. If we wanted to have @code{mystring} objects displayed in a
+different way depending on the output format (default, LaTeX, etc.), we
+would have supplied multiple @code{print_func<>()} options with different
+template parameters (@code{print_dflt}, @code{print_latex}, etc.),
+separated by dots. This is similar to the way options are specified for
+symbolic functions. @xref{Printing}, for a more in-depth description of the
+way expression output is implemented in GiNaC.
+
+The @code{mystring} class can be used in arbitrary expressions:
@example
e += mystring("GiNaC rulez");
// -> 3*"wow"
@end example
-@subsection Other member functions
+@subsection Optional member functions
We have implemented only a small set of member functions to make the class
-work in the GiNaC framework. For a real algebraic class, there are probably
-some more functions that you will want to re-implement, such as
-@code{evalf()}, @code{series()} or @code{op()}. Have a look at @file{basic.h}
-or the header file of the class you want to make a subclass of to see
-what's there. One member function that you will most likely want to
-implement for terminal classes like the described string class is
-@code{calcchash()} that returns an @code{unsigned} hash value for the object
-which will allow GiNaC to compare and canonicalize expressions much more
-efficiently.
-
-You can, of course, also add your own new member functions. Remember,
+work in the GiNaC framework. There are two functions that are not strictly
+required but will make operations with objects of the class more efficient:
+
+@cindex @code{calchash()}
+@cindex @code{is_equal_same_type()}
+@example
+unsigned calchash() const;
+bool is_equal_same_type(const basic &other) const;
+@end example
+
+The @code{calchash()} method returns an @code{unsigned} hash value for the
+object which will allow GiNaC to compare and canonicalize expressions much
+more efficiently. You should consult the implementation of some of the built-in
+GiNaC classes for examples of hash functions. The default implementation of
+@code{calchash()} calculates a hash value out of the @code{tinfo_key} of the
+class and all subexpressions that are accessible via @code{op()}.
+
+@code{is_equal_same_type()} works like @code{compare_same_type()} but only
+tests for equality without establishing an ordering relation, which is often
+faster. The default implementation of @code{is_equal_same_type()} just calls
+@code{compare_same_type()} and tests its result for zero.
+
+@subsection Other member functions
+
+For a real algebraic class, there are probably some more functions that you
+might want to provide:
+
+@example
+bool info(unsigned inf) const;
+ex evalf(int level = 0) const;
+ex series(const relational & r, int order, unsigned options = 0) const;
+ex derivative(const symbol & s) const;
+@end example
+
+If your class stores sub-expressions (see the scalar product example in the
+previous section) you will probably want to override
+
+@cindex @code{let_op()}
+@example
+size_t nops() cont;
+ex op(size_t i) const;
+ex & let_op(size_t i);
+ex subs(const lst & ls, const lst & lr, unsigned options = 0) const;
+ex map(map_function & f) const;
+@end example
+
+@code{let_op()} is a variant of @code{op()} that allows write access. The
+default implementations of @code{subs()} and @code{map()} use it, so you have
+to implement either @code{let_op()}, or @code{subs()} and @code{map()}.
+
+You can, of course, also add your own new member functions. Remember
that the RTTI may be used to get information about what kinds of objects
you are dealing with (the position in the class hierarchy) and that you
can always extract the bare object from an @code{ex} by stripping the
@cindex reference counting
@cindex copy-on-write
@cindex garbage collection
-An expression is extremely light-weight since internally it works like a
-handle to the actual representation and really holds nothing more than a
-pointer to some other object. What this means in practice is that
-whenever you create two @code{ex} and set the second equal to the first
-no copying process is involved. Instead, the copying takes place as soon
-as you try to change the second. Consider the simple sequence of code:
+In GiNaC, there is an @emph{intrusive reference-counting} mechanism at work
+where the counter belongs to the algebraic objects derived from class
+@code{basic} but is maintained by the smart pointer class @code{ptr}, of
+which @code{ex} contains an instance. If you understood that, you can safely
+skip the rest of this passage.
+
+Expressions are extremely light-weight since internally they work like
+handles to the actual representation. They really hold nothing more
+than a pointer to some other object. What this means in practice is
+that whenever you create two @code{ex} and set the second equal to the
+first no copying process is involved. Instead, the copying takes place
+as soon as you try to change the second. Consider the simple sequence
+of code:
@example
#include <iostream>
happens. Knowing this will enable you to write much more efficient
code. If you still have an uncertain feeling with copy-on-write
semantics, we recommend you have a look at the
-@uref{http://www.cerfnet.com/~mpcline/c++-faq-lite/, C++-FAQ lite} by
+@uref{http://www.parashift.com/c++-faq-lite/, C++-FAQ lite} by
Marshall Cline. Chapter 16 covers this issue and presents an
implementation which is pretty close to the one in GiNaC.
and the @samp{AM_PATH_GINAC} macro. The program used here is @file{simple.cpp}:
@example
+#include <iostream>
#include <ginac/ginac.h>
-int main(void)
+int main()
@{
GiNaC::symbol x("x");
GiNaC::ex a = GiNaC::sin(x);
@item
@cite{Computer Algebra: Systems and Algorithms for Algebraic Computation},
-James H. Davenport, Yvon Siret, and Evelyne Tournier, ISBN 0-12-204230-1, 1988,
+James H. Davenport, Yvon Siret and Evelyne Tournier, ISBN 0-12-204230-1, 1988,
Academic Press, London
@item
@cite{The Art of Computer Programming, Vol 2: Seminumerical Algorithms},
Donald E. Knuth, ISBN 0-201-89684-2, 1998, Addison Wesley
+@item
+@cite{Pi Unleashed}, J@"org Arndt and Christoph Haenel,
+ISBN 3-540-66572-2, 2001, Springer, Heidelberg
+
@item
@cite{The Role of gamma5 in Dimensional Regularization}, Dirk Kreimer, hep-ph/9401354
git://www.ginac.de / ginac.git / blobdiff