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
basic that can be used in GiNaC expressions).
GiNaC offers two ways of accomplishing this: either by using the
structure<T> template class, or by rolling your own class from
scratch. This section will discuss the 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 structure<T> template, defined in the GiNaC header file
structure.h, wraps a type that you supply (usually a C++ struct
or class) into a GiNaC object that can be used in expressions.
Let's suppose that we need a way to handle some kind of abstract scalar
product of the form ‘<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++ struct:
#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) {}
};
The default constructor is required. Now, to make a GiNaC class out of this data structure, we need only one line:
typedef structure<sprod_s> sprod;
That's it. This line constructs an algebraic class sprod which
contains objects of type sprod_s. We can now use sprod in
expressions like any other GiNaC class:
...
symbol a("a"), b("b");
ex e = sprod(sprod_s(a, b));
...
Note the difference between sprod which is the algebraic class, and
sprod_s which is the unadorned C++ structure containing the left
and right data members. As shown above, an sprod can be
constructed from an sprod_s object.
If you find the nested sprod(sprod_s()) constructor too unwieldy,
you could define a little wrapper function like this:
inline ex make_sprod(ex left, ex right)
{
return sprod(sprod_s(left, right));
}
The sprod_s object contained in sprod can be accessed with
the GiNaC ex_to<>() function followed by the -> operator or
get_struct():
...
cout << ex_to<sprod>(e)->left << endl;
// -> a
cout << ex_to<sprod>(e).get_struct().right << endl;
// -> b
...
You only have read access to the members of sprod_s.
The type definition of sprod is enough to write your own algorithms
that deal with scalar products, for example:
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;
}
...
f = swap_sprod(e);
// f is now <b|a>
...
While the sprod type is useable it still leaves something to be
desired, most notably proper output:
...
cout << e << endl;
// -> [structure object]
...
By default, any structure types you define will be printed as
‘[structure object]’. To override this you can either specialize the
template's print() member function, or specify print methods with
set_print_func<>(), as described in Printing. Unfortunately,
it's not possible to supply class options like print_func<>() to
structures, so for a self-contained structure type you need to resort to
overriding the print() function, which is also what we will do here.
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:
void sprod::print(const print_context & c, unsigned level) const
{
// tree debug output handled by superclass
if (is_a<print_tree>(c))
inherited::print(c, level);
// get the contained sprod_s object
const sprod_s & sp = get_struct();
// print_context::s is a reference to an ostream
c.s << "<" << sp.left << "|" << sp.right << ">";
}
Now we can print expressions containing scalar products:
...
cout << e << endl;
// -> <a|b>
cout << swap_sprod(e) << endl;
// -> <b|a>
...
The 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 sprod_s. This can lead to some confusing
and undesired behavior:
...
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)
...
To remedy this, we first need to define the operators == and <
for objects of type sprod_s:
inline bool operator==(const sprod_s & lhs, const sprod_s & rhs)
{
return lhs.left.is_equal(rhs.left) && lhs.right.is_equal(rhs.right);
}
inline bool operator<(const sprod_s & lhs, const sprod_s & rhs)
{
return lhs.left.compare(rhs.left) < 0
? true : lhs.right.compare(rhs.right) < 0;
}
The ordering established by the < operator doesn't have to make any
algebraic sense, but it needs to be well defined. Note that we can't use
expressions like lhs.left == rhs.left or lhs.left < rhs.left
in the implementation of these operators because they would construct
GiNaC relational objects which in the case of < 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 sprod type to let
GiNaC know that an ordering relation exists for the embedded objects:
typedef structure<sprod_s, compare_std_less> sprod;
sprod objects then behave as expected:
...
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>
...
The compare_std_less policy parameter tells GiNaC to use the
std::less and std::equal_to functors to compare objects of
type sprod_s. By default, these functors forward their work to the
standard < and == operators, which we have overloaded.
Alternatively, we could have specialized std::less and
std::equal_to for class sprod_s.
GiNaC provides two other comparison policies for structure<T>
objects: the default compare_all_equal, and compare_bitwise
which does a bit-wise comparison of the contained 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 T class might have.
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
nops() and op() methods:
size_t sprod::nops() const
{
return 2;
}
ex sprod::op(size_t i) const
{
switch (i) {
case 0:
return get_struct().left;
case 1:
return get_struct().right;
default:
throw std::range_error("sprod::op(): no such operand");
}
}
Implementing nops() and op() for container types such as
sprod has two other nice side effects:
has() works as expected
calchash() takes subexpressions into account)
There is a non-const variant of op() called let_op() that
allows replacing subexpressions:
ex & sprod::let_op(size_t i)
{
// 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");
}
}
Once we have provided let_op() we also get subs() and
map() for free. In fact, every container class that returns a non-null
nops() value must either implement let_op() or provide custom
implementations of subs() and map().
In turn, the availability of map() enables the recursive behavior of a
couple of other default method implementations, in particular evalf(),
evalm(), normal(), diff() and expand(). Although
we probably want to provide our own version of expand() for scalar
products that turns expressions like ‘<a+b|c>’ into ‘<a|c>+<b|c>’.
This is left as an exercise for the reader.
The 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 eval()). There is, however, one topic
that shall be addressed here, as it demonstrates one peculiarity of the
structure<T> template: archiving.
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 archive() member function and the
unarchiving constructor (the 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 class_name() member
function) is ‘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
get_class_name() member function. Here is a sample implementation
for enabling archiving of the scalar product type defined above:
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);
}
Note that the unarchiving constructor is sprod::structure and not
sprod::sprod, and that we don't need to supply an
sprod::unarchive() function.